Implementing Discrete Gaussian Kernel in Python?

Question

I'm looking to implement the discrete Gaussian kernel as defined by Lindeberg in his work about scale space theory.

It is defined as T(n,t) = exp(-t)*I_n(t) where I_n is the modified Bessel function of the first kind.

I am trying to implement this in Python using Numpy and Scipy but am running into some trouble.

def discrete_gaussian_kernel(t, n):
    return math.exp(-t) * scipy.special.iv(n, t)

I try plotting with:

import math
import numpy as np
import scipy
from matplotlib import pyplot as plt

def kernel(t, n):
    return math.exp(-t) * scipy.special.iv(n, t)
ns = np.linspace(-5, 5, 1000)
y0 = discrete_gaussian_kernel(0.5, ns)
y1 = discrete_gaussian_kernel(1, ns)
y2 = discrete_gaussian_kernel(2, ns)
y3 = discrete_gaussian_kernel(4, ns)
plt.plot(ns, y0, ns, y1, ns, y2, ns, y3)
plt.xlim([-4, 4])
plt.ylim([0, 0.7])

The output looks like:

From the Wikipedia article, it should look like:

I assume I'm making some really trivial mistake. :/ Any thoughts? Thanks!

EDIT: What I wrote is equivalent to scipy.special.ive(n, t). I'm pretty sure it's supposed to be the modified Bessel function of the first kind, not the second, but can someone confirm?

Answer 1

If you want to get the wikipedia plot, replace

ns = np.linspace(-5, 5, 1000)

with

ns = np.arange(-5, 5+1)

Namely, the formula you use only makes sense at the integer points. The Bessel function as a function of negative order is an oscillating function, http://dlmf.nist.gov/10.27#E2 so the plot looks fine to me.