von mises distribution in numpy.random biased
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)
participants (2)
-
killian koepsell
-
Travis E. Oliphant