MaxExtraPrecision = 20000000;
LengthOfCF = 10000;
NumberOfCF = 1000;
MaxDigitCount = 101;
\(\(minpoly = \  (*actually\ minimal\ polynomial\ or\ irreducible\ polynomial*) 
\(CFTable = 
          Table[ContinuedFraction[Root[minpoly, 1], LengthOfCF], {i, 
              NumberOfCF}]\);\)
 (*The\ above\ line\ simply\ generates\ the\ desired\ continued\ fraction\ using\ Mathematica' 
        s\ algorithm\ for\ generating\ an\ arbitrary\ finite\ number\ of\ digits\ of\ a\ continued\ fraction\ representation\ of\ a\ number\ and\ stores\ the\ result\ in\ a\ table*) \)
MatrixForm[CFTable];
a = 
        Table[Count[CFTable[[m]], n], {m, NumberOfCF}, {n, 
            MaxDigitCount}];
  
MatrixForm[a]

<< Statistics`NormalDistribution`
<< Graphics`MultipleListPlot`
Clear[i, j, m, n];
DigitMatrix = 
        Transpose[a];

ScaledDigitFrequency = 
        Table[N[Mean[DigitMatrix[[j]]]/LengthOfCF, 5], {j, 
            MaxDigitCount}];
\(\(Kuzmin = 
        Table[N[Log[1 + 1/\((j*\((j + 2)\))\)]/Log[2], LengthOfCF], {j, 
            MaxDigitCount}];\)
\)
PlotTest1 = 
      MultipleListPlot[Kuzmin, ScaledDigitFrequency, 
        PlotJoined -> {True, False}, 
        SymbolStyle -> {Hue[1], Hue[ .6]}, 
        PlotLegend -> {"Theoretical Frequency", "Tabulated \
Frequency"}, LegendPosition -> {0, 0}, 
        Frame -> 
          True]
  
PlotTest2 = 
      MultipleListPlot[
        Table[N[Variance[
              DigitMatrix[[j]]/
                  LengthOfCF - (Log[1 + 1/(j*(j + 2))])/
                  Log[2]]], {j, MaxDigitCount}], 
        PlotLegend -> {"Variance from Kuz'min's theorem"}, 
        LegendPosition -> {0, 0}, LegendSize -> {1.2,  .5}, 
        Frame -> 
          True]