Repository navigation
numpy.random.vonmises documentation contains incorrect code example. #5811
Description
The numpy.random reference for the von Mises distribution contains an incorrect code sample (pasted below).
mu, kappa = 0.0, 4.0 # mean and dispersion
s = np.random.vonmises(mu, kappa, 1000)
import matplotlib.pyplot as plt
import scipy.special as sps
count, bins, ignored = plt.hist(s, 50, normed=True)
x = np.arange(-np.pi, np.pi, 2*np.pi/50.)
y = -np.exp(kappa*np.cos(x-mu))/(2*np.pi*sps.jn(0,kappa))
plt.plot(x, y/max(y), linewidth=2, color='r')
plt.show()While this works for the example, line 5 doesn't concur with the definition of the von Mises probability density function, given above the example. According to documentation about the special functions, the modified bessel function of order real order is given by scipy.special.iv(). Using that function instead means the sign doesn't need to be changed.
I noticed the problem when I tried to emulate the code and I found that area under the pdf didn't sum to 1 when I used the trapezium rule and would sometimes even go below the x axis, when neither the exponential function or the bessel function should ever be negative. Compare:
#incorrect, according to documentation
import numpy as np
import scipy.special as sps
import matplotlib.pyplot as plt
tau = np.pi*2
fig, ax = plt.subplots(1,1, figsize = (10,8))
n = 10
xs = np.linspace(0, 4, n) #create a range of kappas
mu = 0
m = 50
width = tau/float(m)
sums = np.empty(n) #for storing the integrals under the curves for different kappas
for i in range(n):
kappa = xs[i]
x = np.linspace(-np.pi, np.pi, m)
y = -np.exp(kappa*np.cos(x-mu))/(tau*sps.jn(0,kappa))
ax.plot(x, y, label="kappa = {}".format(kappa))
sums[i] = sum([(y[j]+y[j+1])/2*width for j in range(m-1)])
ax.set_axes("tight")
ax.legend()
print "the area under the pdf for different kappas: "
for i in range(n):
print "kappa =", xs[i], ", integral =", sums[i]
fig.show()#correct, using the correct bessel function
import numpy as np
import scipy.special as sps
import matplotlib.pyplot as plt
tau = np.pi*2
fig, ax = plt.subplots(1,1, figsize = (10,8))
n = 10
xs = np.linspace(0, 4, n) #create a range of kappas
mu = 0
m = 50
width = tau/float(m)
sums = np.empty(n) #for storing the integrals under the curves for different kappas
for i in range(n):
kappa = xs[i]
x = np.linspace(-np.pi, np.pi, m)
y = np.exp(kappa*np.cos(x-mu))/(tau*sps.iv(0,kappa))
ax.plot(x, y, label="kappa = {}".format(kappa))
sums[i] = sum([(y[j]+y[j+1])/2*width for j in range(m-1)])
ax.set_axes("tight")
ax.legend()
print "the area under the pdf for different kappas: "
for i in range(n):
print "kappa =", xs[i], ", integral =", sums[i]
fig.show()Excuse my use of tau.
Using: Python 2.7.9 :: Anaconda 2.2.0 (64-bit), IPython 3.1.0, numpy 1.9.2
Activity
argriffing commentedon Apr 28, 2015
For what it's worth,
scipy.stats.vonmisesusesi0for the pdf.https://github.com/scipy/scipy/blob/maintenance/0.15.x/scipy/stats/_continuous_distns.py#L3956
The framework for these 'circular' or 'wrapped' distributions could be improved in scipy; I think they are implemented as hacks on top of continuous location-scale distributions on the real line.
Ogaday commentedon Apr 28, 2015
Oh, interesting, thanks for that. There do seem to be an overload options between numpy and scipy and I often don't know which to use.
I was fairly certain that the sampling method itself was fine, it's the just example code. Also, the links at the bottom of that page are dead (the source, png, pdf links).
argriffing commentedon Apr 28, 2015
@Ogaday Just so I understand, changing
jn(0,kappa)to-i0(kappa)in the example code and updating the plot and fixing the dead links would be enough to close this issue?Ogaday commentedon Apr 29, 2015
Yes, and the - sign there cancels with the one at the start of the term, so replace
jn(0,kappa)withi0(kappa)and delete the-sign aftery =.Technically, in terms of best practice, I think it should also read:
Which does the same thing but with
linspaceinstead ofarange, as suggested in the docstring ofarange, and it includes the end of the interval, which should be defined under the definition of the function. That is a little pedantic though, but I think it's better practice.argriffing commentedon Apr 29, 2015
@Ogaday could you check that the new docstring in #5814 works for you locally? I think the plot image in the docs would be automatically generated from this docstring code, and I've opened a separate issue for the problem with the dead links.
Ogaday commentedon Apr 30, 2015
Works great for me locally, thanks. I was also wondering about the normalising thing as well, looks much better now.