Skewed maxent

getdist[n_,var_,skew_]:=Module[{},
h[x_]=-Piecewise[{{x Log[x],x>0}},0];
SetAttributes[h,Listable];
ent[x_]:=Total[h[x]];
vals=Table[i,{i,1/n,1,1/n}];
vals=Reverse[-vals]~Join~vals;
cons={vals^0,vals^1,vals^2,vals^3}//N;
ns=Orthogonalize[NullSpace[cons]]//Transpose;
vars=Table[(Subscript[x, i]),{i,1,Dimensions[ns]//Last}];
sol=LinearSolve[cons,{1,0,var,skew}];
realsol=LinearProgramming[ConstantArray[1,Length[vars]],ns,-sol,-Infinity]//Chop;
distr=(ns.vars+sol);
poscons=And@@(#>=0&/@distr);
obj=ent[distr];
distsol=distr/.(FindMaximum[{obj,poscons},Thread[{vars,realsol}]]//Last)
];
dd=getdist[10,.4,.23];
ListPlot[dd,Filling->Axis,PlotRange->{{-1.1,1.1},{0,.3}},DataRange->{-1,1},PlotStyle->PointSize[Medium]]

Space of solutions for 6 bin problem

(* enable subscripts *)
<<Notation`;
Symbolize[(Subscript[_, _])];
genplot[var_,skew_]:=Module[{},
ent[x_]:=-Total[x Log[2,x]];
n=3;
vals=Table[i,{i,1/n,1,1/n}];
vals=Reverse[-vals]~Join~vals;
cons={vals^0,vals^1,vals^2,vals^3};
ns=Orthogonalize[NullSpace[cons]]//Transpose;
vars={(Subscript[x, 1]),(Subscript[x, 2])};
sol=LinearSolve[cons,{1,0,var,skew}];
distr=(ns.vars+sol);
poscons=And@@(#>0&/@distr);
obj=ent[distr];
{maxx,maxy}={(Subscript[x, 1]),(Subscript[x, 2])}/.(NMaximize[{If[poscons,obj,0],poscons},vars]//Last);
gs[c_]:=LinearProgramming[c,ns,-sol,-Infinity];
{{x1,x2},{y1,y2}}={{gs[{1,0}],gs[{-1,0}]}[[All,1]],{gs[{0,1}],gs[{0,-1}]}[[All,2]]};
plot=ContourPlot[obj,{(Subscript[x, 1]),x1,x2},{(Subscript[x, 2]),y1,y2},RegionFunction->Function[{(Subscript[x, 1]),(Subscript[x, 2])},poscons],Contours->5,ContourLabels->Automatic];
Show[plot,ListPlot[{{maxx,maxy}},PlotStyle->PointSize[Large]],Frame->None]
];
genplot[.3,0]
genplot[.3,.15]