The convexity of beta density

In this document, we prove that (Subscript[f, 2]^(-1))(f) is convex, where f(α,β,x) is the beta density function, (Subscript[f, 2]^(-1)) is the inverse of the right branch, assuming that β>α>1.

Beta density

f[α_,β_,x_]=PDF[BetaDistribution[α,β],x]

We only need f for x ∈(0,1).

f0[α_,β_,x_]=Evaluate[f[α,β,x]//Simplify[#,0<x<1]&]

df=D[f0[α,β,x],x]

Let’s find m.

df*Beta[α,β]==0//Simplify

msol=Reduce[%&&0<x<1&&β>α>1,x,Reals]

This is the peak of the density.

m[α_,β_]=Evaluate[msol[[3,-1]]]

An example

We try an example. This can be skipped.

a0=3;b0=7;

f1[x_]:=f0[a0,b0,x];

m1=m[a0,b0]

y1=f1[m1]

Plot[f1[x],{x,0,1},GridLines->{{m1}, {y1}},PlotLabel->"Density function of Beta[3,7]"]

Looks right. Let’s find the inverse on both side.

sol1=Reduce[f1[x]==y&&0<x<1,x,Reals]

We get the two inverse, (Subscript[f, 1]^(-1)) for the left, (Subscript[f, 2]^(-1)) for the right.

fi1[y_]=sol1[[1,2,1,-1]]

fi2[y_]=sol1[[1,2,2,-1]]

Plot[{fi1[y],fi2[y]},{y,0,y1},PlotLegends->{"(Subscript[f, 1]^(-1))[x]","(Subscript[f, 2]^(-1))[x]"},AspectRatio->1,PlotLabel->"The inverse of the density function of Beta[3,7]"]

So we want this function to be convex on (0,(1/4)).

r1[x_]:=Piecewise[{{fi1[f1[x]],1≥x≥m1},{fi2[f1[x]],0≤x<m1}}];

Plot[r1[x]//N[#,30]&,{x,0,1},PlotRange->All,AspectRatio->1,GridLines->{{m1}, {m1}}]

It’s obviously convex.

Let’s prove this symbolically

The derivative

ff[x_]:=f0[α,β,x]

The symbolic inverse is

gg=InverseFunction[ff]

Take 2nd derivative.

dr[0]=D[gg[ff[x]],{x,2}];

Proof for α=3 and β=7

Can be skipped. Let’s just do the simple one.

dr[1]=dr[0]/.{α->3,β->7}

In this case we have the inverse functions, so we just plug them in.

dr[3][x_]=Piecewise[{{dr[1]/.InverseFunction[ff]->fi1,1≥x≥m1},{dr[1]/.InverseFunction[ff]->fi2,0≤x<m1}}];

Let’s first look at the left side.

dr[4]=dr[3][x]//Simplify[#,1≥x≥1/4]&

So the 2nd derivative is positive.

Reduce[dr[4]>0&&1>x>1/4,x,Reals]

Let’s look at the right side.

dr[5]=dr[3][x]//Simplify[#,1/4>x>0]&

Reduce[dr[5]>0&&1/4>x>0,x,Reals]

So we are done.

Module[{p,l,l1,pl,labels},
pl={{1/8,0},{1/8,f1[1/8]},{r1[1/8],f1[1/8]},{r1[1/8],0}};
p=Plot[f1[x],{x,0,1},PlotLabel->"r[x] reflects x through the density function of Beta[3,7]"];
l = Graphics[{Arrow[pl]}];
l1=ListPlot[pl,PlotStyle->Directive[Red,PointSize[.02]]];
labels={Text["x",pl[[1]]+{0.03,0.15}],Text["r[x]",pl[[4]]+{0.05,0.15}]};
Show[p,Graphics[{labels}],l,l1]]

Genera Proof

The first phase.

We reduce the problem to show that x+x1≥2 m[α,β], where x1 is the reflection of x.

We want the 2nd derivative to be positive when α<β. First simplify it a bit.

dd[0]=dr[0]//Simplify

Break it to terms and remove these that are positive. Note that (ff^((-1)))[(((((1-x))^(-1+β)) (x^(-1+α)))/(Beta[α,β]))] is simply the reflection point of x, we replace it by x1 to make things look better.

dd[1]=dd[0]/.(ff^((-1)))[(((((1-x))^(-1+β)) (x^(-1+α)))/(Beta[α,β]))]->x1;

Break it to list of terms

List @@ dd[1]

Get rid of some obviously positive terms (the first 4)

dd[3]=%[[5;;]]

The first term of the left over is also positive.

assump=0<x<m[α,β]<x1<1&&β>α>1;

Reduce[dd[3][[1]]>0&&assump,x1,Reals]

Now only one left. We want it to be positive.

dd[4]=((dd[3][[2]])/((1-x) x))//Simplify

This is equivalent to say

ieq1=dd[4][[-1]]>-dd[4][[1]]//Simplify

Note that the following is true since x and x1 are reflection points.

eq1=(((1-x))^(-1+β)) (x^(-1+α))==(((1-x1))^(-1+β)) (x1^(-1+α));

So ieq1 simplifies to

ieq2=ieq1[[1,3;;]]≥ieq1[[2,3;;]]//FullSimplify[#,assump]&

ieq3=Map[FullSimplify[((#)/(-2+α+β))]&,ieq2]//FullSimplify

Note the right hand side is just 2 m.

ieq3[[2]]==2m[α,β]

Show that x+x1≥2m

When x=0

ieq3/.{x->0,x1->1}//Simplify[#,assump]&

When x=m

ieq3/.{x->m[α,β],x1->m[α,β]}//Simplify

If we look at the picture, it seems obvious that ieq3 is true.

Plot[f1[x],{x,0,1},GridLines->{{m1}, {y1}}]

Note that

x+(Subscript[x, 1])≥2m => x + (Subscript[f, 2]^(-1))[f[x]]≥2m=>(Subscript[f, 2]^(-1))[f[x]]≥2m-x

Since (Subscript[f, 2]^(-1)) is decreasing, this is equivalent to

f[x]≥f[2 m-x]

Let

h1=f0[α,β,2 m[α,β] - x]

h2=f0[α,β,x]

We want the following to be positive

h=h1-h2

To get rid of the powers, we take log of the two terms.

lh=Log[h1]-Log[h2]

If this remains positive, then we are good.

lh1=PowerExpand[lh]

A picture helps.

lh1/.α->3/.β->7

Plot[%,{x,0,m1},PlotRange->Full]

Take derivative.

dlh=D[lh1,x]//Simplify

Reduce[dlh<0&&assump,x,Reals]

So the derivative is always <0.

And

Limit[lh1,x->0,Assumptions->assump]

lh1/.x->m[α,β]

So there can be only one zero of lh1 on [0,m]. And we are done.