[Numpy-discussion] von mises distribution in numpy.random biased

killian koepsell koepsell@gmail....
Sat Sep 8 17:39:23 CDT 2007


hi,

the von mises distribution in numpy.random seems to be biased towards
a higher concentration (kappa). given a concentration of 2, it
produces data that has a concentration of 2.36. i compared the
distribution to the one produced by the CircStats[1] package of R[2]
using RPy [3] and created a figure here:

  http://redwood.berkeley.edu/kilian/vonmises.png

the script i used is attached to this email. i don't know what
algorithm NumPy uses, so i can't tell if it is a real bug or some sort
of rounding error. the CircStats package uses the algorithm by Best
and Fisher [4].

kilian

[1] http://cran.r-project.org/src/contrib/Descriptions/CircStats.html
[2] http://www.r-project.org/
[3] http://rpy.sf.net/
[4] Best, D. and Fisher, N. (1979). Efficient simulation of the von
Mises distribution. Applied Statistics, 24, 152-157.

----------------------
vonmises.py
----------------------
from numpy import array,pi,exp,mod,random
from pylab import hist,linspace,figure,clf,plot,subplot,xlim,title,legend
from rpy import r as R
R.library("CircStats")

# compute von Mises distribuion fo two different values of kappa
size = 1e5
kappa1 = 1
x1 = random.vonmises(0,kappa1,size=size)
y1 = mod(R.rvm(size,0,kappa1),2*pi)
kappa2 = 2
x2 = random.vonmises(0,kappa2,size=size)
y2 = mod(R.rvm(size,0,kappa2),2*pi)

def phasedist(x,kappa):
    ph = linspace(-pi,pi,50)
    vM = R.dvm(ph,0,kappa)
    kest = R.A1inv(abs(exp(1j*array(x)).mean()))
    vMest = R.dvm(ph,0,kest)
    plot(ph,vM,'r',linewidth=2)
    plot(ph,vMest,'k:',linewidth=2)
    legend(['$\\kappa=%d$'%kappa,'$\\hat{\\kappa}=%1.2f$'%kest])
    hist(x,ph,normed=1)
    xlim((-pi,pi))

# plot figure
fig=figure(1)
clf()
ax1=subplot(2,2,1)
title('NumPy.random')
phasedist(x1,kappa1)
ax2=subplot(2,2,3)
phasedist(x2,kappa2)
subplot(2,2,2,sharey=ax1)
title('RPy.CircStats')
phasedist(y1,kappa1)
subplot(2,2,4,sharey=ax2)
phasedist(y2,kappa2)


More information about the Numpy-discussion mailing list