Gauss-Kuzmin Distribution

Author

Eric W. Weisstein
September 6, 2008

This notebook downloaded from http://mathworld.wolfram.com/notebooks/Statistics/Gauss-KuzminDistribution.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/Gauss-KuzminDistribution.html.

©2008 Wolfram Research, Inc. except for portions noted otherwise

Code

GaussKuzminList[x_,n_,d_:1000]:=Table[Count[Rest[ContinuedFraction[x,d+1]],k],{k,n}]/d

Calculation of x(n)

x[-1,r_]:=r
a[n_,r_]:=a[n,r]=Floor[1/x[n-1,r]]
x[n_,r_]:=x[n,r]=1/x[n-1,r]-a[n,r]

...this is very inefficient...

Xi[r_]:=Xi[r]=ContinuedFraction[r,5000]
Xi[r_,n_]:=FromContinuedFraction[Drop[Xi[r],n-1]]//N

ContinuedFraction[Pi,10]

Table[a[n,Pi],{n,5}]

1/(a[1,Pi]+x[1,Pi])

1/(a[1,Pi]+1/(a[2,Pi]+x[2,Pi]))

Table[x[n,Pi],{n,5}]

FractionalPart[1/Table[x[n,Pi],{n,0,10}]]//N

FractionalPart[Table[x[n,Pi],{n,10}]]//N

FractionalPart/@(Xi[Pi,#]&/@Range[10])

Probability Function

D[Log[2,1+x],x]

Integrate[(1/((1+x) Log[2])),{x,0,1}]

Plot[(1/((1+x) Log[2])),{x,0,1},PlotStyle->Red,AxesLabel->TraditionalForm/@{x,lg x},PlotRange->{0,1.5}]

pidata=FractionalPart/@(Xi[Pi,#]&/@Range[2000]);//Timing

pidata=FractionalPart/@(x[#,N[Pi,10000]]&/@Range[5000]);//Timing

eulerdata=FractionalPart/@(x[#,N[EulerGamma,10000]]&/@Range[5000]);//Timing

catalandata=FractionalPart/@(x[#,N[Catalan,10000]]&/@Range[5000]);//Timing

lndata=FractionalPart/@(x[#,N[Log[2],10000]]&/@Range[5000]);//Timing

Bins[data_,n_:10]:=Module[{l={First[#],Length[#]}&/@Split[Sort[Floor[n data]]]},
Last/@Sort[Join[l,{#,0}&/@Complement[Range[n-1],First/@l]]]
]

labels=#[[1]]<>"–"<>#[[2]]&/@Partition[ToString[NumberForm[#,{2,1},NumberPadding->{"","0"}]]&/@Range[0.,1,.1],2,1]

<<BarCharts`

With[{colors={Red,Yellow,Green,Blue}},BarChart[Sequence@@(Bins[#,10]&/@{pidata,eulerdata,catalandata,lndata}),BarSpacing->0,BarLabels->labels,BarStyle->colors,
AxesLabel->TraditionalForm/@{"",StyleForm[#[[1]],FontColor->#[[2]]]&/@Transpose[{{Pi,EulerGamma,Catalan,Log[2]},{Red,Yellow,Green,Blue}}]},ImageSize->500]
]

Distribution Function

Plot[Log[2,1+x],{x,0,1},PlotStyle->Red,AxesLabel->TraditionalForm/@{x,lg x}]

Log[2,1+x]/.x->1

Series[Log[2,1+x],{x,0,2}]

Normalization

-Sum[Log[2,1-(1/(((k+1))^2))],{k,∞}]

C.F. Terms

<<Graphics`Graphics`

data=FromContinuedFraction[Prepend[#,0]]&/@Table[Random[Integer,{1,1000}],{1000},{10}];//N

{First[#],Length[#]}&/@Split[Sort[Floor[100data]]]

Histogram[Length/@Split[Sort[Floor[10data]]],FrequencyData->True]

Distribution of (Subscript[a, k])

Bailey

Sum[-Log[2,1-1/(n+1)^2],{n,∞}]

FullSimplify[Subtract@@{-Log[2,1-1/(n+1)^2],Log[2,1+(1/(n(n+2)))]},n>0]

Havil

Sum[Log[2,1+(1/(k(k+2)))],{k,∞}]//FullSimplify

N[%]

NSum[Log[2,1+(1/(k(k+2)))],{k,∞}]//Chop

Plot

<<Utilities`Plot`

c=FromDigits[{Flatten[IntegerDigits/@Prime[1000]],0}];

Show[GraphicsArray[{
IntegerListPlot[d=Table[-Log[2,1-1/(n+1)^2],{n,20}],
PlotRange->All,
AxesLabel->TraditionalForm/@{n,-HoldForm[Log[2,1-InlineFraction[1,(n+1)^2]]]},
DisplayFunction->Identity],
Graphics[ListPlot[GaussKuzminList[#[[1]],20,500],
PlotStyle->{PointSize[.03],#[[2]]},
DisplayFunction->Identity][[1]]&/@Transpose[{{Pi,Sin[1],EulerGamma,c},{Red,Green,Blue,Yellow}}],
Axes->True,
Epilog->{IntegerListPlot[d,PlotStyle->Black,DisplayFunction->Identity][[1]]},
AxesLabel->TraditionalForm/@{n,HoldForm[P[(Subscript[a, k])=n]]}]
}
]
]

Calculation of (Subscript[F, n])(x)

Off[Divide::"infy",∞::"indet"]

With[{its=100,n=50},Table[Count[Table[x[n,i+10^(-40)Random[Real,{0,1},100]],{its}],_?(#≤i&)]/its,{i,0,1,1/10}]
]//N//Timing

With[{its=200,n=50},Table[Count[Table[x[n,i+10^(-40)Random[Real,{0,1},100]],{its}],_?(#≤i&)]/its,{i,0,1,1/10}]
]//N//Timing

Table[Log[2,1+i],{i,0,1,1/10}]//N