0
$\begingroup$

I've been playing with the Von Mises distribution for a project I'm doing in python and I'm confused about it.

I'm drawing the pdf, which is defined by wikipedia here as $p(x|\mu, k) = \frac{\exp{(k \cos(x-\mu))}}{\tau I_0(k)}$ for the angle $x$, centre $\mu$ and dispersion or concentration $k$ or $kappa$ (and $\tau = 2\pi$).

In python, I attempted to plot the pdf for different kappas to see how the variable $k$ affects the distribution, using the trapezium rule to integrate the area under the curve. When I did this, I found some funny results, specifically: The area under the curve only sums to $1$ for the case $k = 0$, sometimes the pdf goes under the x axis, the maximum value periodically spikes to very high figures too.

I think I am missing something, probably to do with the circularity and the wrapping, but I can't figure out what a good range for kappa is and how it will affect the results I draw from the distribution. Another reason why I'm confused is that wkipedia also confirms here, in the first paragraph that the area under a pdf curve should be equal to $1$. Is it because the von mises is only defined for a segment of length $\tau$? I know I shouldn't be relying on wikipedia so much, so I'm open to redirection towards better resources though I'd rather not have to deal with anything too heavy, this isn't my area of expertise.

My code is below.

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, 2, 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):
    #plot a pdf for each kappa and sum under the curve.
    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()

Thanks for any replies shedding light on this!

Edit: It turns out I was correct about my beliefs about the area under a pdf being $1$ and the equation for the von Mises distribution. My mistake was in following some sample code in the documentation at numpy.

$\endgroup$
  • 1
    $\begingroup$ The area under any pdf is always 1 and it never goes below the $x$ axis. $\endgroup$ – Matt Samuel Apr 28 '15 at 13:34
  • $\begingroup$ $\exp(\cdot)$ is positive and so the pdf couldn't possibly dip below the $x$ axis. Check your code. $\endgroup$ – Dilip Sarwate Apr 28 '15 at 14:03
  • 1
    $\begingroup$ $\exp(\cdot)$ and $I_0(\cdot)$ are both positive everywhere, so if you're seeing a negative density value then there's definitely an error in the program. Maybe you should plot some parts of your formula separately to see if they work as desired. (In particular, is sps.jn(0,kappa) really the modified Bessel function of order $0$?) $\endgroup$ – David K Apr 28 '15 at 14:05
  • $\begingroup$ @MattSamuel that's what I thought. $\endgroup$ – Ogaday Apr 28 '15 at 14:17
  • $\begingroup$ @DavidK I didn't know that about $I_0(\dot)$. sps.jn(0,kappa) does sometimes return negative values so it's not correct (eg. sps.jn(0,4.0). Thanks. Further investigation shows that I should be using sps.iv(n, kappa). Thanks for helping me out. I was following the docs and so there's obviously a mistake there. Thanks again for your time. $\endgroup$ – Ogaday Apr 28 '15 at 14:21
2
$\begingroup$

Thanks to @Davidk for clarifying the range of the modified bessel function. The documentation I was following is actually incorrect.

The modified bessel function is actually implemented by scipy.special.iv(n, k) for order $n$ and kappa $k$, instead of the function I was using.

Replacing one for the other fixes all problems I had.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.