Zipf's law and sequence A061017
Yaroslav Bulatov
http://www.research.att.com/~njas/sequences/A061017 well?
This is an interesting sequence because it's reciprocal describes
distribution of random bigrams when individual words follow Zipf's law
Rainer Rosenthal
There is a comment saying: "a(n) is also an inverse of A006218."
If Benoit Cloitre is right, then
A006218(n)=2*sum(i=1, floor(sqrt(n)), floor(n/i))-floor(sqrt(n))^2
HTH
Rainer Rosenthal
r.ros...@web.de
David W. Cantrell
> Yaroslav Bulatov schrieb:
> > Can anyone see a closed form expression that approximates
> > http://www.research.att.com/~njas/sequences/A061017 well?
> >
> > This is an interesting sequence because it's reciprocal describes
> > distribution of random bigrams when individual words follow Zipf's law
> >
>
> There is a comment saying: "a(n) is also an inverse of A006218."
Yes, I see the comment, but I don't understand it. Can anyone explain to me
the sense in which A061017 is "an inverse of A006218" ?
David
David W. Cantrell
OK, I now understand. That allows me to give a crude approximation which
might satisfy Yaroslav:
For A061017, a(n) is approximately
n / W(e^(2g - 1) n)
where W denotes the principal branch of the Lambert W function (see
<http://mathworld.wolfram.com/LambertW-Function.html>, for example)
and g denotes the Euler-Mascheroni constant.
I do not know how good the approximation is for large n. But for the 77
entries given in A061017, we have -0.70 < error < 0.91 .
David W. Cantrell
Rainer Rosenthal
>>> There is a comment saying: "a(n) is also an inverse of A006218."
>> Yes, I see the comment, but I don't understand it. Can anyone explain to
>> me the sense in which A061017 is "an inverse of A006218" ?
>
> OK, I now understand. That allows me to give a crude approximation which
> might satisfy Yaroslav:
>
> For A061017, a(n) is approximately
>
> n / W(e^(2g - 1) n)
>
> where W denotes the principal branch of the Lambert W function (see
> <http://mathworld.wolfram.com/LambertW-Function.html>, for example)
> and g denotes the Euler-Mascheroni constant.
>
Bravo!
> I do not know how good the approximation is for large n. But for the 77
> entries given in A061017, we have -0.70 < error < 0.91 .
So Neill might find another comment after returning from his
holiday?
Cheers,
Rainer Rosenthal
r.ros...@web.de
David W. Cantrell
> David W. Cantrell schrieb:
>
> >>> There is a comment saying: "a(n) is also an inverse of A006218."
> >> Yes, I see the comment, but I don't understand it. Can anyone explain
> >> to me the sense in which A061017 is "an inverse of A006218" ?
> >
> > OK, I now understand. That allows me to give a crude approximation
> > which might satisfy Yaroslav:
> >
> > For A061017, a(n) is approximately
> >
> > n / W(e^(2g - 1) n)
> >
> > where W denotes the principal branch of the Lambert W function (see
> > <http://mathworld.wolfram.com/LambertW-Function.html>, for example)
> > and g denotes the Euler-Mascheroni constant.
>
> Bravo!
Thanks, Rainer. But my approximation does not deserve a bravo. Just look at
its error bounds.
> > I do not know how good the approximation is for large n. But for the 77
> > entries given in A061017, we have -0.70 < error < 0.91 .
>
> So Neill might find another comment after returning from his holiday?
Yes. However, perhaps someone can find a better approximation...
Best regards,
David
David W. Cantrell
> David W. Cantrell <DWCan...@sigmaxi.net> wrote:
> > Rainer Rosenthal <r.ros...@web.de> wrote:
> > > Yaroslav Bulatov schrieb:
> > > > Can anyone see a closed form expression that approximates
> > > > http://www.research.att.com/~njas/sequences/A061017 well?
> > > >
> > > > This is an interesting sequence because it's reciprocal describes
> > > > distribution of random bigrams when individual words follow
> > > > Zipf's law
> > >
> > > There is a comment saying: "a(n) is also an inverse of A006218."
> >
> > Yes, I see the comment, but I don't understand it. Can anyone explain
> > to me the sense in which A061017 is "an inverse of A006218" ?
>
> OK, I now understand.
In a private email, Rainer suggested that I explain the sense in which
A061017 is an inverse of A006218. Here's the explanation:
A006218 <http://www.research.att.com/~njas/sequences/A006218>,
after omitting its 0th element, is 1, 3, 5, 8, 10, 14, ...
Its inverse function, shown in terms of ordered pairs, is
(1,1), (3,2), (5,3), (8,4), (10,5), (14,6), ...
Let's compare. On the first line below, I show the inverse of A006218,
now as a "sequence with some missing elements"; on the second line, I
show A061017:
1, _, 2, _, 3, _, _, 4, _, 5, _, _, _, 6, _, ...
1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, ...
I trust that the comparison is obvious.
Thus, the inverse of A006218 is a particular _sub_sequence of A061017.
In any event, I don't like saying that A061017 is an inverse of A006218; I
found that confusing. Rather, I think it would be better if the OEIS entry
said something like "A subsequence of A061017 is the inverse of A006218."
And it would be better still if a description of that particular
subsequence were given.
Best regards,
David
Rainer Rosenthal
> Let's compare. On the first line below, I show the inverse of A006218,
> now as a "sequence with some missing elements"; on the second line, I
> show A061017:
>
> 1, _, 2, _, 3, _, _, 4, _, 5, _, _, _, 6, _, ...
> 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, ...
>
> I trust that the comparison is obvious.
> Thus, the inverse of A006218 is a particular _sub_sequence of A061017.
Wonderful, thank you! It seems as if
A006218(n) = A061017(k)
implies (1)
A006218(n+1) = A06107(k+1)
which would be nice. But (1) could be a case of the law of small numbers.
Maybe there is a solution manual for it :-/
Cheers,
Rainer Rosenthal
r.ros...@web.de
Rainer Rosenthal
>
> A006218(n) = A061017(k)
> implies (1)
> A006218(n+1) = A06107(k+1)
>
Ooops, now I am officially embarrassed. I will have to think it over,
what I was trying to say :-(
Cheers,
RR
Rainer Rosenthal
> Let's compare. On the first line below, I show the inverse of A006218,
> now as a "sequence with some missing elements"; on the second line, I
> show A061017:
>
> 1, _, 2, _, 3, _, _, 4, _, 5, _, _, _, 6, _, ...
> 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, ...
>
> I trust that the comparison is obvious.
> Thus, the inverse of A006218 is a particular _sub_sequence of A061017.
These are the definitions from the OEIS:
n: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
------------------------------------------------------------------------
f = A006218: 1 3 5 8 10 14 16 20 23 27 29 35 37 41 45 50 ...
g = A061017: 1 2 2 3 3 4 4 4 5 5 6 6 6 6 7 7 ...
f(n) is the largest k such that g(k) = n. I.e. f(n) = max g^(-1)[{n}].
Therefore g(f(n)) = n.
And from the "largest" property and the slow rise of g we get:
g(f(n)+1) = n+1. If we define f'(n) = f(n-1)+1 then g(f'(n)) = n also.
So A061017 inverts A006218 but not vice versa.
Cheers,
Rainer Rosenthal
r.ros...@web.de
David W. Cantrell
> David W. Cantrell wrote:
>
> > Let's compare. On the first line below, I show the inverse of A006218,
> > now as a "sequence with some missing elements"; on the second line, I
> > show A061017:
> >
> > 1, _, 2, _, 3, _, _, 4, _, 5, _, _, _, 6, _, ...
> > 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, ...
> >
> > I trust that the comparison is obvious.
> > Thus, the inverse of A006218 is a particular _sub_sequence of A061017.
>
> These are the definitions from the OEIS:
>
> n: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
> ------------------------------------------------------------------------
> f = A006218: 1 3 5 8 10 14 16 20 23 27 29 35 37 41 45 50 ...
> g = A061017: 1 2 2 3 3 4 4 4 5 5 6 6 6 6 7 7 ...
Note that I added another row above.
> f(n) is the largest k such that g(k) = n.
Perfectly clearly said!
> I.e. f(n) = max g^(-1)[{n}]. Therefore g(f(n)) = n.
Recalling that Bottomley made the comment that A061017(n) "is also an
inverse of A006218." and noting that he is also the author of A054519
<http://www.research.att.com/~njas/sequences/A054519>, I'm slightly
surprised that he didn't make a similar comment about it. That is, he might
have said A061017(n) is also an inverse of A054519. But the clearest thing
to have said would have been, following what you said earlier,
A054519(n) is the smallest k such that A061017(k) = n.
Cheers,
David
Rainer Rosenthal
>> These are the definitions from the OEIS:
>>
>> n: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...
>> ------------------------------------------------------------------------
>> f = A006218: 1 3 5 8 10 14 16 20 23 27 29 35 37 41 45 50 ...
>> g = A061017: 1 2 2 3 3 4 4 4 5 5 6 6 6 6 7 7 ...
>| h = A054519: 1 2 4 6 9 11 15 17 21 24 28 30 36 38 42 46 ...
>
> Note that I added another row (h=A054519) above.
>
>> f(n) is the largest k such that g(k) = n.
>> g(f(n)) = n.
>
> Recalling that Bottomley made the comment that A061017(n) "is also an
> inverse of A006218." and noting that he is also the author of A054519
> ... he might have said A061017(n) is also an inverse of A054519.
>
> A054519(n) is the smallest k such that A061017(k) = n.
Nice detective work.
>> And from the "largest" property and the slow rise of g we get:
>> g(f(n)+1) = n+1. If we define f'(n) = f(n-1)+1 then g(f'(n)) = n also.
This was just a funny remark of me. If I hadn't been lazy I would have
looked up the OEIS for this function f'. And then I would have been
very happy about f'(n) = f(n-1)+1 or, to put it more bombastically,
-1
f' = A o f o A where A(x) = x+1. Anyhow this is a good observation:
A054519(n) = A006218(n-1) + 1
and
A006218(n) = A054519(n+1) - 1
Back to the question of the OP: it seems as if this other inverse function
f' = h = A054519 could be helpful in computing A061017. I am not at all
good in Lambert W and such things but I feel that we have "the other end"
now.
Best regards,
Rainer Rosenthal
r.ros...@web.de
Yaroslav Bulatov
errors of this approximation have a peculiar "ringed" appearance
http://yaroslavvb.com/temp/lambert-errors.png
ll = Flatten[Table[Table[n, {Length[Divisors[n]]}], {n, 1, 10000}]];
cantrell[n_] := n/ProductLog[Exp[(2 EulerGamma - 1)] n]
errors = cantrell[#] - ll[[#]] & /@ Range[1, 10000] // N;
ListPlot[errors]
Rainer Rosenthal
This looks like a success for your question. Congratulation!
Do you see a refinement in using the other 'inverse function'?
As I said already, there might be some potential in using
the function h(n), given by David Cantrell:
h = A054519: 1 2 4 6 9 11 15 17 21 24 28 30 36 38 42 46 ...
This function is very much like
f = A006218: 1 3 5 8 10 14 16 20 23 27 29 35 37 41 45 50 ...
because h(n) = f(n+1)-1. The inversion of f led David to the
formula that helped you. It should be helpful to invert h
as well and use it, too. As I said, I am not good in these
manipulations but I /feel/ that the inversion of h should be
a straightforword matter for someone who managed to invert f.
David W. Cantrell
> Thanks, that seems like a good approximation. Curiously, plot of the
> errors of this approximation have a peculiar "ringed" appearance
>
> http://yaroslavvb.com/temp/lambert-errors.png
Yes, the appearance is peculiar, but I'm not surprised.
> ll = Flatten[Table[Table[n, {Length[Divisors[n]]}], {n, 1, 10000}]];
> cantrell[n_] := n/ProductLog[Exp[(2 EulerGamma - 1)] n]
> errors = cantrell[#] - ll[[#]] & /@ Range[1, 10000] // N;
> ListPlot[errors]
Earlier in this thread, I had said that I didn't know how my approximation
behaves for large n. But now, a final comment:
Although error surely increases without bound, relative error surely goes
to 0. In fact, I now conjecture, crudely, that
| relative error | < 1/(n^(1/3) log(n+1))
for all n.
David W. Cantrell