A241773 A sequence constructed so that the probability of occurrence of integer i > 0 is given by p(i) = log_2[(i+1)^2/(i^2+2*i)], which is the Gauss-Kuzmin distribution.

1, 2, 3, 1, 4, 1, 5, 2, 1, 6, 1, 7, 1, 2, 3, 1, 8, 1, 9, 2, 1, 4, 1, 10, 1, 3, 2, 1, 11, 1, 2, 5, 1, 12, 1, 3, 1, 2, 4, 1, 13, 1, 2, 6, 1, 14, 1, 3, 1, 2, 15, 1, 16, 1, 2, 4, 1, 3, 1, 5, 7, 1, 2, 1, 17, 1, 2, 3, 1, 18, 1, 8, 2, 1, 4, 1, 6, 1, 2, 3, 1, 5, 1, 19

Offset

1,2

Comments

Let the sequence be A = {a(i)}, i = 1, 2, 3,... and define p(i) =

log_2[(i+1)^2/(i^2+2*i)]. Additionally, define u(j, k) = k*p(j) - N(j, k), where N(j, k) is the number of occurrences of j in {a(i)}, i = 1,..., k-1. Refer to the first argument of u as the "index" of u. Then A is defined by a(1) = 1 and, for i > 1, a(i) = m, where m is the index of the maximal element of the set {u(j, i)}, j = 1, 2, 3,... That there is a single maximal element for all i is guaranteed by the fact that p(i) - p(j) is irrational for all i not equal to j.

Interpreting sequence A as the partial coefficients of the continued fraction expansion of a real number C, say, then C = 1.44224780173510148... which is, by construction, normal (in the continued fraction sense).

The geometric mean of the sequence equals Khintchine's constant K=2.685452001 = A002210 since the frequency of the integers agrees with the Gauss-Kuzmin distribution. - Jwalin Bhatt, Feb 11 2025

Links

Jonathan Deane, Table of n, a(n) for n = 1..10000

Jonathan Deane, An integer sequence whose members obey a given p.d.f.

Maple

pdf := i -> -log[2](1 - 1/(i+1)^2);
gen_seq := proc(n)
local i, j, N, A, u, mm, ndig;
ndig := 40; N := 'N';
for i from 1 to n do    N[i] := 0;      end do;
A := 'A'; A[1] := 1; N[1] := 1;
for i from 2 to n do
        u := 'u';
        for j from 1 to n do
                u[j] := i*pdf(j) - N[j];
        end do;
        mm := max_maxind(evalf(convert(u, list), ndig));
        if mm[3] then
                A[i] := mm[1];
                N[mm[1]] := N[mm[1]] + 1;
        else
                return();
        end if;
end do;
return(convert(A, list));
end:
max_maxind := proc(inl)
local uniq, mxind, mx, i;
uniq := `true`;
if nops(inl) = 1 then return([1, inl[1], uniq]); end if;
mxind := 1; mx := inl[1];
for i from 2 to nops(inl) do
        if inl[i] > mx then
                mxind := i;
                mx := inl[i];
                uniq := `true`;
        elif inl[i] = mx then
                uniq := `false`;
        end if;
end do;
return([mxind, mx, uniq]);
end:
gen_seq(100);

Prog

(Python)
from fractions import Fraction
def prob_count_diff(j, k, count):
    return  - ((1 - Fraction(1, (j+1)*(j+1)))**k) * (2**count)
def sample_pdf(n):
    coeffs, unreached_val, counts = [], 1, {}
    for k in range(1, n+1):
        prob_count_diffs = [prob_count_diff(i, k, counts.get(i, 0)) for i in range(1, unreached_val+1)]
        most_probable = prob_count_diffs.index(max(prob_count_diffs)) + 1
        unreached_val += most_probable == unreached_val
        coeffs.append(most_probable)
        counts[most_probable] = counts.get(most_probable, 0) + 1
    return coeffs
A241773 = sample_pdf(120)  # Jwalin Bhatt, Feb 09 2025

Crossrefs

Cf. A002210, A084580, A089618.

Sequence in context: A286477 A277230 A218534 * A348388 A205790 A279820

Adjacent sequences: A241770 A241771 A241772 * A241774 A241775 A241776

Keyword

cofr,easy,nonn

Author

Jonathan Deane, Apr 28 2014

Status

approved