# [Numpy-discussion] numpy.random.logseries - incorrect convergence for k=1, k=2

joep josef.pktd@gmail....
Fri Oct 3 00:17:28 CDT 2008

Filed as http://scipy.org/scipy/numpy/ticket/923

and I think i finally tracked down the source of the incorrect random
numbers, a reversed inequality in
http://scipy.org/scipy/numpy/browser/trunk/numpy/random/mtrand/distributions.c
line 871, see my last comment to the trac ticket.

Josef

On Sep 27, 2:12 pm, joep <josef.p...@gmail.com> wrote:
> random numbers generated by numpy.random.logseries do not converge to
> theoretical distribution:
>
> for probability paramater pr = 0.8, the random number generator
> converges to a
> frequency for k=1 at 39.8 %, while the theoretical probability mass is
> 49.71
> k=2 is oversampled, other k's look ok
>
> check frequency of k=1 and k=2 at N =  1000000
> 0.398406 0.296465
> pmf at k = 1 and k=2 with formula
> [ 0.4971  0.1988]
>
> for probability paramater pr = 0.3, the results are not as bad, but
> still off:
> frequency for k=1 at 82.6 %, while the theoretical probability mass is
> 84.11
>
> check frequency of k=1 and k=2 at N =  1000000
> 0.826006 0.141244
> pmf at k = 1 and k=2 with formula
> [ 0.8411  0.1262]
>
> below is a quick script for checking this
>
> Josef
>
> {{{
> import numpy as np
> from scipy import stats
>
> pr = 0.8
> np.set_printoptions(precision=2, suppress=True)
>
> # calculation for N=1million takes some time
> for N in [1000, 10000, 10000, 1000000]:
>     rvsn=np.random.logseries(pr,size=N)
>     fr=stats.itemfreq(rvsn)
>     pmfs=stats.logser.pmf(fr[:,0],pr)*100
>     print 'log series sample frequency and pmf (in %) with N = ', N
>     print np.column_stack((fr[:,0],fr[:,1]*100.0/N,pmfs))
>
> np.set_printoptions(precision=4, suppress=True)
>
> print 'check frequency of k=1 and k=2 at N = ', N
> print np.sum(rvsn==1)/float(N),
> print np.sum(rvsn==2)/float(N)
>
> k = np.array([1,2])
> print 'pmf at k = 1 and k=2 with formula'
> print -pr**k * 1.0 / k / np.log(1-pr)}}}
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discuss...@scipy.orghttp://projects.scipy.org/mailman/listinfo/numpy-discussion