Why is the area under the pdf for the Von Mises distribution not one?

Question

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.

Answer 1

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.