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).

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);

Crossrefs

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