[SciPy-User] How do I use vonmises.fit()?
josef.pktd at gmail.com
josef.pktd at gmail.com
Tue Apr 13 15:13:18 EDT 2010
On Tue, Apr 13, 2010 at 2:52 PM, David Ho <itsdho at ucla.edu> wrote:
> 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. =)
Just a quick reply, I haven't read the script very carefully
loc is shifting the location, scale is extending the support which
will be 2*pi*scale.
in the way I did it originally with starting values
stats.vonmises.fit(vm_rvs,1,loc=vm_rvs.mean(),scale=1)
the estimation is supposed to do the location change automatically. I
think your shifting is essentially the same as what the distributions
do internally with loc.
In the current fit version, loc (and with it the support of the
distribution) and scale are always estimated. In some cases this is
not desired.
You are transforming the original data to fit into the standard
distribution with loc=0, scale=1 . Do you get reasonable estimates for
loc and scale in this case?
If not, then there is another patch or enhanced fit function that
could take loc and scale as fixed.
I will look at some details in your function later, especially I'm
curious how the circular statistics works.
Thanks for the example, maybe I can stop ignoring vonmises.
Josef
>
> --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
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
More information about the SciPy-User
mailing list