[SciPy-User] How do I use vonmises.fit()?

Tue Apr 13 14:52:35 EDT 2010

Setting scipy.stats.distributions.vonmises.a and
scipy.stats.distributions.vonmises.b seems to work well; thanks for your
help!
For the record, here's what I ended up doing.
I had to manually shift the data so that the mean was 0, perform the
fitting, and then shift everything back.
(Actually, I originally expected vonmises.fit() to do something like that in
the first place.)

>>> import matplotlib.pyplot as mp
>>> import scipy.stats
>>> scipy.stats.distributions.vonmises.a = -numpy.pi
>>> scipy.stats.distributions.vonmises.b = numpy.pi
>>>
>>> def fit_vonmises(spike_phases):
>>>
>>>     # the values in spike_phases originally lie between 0 and 2pi.
>>>     # we want the values in shifted_spike_phases to lie between -pi and
pi, with the mean at 0.
>>>     mean_spike_phase = scipy.stats.circmean(spike_phases)
>>>     shifted_spike_phases = numpy.angle(numpy.exp(1j*(spike_phases -
mean_spike_phase)))
>>>     # that's just the "circular difference" between spike_phases and
mean_spike_phase.
>>>
>>>     # plot spike_phases and shifted_spike_phases.
>>>     mp.figure()
>>>     mp.subplot(211)
>>>     mp.hist(spike_phases, numpy.linspace(0, 2*numpy.pi, 37),
facecolor='b')
>>>     mp.axvline(mean_spike_phase, color='r')
>>>     mp.axis(xmin=-numpy.pi, xmax=2*numpy.pi)
>>>     mp.subplot(212)
>>>     mp.hist(shifted_spike_phases, numpy.linspace(-numpy.pi, numpy.pi,
37), facecolor='g')
>>>     mp.axvline(0, color='r')
>>>     mp.axis(xmin=-numpy.pi, xmax=2*numpy.pi)
>>>
>>>     # fit a von mises distribution to spike_phases_shifted.
>>>     pars = scipy.stats.distributions.vonmises.fit(shifted_spike_phases)
>>>     xs = numpy.linspace(-numpy.pi, numpy.pi, 361)
>>>     ys = scipy.stats.distributions.vonmises.pdf(xs, pars[0],
loc=pars[1], scale=pars[2])
>>>
>>>     # plot the fitted distribution on top of shifted_spike_phases.
>>>     mp.figure()
>>>     mp.hist(shifted_spike_phases, numpy.linspace(-numpy.pi, numpy.pi,
37), facecolor='g')
>>>     mp.axis(xmin=-numpy.pi, xmax=numpy.pi)
>>>     mp.twinx()
>>>     mp.axis(xmin=-numpy.pi, xmax=numpy.pi)
>>>     mp.plot(xs, ys, 'r') # but now, the indices of ys match up with xs,
which goes from -pi to pi.
>>>
>>>     # since the original data goes from 0 to 2pi, we have to shift
everything back by mean_spike_phase.
>>>     shifted_xs = numpy.mod(xs + mean_spike_phase, 2*numpy.pi) # use
mod() to shift around a circle.
>>>     shifted_indices = numpy.argsort(shifted_xs)
>>>     shifted_ys = ys[shifted_indices]
>>>
>>>     # now, plot the fitted distribution on top of the original
spike_phases.
>>>     mp.figure()
>>>     mp.hist(spike_phases, numpy.linspace(0, 2*numpy.pi, 37),
facecolor='b')
>>>     mp.axis(xmin=0, xmax=2*numpy.pi)
>>>     mp.twinx()
>>>     mp.plot(shifted_xs[shifted_indices], shifted_ys, 'r')
>>>     mp.axis(xmin=0, xmax=2*numpy.pi)
>>>
>>>     return (pars[0], mean_spike_phase, pars[2])
>>>
>>> spike_phases = numpy.mod(scipy.stats.distributions.vonmises.rvs(1.23,
loc=5.67, scale=1, size=1000), 2*numpy.pi)
>>> fit_vonmises(spike_phases)

It still seems a little bit hacky, so let me know if you can think of a
better way. =)

--David Ho

On Tue, Mar 30, 2010 at 1:51 PM, <josef.pktd at gmail.com> wrote:
>
> On Tue, Mar 30, 2010 at 4:34 PM,  <josef.pktd at gmail.com> wrote:
> > On Tue, Mar 30, 2010 at 3:10 PM, David Ho <itsdho at ucla.edu> wrote:
> >> Hi all!
> >>
> >> I want to fit some data with a Von Mises distribution, and get the mean
and
> >> "kappa" parameters for that distribution.
> >> I'm a little confused about scipy.stats.distributions.vonmises.fit().
> >>
> >> Just as a test, I tried this:
> >>
> >>>>> vm_rvs = scipy.stats.vonmises.rvs(1.2, 2.3, size=10000) # the first
> >>>>> argument appears to be "kappa", and the second argument is the mean.
> >>>>> scipy.stats.distributions.vonmises.fit(vm_rvs)
> >> array([  1.17643696e-01,   3.38956854e-03,   1.27331662e-27])
> >>
> >> I got an array of 3 values, none of which seem to be equal to the mean
or
> >> kappa of the distribution.
> >> I looked around in some of the docstrings, but I couldn't find any
clear
> >> documentation of what these values are supposed to correspond to.
> >>
> >> I would've expected vonmises.fit() to return something like:
> >> array([  1.2,   2.3])
> >>
> >> What am I doing wrong?
> >> Thanks for your help,
> >
> > When I did I check of the fit method for all distributions, then
> > vonmises most of the time didn't produce any useful results,
> > especially estimated scale of almost zero.
> >
> > The problem is that vonmises doesn't have a well defined pdf on the real
line
> >
> >>>> import matplotlib.pyplot as plt
> >>>> plt.plot(scipy.stats.vonmises.pdf(np.linspace(-10,10),1.2, loc=2.3))
> >>>> plt.show()
> >
> > Since vonmises is intended for circular data, and I don't know much
> > about circular statistics (except for what Anne explained on the
> > mailing list), I never tried to get it to work like the other
> > distributions.
> >
> > It might be possible to patch vonmises to work on the real line but I
> > never tried.
>
> Here is a way that seems to work
> * set bounds a,b
> * use mean of random variable as starting value to start inside the
> support of the distribution
>
> >>> vm_rvs = scipy.stats.vonmises.rvs(1.2, 2.3, size=10000)
> >>> stats.vonmises.a = -np.pi
> >>> stats.vonmises.b = np.pi
> >>> stats.vonmises.fit(vm_rvs,1,loc=vm_rvs.mean(),scale=1)
> array([ 1.19767227,  2.30077064,  0.99916372])
>
> Josef
>
>
> >
> > Josef
> >
> >>
> >> --David Ho
> >>
> >> PS: I also don't fully understand the "scale" and "loc" parameters for
all
> >> of the continuous random variables.
> >> Does "loc" always correspond to the mean, and "scale" always correspond
to
> >> the square root of the variance, or something else?
> >
> > loc=0, scale=1 is the standard distribution, the loc and scale
> > parameters transform the distribution of the transformation x~ =
> > (x-loc)/scale
> >
> > If the standard distribution has mean or variance different from 0, 1
> > then the transformed distribution has mean and variance different from
> > loc, scale.
> > the stats method shows the implied variance scale:
> >
> >>>> scipy.stats.vonmises.stats(1.2, loc=2.3, scale=2)
> > (array(2.2999999999999998), array(5.4900668161122423))
> >>>> scipy.stats.vonmises.stats(1.2, loc=0, scale=1)
> > (array(1.5832118030775403e-017), array(1.3725167040280606))
> >
> >
> >>
> >> _______________________________________________
> >> SciPy-User mailing list
> >> SciPy-User at scipy.org
> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>
> >>
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100413/cc458367/attachment.html>