[SciPy-User] How do I use vonmises.fit()?
josef.pktd@gmai...
josef.pktd@gmai...
Tue Apr 13 14:59:32 CDT 2010
On Tue, Apr 13, 2010 at 3:13 PM, <josef.pktd@gmail.com> wrote:
> On Tue, Apr 13, 2010 at 2:52 PM, David Ho <itsdho@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.
No, I don't think this is correct, you are using circular
transformation to shift your original data, and I don't think the
distribution would be able to do this correctly.
The results in the example look good.
All other ideas, I have, might also not be correct because I don't
have any intuition for circular statistics.
Josef
>
> 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@gmail.com> wrote:
>>>
>>> On Tue, Mar 30, 2010 at 4:34 PM, <josef.pktd@gmail.com> wrote:
>>> > On Tue, Mar 30, 2010 at 3:10 PM, David Ho <itsdho@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@scipy.org
>>> >> http://mail.scipy.org/mailman/listinfo/scipy-user
>>> >>
>>> >>
>>> >
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>
More information about the SciPy-User
mailing list