Subversion scipy.stats irregular problem with source code example
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered. Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like: Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line: Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func and program flow is stopped. In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code. James Phillips ######################################################## import numpy, scipy, scipy.stats # test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data) # test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data) print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print
On Tue, Sep 28, 2010 at 1:12 PM, James Phillips <zunzun@zunzun.com> wrote:
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered.
Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like:
Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract
and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line:
Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func
and program flow is stopped.
In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code.
James Phillips
########################################################
import numpy, scipy, scipy.stats
# test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data)
# test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data)
print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print _______________________________________________
Is there an existing bug ticket for this? If not there probably should be... I think the fitting code should be looked at as experimental. It's good that you caught that no fitting is actually done in these cases. The problem stems (for the most part) from bad starting values (outside the support of the distribution for those with bounded support). I've tried to go through and fix this, giving very naive (but correct) starting values to fit methods, but I haven't gotten much further than that. I don't know if Travis or Josef have gone back to look at this. Hopefully one of these days I will find some more time to look at this and try to give a systematic fix. Skipper
On Tue, Sep 28, 2010 at 7:18 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
I think the fitting code should be looked at as experimental.
Oops. Not knowing this, I naively fit all continuous distributions to a user-entered data set and return the fitted results. Link is: http://zunzun.com/StatisticalDistributions/1/ and push the "Submit" version to fit the example data to all continuous distributions. Takes about 15-20 seconds and shows sorted results with graphical output - this is with the SciPy version in Ubuntu Lucid Lynx, 0.7.0. I was actually trying the repository versions of Numpy and SciPy to see if their fits would work better than in 0.7.0. James
On Tue, Sep 28, 2010 at 1:31 PM, James Phillips <zunzun@zunzun.com> wrote:
On Tue, Sep 28, 2010 at 7:18 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
I think the fitting code should be looked at as experimental.
Just my opinion, so take it for what you will. It should be okay for those distributions that do not have a bounded support, but as you indicated, some do no fitting and just return the start values.
Oops. Not knowing this, I naively fit all continuous distributions to a user-entered data set and return the fitted results. Link is:
http://zunzun.com/StatisticalDistributions/1/
and push the "Submit" version to fit the example data to all continuous distributions. Takes about 15-20 seconds and shows sorted results with graphical output - this is with the SciPy version in Ubuntu Lucid Lynx, 0.7.0.
Very cool. Are the test data sets you use in the public domain?
I was actually trying the repository versions of Numpy and SciPy to see if their fits would work better than in 0.7.0.
James _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Tue, Sep 28, 2010 at 12:42 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
Very cool. Are the test data sets you use in the public domain?
BSD licensed source for all fitting code is at http://code.google.com/p/pythonequations/downloads/list, the made-up fictional example data is surely public domain. James
On 9/28/2010 1:46 PM, James Phillips wrote:
BSD licensed source for all fitting code is at http://code.google.com/p/pythonequations/downloads/list
Is this something that might appropriately be moved into the statsmodels scikit? Alan Isaac
You might consider the rate of change might be too high for that (?) - see http://code.google.com/p/pythonequations/source/list - the source code repository is quite active on updates and commits. James On Tue, Sep 28, 2010 at 12:52 PM, Alan G Isaac <alan.isaac@gmail.com> wrote:
On 9/28/2010 1:46 PM, James Phillips wrote:
BSD licensed source for all fitting code is at http://code.google.com/p/pythonequations/downloads/list
Is this something that might appropriately be moved into the statsmodels scikit?
Alan Isaac
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Tue, Sep 28, 2010 at 2:01 PM, James Phillips <zunzun@zunzun.com> wrote:
You might consider the rate of change might be too high for that (?) - see http://code.google.com/p/pythonequations/source/list - the source code repository is quite active on updates and commits.
If interested, we might be able to mirror your repo(?). Skipper
On Tue, Sep 28, 2010 at 8:04 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Tue, Sep 28, 2010 at 2:01 PM, James Phillips <zunzun@zunzun.com> wrote:
If interested, we might be able to mirror your repo(?).
I would be honored, sir. James
On 9/28/2010 2:01 PM, James Phillips wrote:
You might consider the rate of change might be too high
As long as you can provide unit tests, I don't see a problem. But you and Skipper shd work out the details. Cheers, Alan Isaac
On Tue, Sep 28, 2010 at 2:10 PM, Alan G Isaac <alan.isaac@gmail.com> wrote:
On 9/28/2010 2:01 PM, James Phillips wrote:
You might consider the rate of change might be too high
As long as you can provide unit tests, I don't see a problem.
But you and Skipper shd work out the details.
The part that fits the distributions from what I can see is LGPL License is GNU LGPL v3, see https://launchpad.net/twolumps more in a short time Josef
Cheers, Alan Isaac _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Tue, Sep 28, 2010 at 3:35 PM, <josef.pktd@gmail.com> wrote:
On Tue, Sep 28, 2010 at 2:10 PM, Alan G Isaac <alan.isaac@gmail.com> wrote:
On 9/28/2010 2:01 PM, James Phillips wrote:
You might consider the rate of change might be too high
As long as you can provide unit tests, I don't see a problem.
But you and Skipper shd work out the details.
The part that fits the distributions from what I can see is LGPL
License is GNU LGPL v3, see https://launchpad.net/twolumps
more in a short time
James, I didn't see where in pythonequations you calculate the results statistics, e.g. aic, in the base class? The ranking of distributions, that I have seen are based on kolmogorov-smirnov, chisquare-test or Anderson-Darling. I have a test script where I use kstest for the ranking, and starting values dependent on the support, similar to twolump (but not as clean.) As Skipper said, the fit method works pretty well for the distributions with open support, with bound support setting the starting parameters based the data also works ok for many of them, e.g. like the min in twolump. For several of the bound support distributions an option to fix the support, (e.g. loc) would make more sense. (generalized extreme value and generalized pareto are supposed to have some problems with mle for a large parameter range.) The current scipy trunk is still "experimental", in the sense that it doesn't work for all distributions or might raise some random exceptions. I haven't looked much at the new fit version yet. (I'm still experimenting in statsmodels.) Any bugs with the current code should go into trac tickets. What I did in my test script was to split up the distributions into open, one-sided and bound distributions, and select starting values depending on this. I would restrict to the subset of distributions for which MLE works or put it inside an try except ( I think pymvpa did it this way.) Josef
Josef
Cheers, Alan Isaac _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
I am in the process of eliminating the LGPL code, all I sed it for was to iterate over the distributions and guess initial parameters. It was the investigation of why thiat code would not work correctly that led me to the problem I originally discussed. The "Simple" examples have a StatisticalDistribution.py file, as soon as I dump the LGPL code I'll make a "Complex" example that iterates over all of the continuous distributions. Perhaps we should wait until the LGPL file is gone from my source code distribution before considering the code for mirroring. I do not use the Greedy Programmer's License (GPL) myself. James On Tue, Sep 28, 2010 at 9:35 AM, <josef.pktd@gmail.com> wrote:
On Tue, Sep 28, 2010 at 2:10 PM, Alan G Isaac <alan.isaac@gmail.com> wrote:
On 9/28/2010 2:01 PM, James Phillips wrote:
You might consider the rate of change might be too high
As long as you can provide unit tests, I don't see a problem.
But you and Skipper shd work out the details.
The part that fits the distributions from what I can see is LGPL
License is GNU LGPL v3, see https://launchpad.net/twolumps
more in a short time
Josef
Cheers, Alan Isaac _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Tue, Sep 28, 2010 at 4:28 PM, James Phillips <zunzun@zunzun.com> wrote:
I am in the process of eliminating the LGPL code, all I sed it for was to iterate over the distributions and guess initial parameters. It was the investigation of why thiat code would not work correctly that led me to the problem I originally discussed.
The "Simple" examples have a StatisticalDistribution.py file, as soon as I dump the LGPL code I'll make a "Complex" example that iterates over all of the continuous distributions.
Perhaps we should wait until the LGPL file is gone from my source code distribution before considering the code for mirroring. I do not use the Greedy Programmer's License (GPL) myself.
I just added the beta patch to my old script http://bazaar.launchpad.net/~josef-pktd/statsmodels/statsmodels-josef-experi... This runs with scipy trunk from maybe a month ago. it saves 79 histograms in a sub-directory. I haven't looked at it in a while, and I didn't update it yet with what Skipper found out when he worked on this. I have been working more on getting standard errors for the estimates, using GenericLikelihoodModel which is in development in statsmodel, which should give asymptotic or bootstrap and eventually profile likelihood standard errors and confidence intervals for the parameter estimates, and allows parameters to depend on some explanatory variables. (When I tried to do this this summer, I got stuck in the literature on pareto distributions and power-laws.) One of the main work for fitting is getting good starting values for individual distributions, and I think Skipper went further than what I have. It would be good if we can get this to work (with BSD). Josef
James
On Tue, Sep 28, 2010 at 9:35 AM, <josef.pktd@gmail.com> wrote:
On Tue, Sep 28, 2010 at 2:10 PM, Alan G Isaac <alan.isaac@gmail.com> wrote:
On 9/28/2010 2:01 PM, James Phillips wrote:
You might consider the rate of change might be too high
As long as you can provide unit tests, I don't see a problem.
But you and Skipper shd work out the details.
The part that fits the distributions from what I can see is LGPL
License is GNU LGPL v3, see https://launchpad.net/twolumps
more in a short time
Josef
Cheers, Alan Isaac _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Tue, Sep 28, 2010 at 10:46 AM, <josef.pktd@gmail.com> wrote:
One of the main work for fitting is getting good starting values for individual distributions...
Attached is source code that uses Robert Kern's sandbox version of the Differential Evolution genetic algorithm (diffev.py) to find starting parameters of loc and scale, it works well in my tests. I plan to use this approach on the zunzun.com web site, having no better general options than this one at the moment. Note that I have slightly modified Robert's code so that I can pass in a stopping criteria, that file is also attached. James
On Thu, Sep 30, 2010 at 2:28 PM, James Phillips <zunzun@zunzun.com> wrote:
On Tue, Sep 28, 2010 at 10:46 AM, <josef.pktd@gmail.com> wrote:
One of the main work for fitting is getting good starting values for individual distributions...
Attached is source code that uses Robert Kern's sandbox version of the Differential Evolution genetic algorithm (diffev.py) to find starting parameters of loc and scale, it works well in my tests. I plan to use this approach on the zunzun.com web site, having no better general options than this one at the moment. Note that I have slightly modified Robert's code so that I can pass in a stopping criteria, that file is also attached.
interesting, Why did you choose to minimize the squared difference of quantiles, instead of the negative log-likelihood in the diffev ? Josef
James
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Thu, Sep 30, 2010 at 10:05 AM, <josef.pktd@gmail.com> wrote:
Why did you choose to minimize the squared difference of quantiles, instead of the negative log-likelihood in the diffev ?
For my use in estimating initial parameters, the genetic algorithm needs a cutoff point at which a sufficient solution has been reached and I have useful starting parameters. While I would not know in advance what log likelihood value to use as a stopping parameter, I do know in advance that a residual sum-of-squares near zero should have good initial parameter estimates. That means I can pass a small value to the genetic algorithm and have it stop if it finds such parameters, there would no need to continue past that point. I did notice that it really drags on for the beta distribution, I would need to tune my GA parameters; for example, to give up after fewer than the 500 generations I used in the initial test code. James
On Thu, Sep 30, 2010 at 4:27 PM, James Phillips <zunzun@zunzun.com> wrote:
On Thu, Sep 30, 2010 at 10:05 AM, <josef.pktd@gmail.com> wrote:
Why did you choose to minimize the squared difference of quantiles, instead of the negative log-likelihood in the diffev ?
For my use in estimating initial parameters, the genetic algorithm needs a cutoff point at which a sufficient solution has been reached and I have useful starting parameters. While I would not know in advance what log likelihood value to use as a stopping parameter, I do know in advance that a residual sum-of-squares near zero should have good initial parameter estimates. That means I can pass a small value to the genetic algorithm and have it stop if it finds such parameters, there would no need to continue past that point.
ppf is (very) expensive to calculate for some distributions. I tried something similar, but then switched to matching the cdf instead of the pdf. The differential evolution seems pretty slow in this case. What I also did in larger samples is to match only a few quantiles instead of each observation. (I haven't quite figured out yet how to get standard errors when matching quantiles this way.) I tried a few different starting values for your powerlaw example and my impression is that it doesn't converge to a unique solution, i.e. different starting values end up with different maximum likellihood fit results. I haven't tried powerlaw specifically before, but for pareto there are some problems with mle if loc and scale are also estimated. It could be that powerlaw is also a distribution that requires special fitting methods. I was using gamma as one of the distributions to play with for fitting.
I did notice that it really drags on for the beta distribution, I would need to tune my GA parameters; for example, to give up after fewer than the 500 generations I used in the initial test code.
If you mainly look for starting values instead of getting a good estimate, then it might be possible just to do a rough randomization. Especially for online/interactive estimation of many distributions, this looks a bit too slow to me. But I think using a global optimizer will be quite a bit more robust for several distributions where the likelihood doesn't have a well behaved shape. Josef
James _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Thu, Sep 30, 2010 at 7:20 PM, <josef.pktd@gmail.com> wrote:
ppf is (very) expensive to calculate for some distributions.
and
different starting values end up with different maximum likellihood fit results.
As we here in Alabama sometimes say, "Well, dang it!" So there is no (even slow) general method to use at the moment. For the time being, it seems the zunzun.com web site will have to cut back on the number of distributions. James
On Thu, Sep 30, 2010 at 2:20 PM, <josef.pktd@gmail.com> wrote:
ppf is (very) expensive to calculate for some distributions.
But I think using a global optimizer will be quite a bit more robust for several distributions where the likelihood doesn't have a well behaved shape.
Version Two of Robert Kern's edited code and my example file are attached, this runs much faster. It seems to work considerably better that the first versions I posted. I have folded optimize.fmin() into Robert's code, among several other changes. The example works for powerlaw, beta, gamma, and pareto, although some of the distros are very parameter sensitive as you might see at the end of the example. This code would work much faster if it did not always solve for loc and scale, in some distributions these are known from the data and should not be solved for. Please give this test code a try. James
On Thu, Sep 30, 2010 at 2:20 PM, <josef.pktd@gmail.com> wrote:
ppf is (very) expensive to calculate for some distributions.
I forgot to mention that the code uses the log likelihood nnlf() method now. Thank you for letting me know about the computational expense. James
On Fri, Oct 1, 2010 at 12:11 PM, James Phillips <zunzun@zunzun.com> wrote:
On Thu, Sep 30, 2010 at 2:20 PM, <josef.pktd@gmail.com> wrote:
ppf is (very) expensive to calculate for some distributions.
I forgot to mention that the code uses the log likelihood nnlf() method now. Thank you for letting me know about the computational expense.
Which version of scipy are you using for testing? I had to fix some python 2.5 incompatibilities (*args), but when I ran the script it seemed to get stuck after printing 4 or 5 iterations. I let it run for a few minutes but then killed the process. Josef
James _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Fri, Oct 1, 2010 at 6:55 AM, <josef.pktd@gmail.com> wrote:
Which version of scipy are you using for testing?
Scipy 0.7.0 on Ubuntu Lucid Lynx with Python 2.6.5.
I had to fix some python 2.5 incompatibilities (*args), but when I ran the script it seemed to get stuck after printing 4 or 5 iterations. I let it run for a few minutes but then killed the process.
At the top of diffev.py's solve() method is a loop that runs per generation of the genetic algorithm, you might consider placing a "print gen" statement there to see if it is running or stuck somehow. The virtual server I'm using has 4 CPUs so I plan to run the distribution fitting in parallel. Here are current timing results fitting both loc and scale, in parallel these would have taken 207 seconds (gamma time): distribution: powerlaw took 69.5095369816 seconds to run de parameters [ 1.31332767 0.55933617 2.45766383] nnlf: 23.3610045864 distribution: beta took 47.2591450214 seconds to run de parameters [ 0.70567183 0.6011756 0.607 2.41 ] nnlf: 7.45846316363 distribution: gamma took 206.991358995 seconds to run de parameters [ 1.40727512e+02 -5.83784812e+00 5.55091421e-02] nnlf: 26.7597491247 distribution: pareto took 86.9141609669 seconds to run de parameters [ 7.75987131e+13 -1.05369726e+14 1.05369726e+14] nnlf: 35.2599606822 James
On Fri, Oct 1, 2010 at 7:21 AM, James Phillips <zunzun@zunzun.com> wrote:
Here are current timing results fitting both loc and scale...
I'm iterating over all continuous distributions now, and the genetic algorithm results are showing which distributions can be run with loc = min(data) and scale = max(data) - min(data). With that information in hand I can then speed up the overall fitting considerably by not fitting those parameters. James
On Fri, Oct 1, 2010 at 1:58 PM, James Phillips <zunzun@zunzun.com> wrote:
On Fri, Oct 1, 2010 at 7:21 AM, James Phillips <zunzun@zunzun.com> wrote:
Here are current timing results fitting both loc and scale...
I'm iterating over all continuous distributions now, and the genetic algorithm results are showing which distributions can be run with loc = min(data) and scale = max(data) - min(data). With that information in hand I can then speed up the overall fitting considerably by not fitting those parameters.
good to hear I used loc = min(data) - 1e-6 to avoid evaluating the log likelihood at the lower bound, in some cases this might be infinite. (I picked 1e-6 arbitrarily)
scipy.stats.gamma.logpdf(np.linspace(0,1,3), 300.5609591140931425, loc=0, scale=0.25) array([ -Inf, -1205.37511721, -999.735283 ])
One more: vonmises doesn't define the bounds because it can be used as a circular distribution. If you want to include it as a regular distribution, then you should set the bound a and b. stats.distributions.vonmises.a = -np.pi stats.distributions.vonmises.b = np.pi (Even more fun would be to specify a proper Bayesian prior on all distribution parameters to make the random starting values match the distribution.) Josef
James _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Fri, Oct 1, 2010 at 8:24 AM, <josef.pktd@gmail.com> wrote:
I used loc = min(data) - 1e-6 to avoid evaluating the log likelihood at the lower bound, in some cases this might be infinite. (I picked 1e-6 arbitrarily)
scipy.stats.gamma.logpdf(np.linspace(0,1,3), 300.5609591140931425, loc=0, scale=0.25) array([ -Inf, -1205.37511721, -999.735283 ])
Ah - understood. Thank you.
One more: vonmises doesn't define the bounds because it can be used as a circular distribution. If you want to include it as a regular distribution, then you should set the bound a and b. stats.distributions.vonmises.a = -np.pi stats.distributions.vonmises.b = np.pi
I shall try this. James
The genetic algorithm approach is not working as a general solution to the problem of finding starting parameters for fmin() for statistical distributions, presumably due to extreme parameter sensitivity. I do not see a general solution to the problem given these results. See the attached Python file, also copied below. My results: Digits of precision test for the beta distribution nnlf native = inf nnlf 16 digits = 10.14091764 nnlf 15 digits = 10.3222074111 nnlf 14 digits = 10.977829575 nnlf 13 digits = inf nnlf 12 digits = 13.198954184 James import scipy, scipy.stats data = scipy.array([ 3.017,2.822,2.632,2.287,2.207,2.048, 1.963,1.784,1.712,2.972,2.719,2.495, 2.070,1.969,1.768,1.677,1.479,1.387, 2.843,2.485,2.163,1.687,1.408,1.279, 1.016,0.742,0.607]) # parameters p1 = 7.69589403034175001E-01 p2 = 5.52884409849620395E-01 p3 = 6.06094740472452820E-01 p4 = 2.41090525952754753E+00 print "Digits of precision test for the beta distribution" print "nnlf native =", scipy.stats.beta.nnlf([p1, p2, p3, p3], data) print "nnlf 16 digits =", scipy.stats.beta.nnlf([float("%.16E" % p1), float("%.16E" % p2), float("%.16E" % p3), float("%.16E" % p4)], data) print "nnlf 15 digits =", scipy.stats.beta.nnlf([float("%.15E" % p1), float("%.15E" % p2), float("%.15E" % p3), float("%.15E" % p4)], data) print "nnlf 14 digits =", scipy.stats.beta.nnlf([float("%.14E" % p1), float("%.14E" % p2), float("%.14E" % p3), float("%.14E" % p4)], data) print "nnlf 13 digits =", scipy.stats.beta.nnlf([float("%.13E" % p1), float("%.13E" % p2), float("%.13E" % p3), float("%.13E" % p4)], data) print "nnlf 12 digits =", scipy.stats.beta.nnlf([float("%.12E" % p1), float("%.12E" % p2), float("%.12E" % p3), float("%.12E" % p4)], data) On Fri, Oct 1, 2010 at 7:58 AM, James Phillips <zunzun@zunzun.com> wrote:
On Fri, Oct 1, 2010 at 7:21 AM, James Phillips <zunzun@zunzun.com> wrote:
Here are current timing results fitting both loc and scale...
I'm iterating over all continuous distributions now, and the genetic algorithm results are showing which distributions can be run with loc = min(data) and scale = max(data) - min(data). With that information in hand I can then speed up the overall fitting considerably by not fitting those parameters.
James
On Mon, Oct 11, 2010 at 3:45 PM, James Phillips <zunzun@zunzun.com> wrote:
The genetic algorithm approach is not working as a general solution to the problem of finding starting parameters for fmin() for statistical distributions, presumably due to extreme parameter sensitivity. I do not see a general solution to the problem given these results. See the attached Python file, also copied below.
My results:
Digits of precision test for the beta distribution nnlf native = inf nnlf 16 digits = 10.14091764 nnlf 15 digits = 10.3222074111 nnlf 14 digits = 10.977829575 nnlf 13 digits = inf nnlf 12 digits = 13.198954184
James
import scipy, scipy.stats
data = scipy.array([ 3.017,2.822,2.632,2.287,2.207,2.048, 1.963,1.784,1.712,2.972,2.719,2.495, 2.070,1.969,1.768,1.677,1.479,1.387, 2.843,2.485,2.163,1.687,1.408,1.279, 1.016,0.742,0.607])
# parameters p1 = 7.69589403034175001E-01 p2 = 5.52884409849620395E-01 p3 = 6.06094740472452820E-01 p4 = 2.41090525952754753E+00
print "Digits of precision test for the beta distribution" print "nnlf native =", scipy.stats.beta.nnlf([p1, p2, p3, p3], data)
typo should be p4
print "nnlf 16 digits =", scipy.stats.beta.nnlf([float("%.16E" % p1), float("%.16E" % p2), float("%.16E" % p3), float("%.16E" % p4)], data) print "nnlf 15 digits =", scipy.stats.beta.nnlf([float("%.15E" % p1), float("%.15E" % p2), float("%.15E" % p3), float("%.15E" % p4)], data) print "nnlf 14 digits =", scipy.stats.beta.nnlf([float("%.14E" % p1), float("%.14E" % p2), float("%.14E" % p3), float("%.14E" % p4)], data) print "nnlf 13 digits =", scipy.stats.beta.nnlf([float("%.13E" % p1), float("%.13E" % p2), float("%.13E" % p3), float("%.13E" % p4)], data) print "nnlf 12 digits =", scipy.stats.beta.nnlf([float("%.12E" % p1), float("%.12E" % p2), float("%.12E" % p3), float("%.12E" % p4)], data)
If I remember correctly, you have observations that are too close to the upper boundary.
data.max() -(p3+p4) -4.4408920985006262e-016
If you have an observation at the boundary, the loglikelihood is inf
scipy.stats.beta.nnlf([p1, p2, p3, p4], data) 10.14091764000147
reduce variance a bit, shrinks the support
scipy.stats.beta.nnlf([p1, p2, p3, p4-1e-15], data) 1.#INF
I think in these cases you have to keep the boundary of the support away from the max and min of the data. Similar in other distributions, as I mentioned before. If MLE doesn't work for a distribution then a global optimizer wouldn't help either. In these cases, usually another estimation method is recommended in the literature. For example matching quantiles similar to your initial version. Are all distributions that have support in the entire real line working ? Josef
On Fri, Oct 1, 2010 at 7:58 AM, James Phillips <zunzun@zunzun.com> wrote:
On Fri, Oct 1, 2010 at 7:21 AM, James Phillips <zunzun@zunzun.com> wrote:
Here are current timing results fitting both loc and scale...
I'm iterating over all continuous distributions now, and the genetic algorithm results are showing which distributions can be run with loc = min(data) and scale = max(data) - min(data). With that information in hand I can then speed up the overall fitting considerably by not fitting those parameters.
James
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Mon, Oct 11, 2010 at 10:10 AM, <josef.pktd@gmail.com> wrote:
typo should be p4
Oops - thank you.
If I remember correctly, you have observations that are too close to the upper boundary.
If you have an observation at the boundary, the loglikelihood is inf
I think in these cases you have to keep the boundary of the support away from the max and min of the data. Similar in other distributions, as I mentioned before.
Thank you.
If MLE doesn't work for a distribution then a global optimizer wouldn't help either. In these cases, usually another estimation method is recommended in the literature. For example matching quantiles similar to your initial version.
Thank you again.
Are all distributions that have support in the entire real line working ?
Stopped at beta. I'll pick up with this work again after having read your discussion above. James
On Mon, Oct 11, 2010 at 4:24 PM, James Phillips <zunzun@zunzun.com> wrote:
On Mon, Oct 11, 2010 at 10:10 AM, <josef.pktd@gmail.com> wrote:
typo should be p4
Oops - thank you.
If I remember correctly, you have observations that are too close to the upper boundary.
If you have an observation at the boundary, the loglikelihood is inf
I think in these cases you have to keep the boundary of the support away from the max and min of the data. Similar in other distributions, as I mentioned before.
Thank you.
If MLE doesn't work for a distribution then a global optimizer wouldn't help either. In these cases, usually another estimation method is recommended in the literature. For example matching quantiles similar to your initial version.
Thank you again.
Are all distributions that have support in the entire real line working ?
Stopped at beta. I'll pick up with this work again after having read your discussion above.
I spent 1 to 2 weeks on the pareto family this summer. It's easy to get lost while trying to figure out specific distributions. I would recommend that you concentrate on the "easy" distributions first, and then work your way to the more difficult ones. I hope that eventually I also have enough time again to finish up the other (non-mle) estimators. Josef
James _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Mon, Oct 11, 2010 at 4:24 PM, James Phillips <zunzun@zunzun.com> wrote:
On Mon, Oct 11, 2010 at 10:10 AM, <josef.pktd@gmail.com> wrote:
typo should be p4
Oops - thank you.
If I remember correctly, you have observations that are too close to the upper boundary.
If you have an observation at the boundary, the loglikelihood is inf
I think in these cases you have to keep the boundary of the support away from the max and min of the data. Similar in other distributions, as I mentioned before.
Thank you.
If MLE doesn't work for a distribution then a global optimizer wouldn't help either. In these cases, usually another estimation method is recommended in the literature. For example matching quantiles similar to your initial version.
The maximum product of spacings (MPS) method is a general method of estimating parameters in continuous univariate distributions that in many cases solves this problem. It is especially suited to cases where one of the parameters is an unknown shifted origin. This occurs, for example, in the three-parameter lognormal, gamma, Generalized Extreme Value Generalized Pareto and Weibull models. For such distributions it is known that maximum likelihood (ML) estimation can break down because the likelihood is unbounded and this can lead to inconsistent estimators. In particular MPS is shown to give consistent estimators with asymptotic efficiency equal to ML estimators when these exist. Moreover it gives consistent, asymptoti- cally efficient estimators in situations where ML fails. Finally, as a by-product of the MPS, a goodness of fit statistic, Moran’s statistic, is available for evaluating the fit to the selected distribution. Two years ago I implemented this method + some other enhancements. The source code is available here: http://code.google.com/p/joepython/source/browse/trunk/joepython/scipystats/... The parameters are estimated by minimizing the method nlogps in the rv_continous class. You will find more details on the method in the following references: Estimating Parameters in Continuous Univariate Distributions with a Shifted Origin R. C. H. Cheng and N. A. K. Amin (http://links.jstor.org/sici?sici=0035-9246%281983%2945%3A3%3C394%3AEPICUD%3E...) A Note on the Estimation of Extreme Value Distributions Using Maximum Product of Spacings T. S. T. Wong and W. K. Li Lecture Notes-Monograph Series Vol. 52, Time Series and Related Topics: In Memory of Ching-Zong Wei (2006), pp. 272-283 (article consists of 12 pages) Published by: Institute of Mathematical Statistics Stable URL: http://www.jstor.org/stable/20461444 A Goodness-Of-Fit Test Using Moran's Statistic with Estimated Parameters R. C. H. Cheng and M. A. Stephens (http://links.jstor.org/sici?sici=0006-3444%28198906%2976%3A2%3C385%3AAGTUMS%... )
On Fri, Oct 1, 2010 at 1:21 PM, James Phillips <zunzun@zunzun.com> wrote:
On Fri, Oct 1, 2010 at 6:55 AM, <josef.pktd@gmail.com> wrote:
Which version of scipy are you using for testing?
Scipy 0.7.0 on Ubuntu Lucid Lynx with Python 2.6.5.
I had to fix some python 2.5 incompatibilities (*args), but when I ran the script it seemed to get stuck after printing 4 or 5 iterations. I let it run for a few minutes but then killed the process.
At the top of diffev.py's solve() method is a loop that runs per generation of the genetic algorithm, you might consider placing a "print gen" statement there to see if it is running or stuck somehow.
The virtual server I'm using has 4 CPUs so I plan to run the distribution fitting in parallel. Here are current timing results fitting both loc and scale, in parallel these would have taken 207 seconds (gamma time):
much better than my old notebook
distribution: powerlaw took 69.5095369816 seconds to run de parameters [ 1.31332767 0.55933617 2.45766383] nnlf: 23.3610045864
distribution: beta took 47.2591450214 seconds to run de parameters [ 0.70567183 0.6011756 0.607 2.41 ] nnlf: 7.45846316363
distribution: gamma took 206.991358995 seconds to run de parameters [ 1.40727512e+02 -5.83784812e+00 5.55091421e-02] nnlf: 26.7597491247
distribution: pareto took 86.9141609669 seconds to run de parameters [ 7.75987131e+13 -1.05369726e+14 1.05369726e+14] nnlf: 35.2599606822
I did manage to run gamma, the parameters change a lot from generation to generation, it could also be the case that with your small sample some parameters are not well identified. However, with my scipy version, I'm getting some nans for which I don't see any reason at all.
scipy.stats.gamma.pdf(np.linspace(0,5,11), 300.5609591140931425, loc=0, scale=0.5) array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) scipy.stats.gamma.pdf(np.linspace(0,5,11), 300.5609591140931425, loc=0, scale=0.25) array([ 0., 0., 0., 0., 0., 0., NaN, NaN, NaN, NaN, NaN])
Given that you are doing a randomized search, there might not be enough restrictions that the parameters make sense. I did some fuzz testing when I got started with distributions but haven't run them in a long time. Your random parameters might hit some range where the numerical accuracy and correct results have never been checked. It's possible that fmin could get stuck with some nans in the results. I don't know how many dark corners are left. Josef
James _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Fri, Oct 1, 2010 at 8:07 AM, <josef.pktd@gmail.com> wrote:
On Fri, Oct 1, 2010 at 1:21 PM, James Phillips <zunzun@zunzun.com> wrote:
distribution: gamma took 206.991358995 seconds to run de parameters [ 1.40727512e+02 -5.83784812e+00 5.55091421e-02] nnlf: 26.7597491247
scipy.stats.gamma.pdf(np.linspace(0,5,11), 300.5609591140931425, loc=0, scale=0.5) array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) scipy.stats.gamma.pdf(np.linspace(0,5,11), 300.5609591140931425, loc=0, scale=0.25) array([ 0., 0., 0., 0., 0., 0., NaN, NaN, NaN, NaN, NaN])
Given that you are doing a randomized search, there might not be enough restrictions that the parameters make sense. I did some fuzz testing when I got started with distributions but haven't run them in a long time. Your random parameters might hit some range where the numerical accuracy and correct results have never been checked. It's possible that fmin could get stuck with some nans in the results. I don't know how many dark corners are left.
You have loc=0 where the GA fit the example data set to -5.8 or so, and I see from your results that the NaNs are dependent on your scale value of 0.5 or 0.25. James
Skipper Seabold <jsseabold <at> gmail.com> writes:
On Tue, Sep 28, 2010 at 1:12 PM, James Phillips <zunzun <at> zunzun.com>
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered.
Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like:
Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract
and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line:
Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-
line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-
line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-
wrote: packages/scipy/stats/distributions.py", packages/scipy/stats/distributions.py", packages/scipy/optimize/minpack.py",
line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func
and program flow is stopped.
In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code.
James Phillips
########################################################
import numpy, scipy, scipy.stats
# test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data)
# test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data)
print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print _______________________________________________
Is there an existing bug ticket for this? If not there probably should be...
I think the fitting code should be looked at as experimental. It's good that you caught that no fitting is actually done in these cases. The problem stems (for the most part) from bad starting values (outside the support of the distribution for those with bounded support). I've tried to go through and fix this, giving very naive (but correct) starting values to fit methods, but I haven't gotten much further than that.
I don't know if Travis or Josef have gone back to look at this. Hopefully one of these days I will find some more time to look at this and try to give a systematic fix.
Skipper
Hi, I'm very sorry for entering the thread like this, but after a long search over the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data: ########################## import scipy as sp import scipy.stats as ss exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data) fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) ######################### I then receive many messages on the command line: Warning: invalid value encountered in subtract Which ends with no fitting of the parameters: (1.0, 0.0, 1.0) With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3" that is causing problem with the fitting in this case. I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution? Many thanks, Charles
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Skipper Seabold <jsseabold <at> gmail.com> writes:
On Tue, Sep 28, 2010 at 1:12 PM, James Phillips <zunzun <at> zunzun.com>
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered.
Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like:
Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract
and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line:
Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-
line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-
line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-
wrote: packages/scipy/stats/distributions.py", packages/scipy/stats/distributions.py", packages/scipy/optimize/minpack.py",
line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func
and program flow is stopped.
In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code.
James Phillips
########################################################
import numpy, scipy, scipy.stats
# test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data)
# test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data)
print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print _______________________________________________
Is there an existing bug ticket for this? If not there probably should be...
I think the fitting code should be looked at as experimental. It's good that you caught that no fitting is actually done in these cases. The problem stems (for the most part) from bad starting values (outside the support of the distribution for those with bounded support). I've tried to go through and fix this, giving very naive (but correct) starting values to fit methods, but I haven't gotten much further than that.
I don't know if Travis or Josef have gone back to look at this. Hopefully one of these days I will find some more time to look at this and try to give a systematic fix.
Skipper
Hi, I'm very sorry for entering the thread like this, but after a long search over the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3" that is causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset In [1]: import scipy as sp In [2]: import scipy.stats as ss In [3]: In [4]: exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] In [5]: data = sp.array(exp_data) In [6]: In [7]: fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) In [8]: print(fit_alpha, fit_loc, fit_beta) (0.37079887324711569, 25.599999999999998, 2459.7323873048508) Skipper
Thanks Skipper for the quick answer. I'm using scipy 0.8.0 with pydev and eclipse. So there is a newer version than 0.8.0? 2010/12/9 Skipper Seabold <jsseabold@gmail.com>
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Skipper Seabold <jsseabold <at> gmail.com> writes:
On Tue, Sep 28, 2010 at 1:12 PM, James Phillips <zunzun <at> zunzun.com
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered.
Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like:
Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract
and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line:
Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-
line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-
line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-
wrote: packages/scipy/stats/distributions.py", packages/scipy/stats/distributions.py", packages/scipy/optimize/minpack.py",
line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func
and program flow is stopped.
In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code.
James Phillips
########################################################
import numpy, scipy, scipy.stats
# test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data)
# test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data)
print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print _______________________________________________
Is there an existing bug ticket for this? If not there probably should be...
I think the fitting code should be looked at as experimental. It's good that you caught that no fitting is actually done in these cases. The problem stems (for the most part) from bad starting values (outside the support of the distribution for those with bounded support). I've tried to go through and fix this, giving very naive (but correct) starting values to fit methods, but I haven't gotten much further than that.
I don't know if Travis or Josef have gone back to look at this. Hopefully one of these days I will find some more time to look at this and try to give a systematic fix.
Skipper
Hi, I'm very sorry for entering the thread like this, but after a long search over the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3" that is causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset
In [1]: import scipy as sp
In [2]: import scipy.stats as ss
In [3]:
In [4]: exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]
In [5]: data = sp.array(exp_data)
In [6]:
In [7]: fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data)
In [8]: print(fit_alpha, fit_loc, fit_beta) (0.37079887324711569, 25.599999999999998, 2459.7323873048508)
Skipper _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Thu, Dec 9, 2010 at 3:30 PM, Charles Moliere <charles.moliere@gmail.com> wrote:
Thanks Skipper for the quick answer. I'm using scipy 0.8.0 with pydev and eclipse. So there is a newer version than 0.8.0?
You can get the development version and install it yourself. http://www.scipy.org/Download#head-d0ec9f4cfbf6ed3029a2c409b6d9899586eee0e3 I think 0.9 is to be released fairly soon (?) if you want to wait for an official release. Skipper
On Thu, Dec 9, 2010 at 3:12 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Skipper Seabold <jsseabold <at> gmail.com> writes:
On Tue, Sep 28, 2010 at 1:12 PM, James Phillips <zunzun <at> zunzun.com>
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered.
Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like:
Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract
and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line:
Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-
line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-
line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-
wrote: packages/scipy/stats/distributions.py", packages/scipy/stats/distributions.py", packages/scipy/optimize/minpack.py",
line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func
and program flow is stopped.
In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code.
James Phillips
########################################################
import numpy, scipy, scipy.stats
# test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data)
# test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data)
print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print _______________________________________________
Is there an existing bug ticket for this? If not there probably should be...
I think the fitting code should be looked at as experimental. It's good that you caught that no fitting is actually done in these cases. The problem stems (for the most part) from bad starting values (outside the support of the distribution for those with bounded support). I've tried to go through and fix this, giving very naive (but correct) starting values to fit methods, but I haven't gotten much further than that.
I don't know if Travis or Josef have gone back to look at this. Hopefully one of these days I will find some more time to look at this and try to give a systematic fix.
Skipper
Hi, I'm very sorry for entering the thread like this, but after a long search over the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3" that is causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset
I don't remember any changes, but in this case choosing the right starting values will be important, and the default ones might work with one scipy version but with an other. I think the warnings mean that the starting values don't make much sense and the likelihood is evaluated at "bad" places Josef
In [1]: import scipy as sp
In [2]: import scipy.stats as ss
In [3]:
In [4]: exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]
In [5]: data = sp.array(exp_data)
In [6]:
In [7]: fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data)
In [8]: print(fit_alpha, fit_loc, fit_beta) (0.37079887324711569, 25.599999999999998, 2459.7323873048508)
Skipper _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Thu, Dec 9, 2010 at 3:53 PM, <josef.pktd@gmail.com> wrote:
On Thu, Dec 9, 2010 at 3:12 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Hi, I'm very sorry for entering the thread like this, but after a long search over the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3" that is causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset
I don't remember any changes, but in this case choosing the right starting values will be important, and the default ones might work with one scipy version but with an other.
Ah, well I had an older trunk installed on this machine and got the same as the OP. I updated and then it worked (ie., returned something)... Skipper
On Thu, Dec 9, 2010 at 3:59 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Thu, Dec 9, 2010 at 3:53 PM, <josef.pktd@gmail.com> wrote:
On Thu, Dec 9, 2010 at 3:12 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Hi, I'm very sorry for entering the thread like this, but after a long search over the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3" that is causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset
I don't remember any changes, but in this case choosing the right starting values will be important, and the default ones might work with one scipy version but with an other.
Ah, well I had an older trunk installed on this machine and got the same as the OP. I updated and then it worked (ie., returned something)...
with the version I have:
stats.gamma._fitstart([25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]) (0.56777181998920023, -422.95393359547484, 1742.449869411083)
stats.gamma.fit([25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]) (0.37079853616614639, 25.599999999999998, 2459.7325273317128)
Charles, you could try starting parameters like (0.5, 0.0, 2000) and see if this works. I still haven't looked at all the changes in the fit Josef
Skipper _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
2010/12/9 <josef.pktd@gmail.com>
On Thu, Dec 9, 2010 at 3:53 PM, <josef.pktd@gmail.com> wrote:
On Thu, Dec 9, 2010 at 3:12 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Hi, I'm very sorry for entering the thread like this, but after a long
search over
the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3"
On Thu, Dec 9, 2010 at 3:59 PM, Skipper Seabold <jsseabold@gmail.com> wrote: that is
causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset
I don't remember any changes, but in this case choosing the right starting values will be important, and the default ones might work with one scipy version but with an other.
Ah, well I had an older trunk installed on this machine and got the same as the OP. I updated and then it worked (ie., returned something)...
with the version I have:
stats.gamma._fitstart([25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]) (0.56777181998920023, -422.95393359547484, 1742.449869411083)
stats.gamma.fit([25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]) (0.37079853616614639, 25.599999999999998, 2459.7325273317128)
Charles, you could try starting parameters like (0.5, 0.0, 2000) and see if this works.
I still haven't looked at all the changes in the fit
Josef
Skipper _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
Thanks a lot for your help to you both. I haven't found any documentation about _fitstart(). Is it a new fitting that uses the data to choose the starting parameters (as mentioned in ticket#808 http://projects.scipy.org/scipy/ticket/808)? In this case, is it available in the new trunk or in 0.9.0? I also noticed you are both using 64bits, perhaps it could have an effect on the convergence of the fitting? I apologize for all the questions.
2010/12/9 Charles Moliere <charles.moliere@gmail.com>
2010/12/9 <josef.pktd@gmail.com>
On Thu, Dec 9, 2010 at 3:59 PM, Skipper Seabold <jsseabold@gmail.com>
On Thu, Dec 9, 2010 at 3:53 PM, <josef.pktd@gmail.com> wrote:
On Thu, Dec 9, 2010 at 3:12 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Thu, Dec 9, 2010 at 2:34 PM, Charles <charles.moliere@gmail.com> wrote:
Hi, I'm very sorry for entering the thread like this, but after a long
search over
the web, this thread is the more relevant to my problem which I'm stuck with. I'm actually trying to fit a gamma distribution on a set of experimental values with gamma.fit() in scipy 0.8.0. Here is the very simple code I'm using with a sample of my data:
########################## import scipy as sp import scipy.stats as ss
exp_data =[25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3] data = sp.array(exp_data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data) print(fit_alpha,fit_loc,fit_beta) #########################
I then receive many messages on the command line: Warning: invalid value encountered in subtract
Which ends with no fitting of the parameters: (1.0, 0.0, 1.0)
With earlier version of scipy (0.7.2), the error message are absent but still no fitting is done. Apparently, it is the extrem value of "4500.3"
wrote: that is
causing problem with the fitting in this case.
I know you metionned earlier that the fitting code should be considered as experimental, however I was wondering if this should be considered as a bug, or if I'm making a mistake. In either case, is there a fix for the fit method to work with a gamma distribution?
It looks like Josef's recent changes have got this working. Using the most recent trunk, so you might want to upgrade or see the changeset
I don't remember any changes, but in this case choosing the right starting values will be important, and the default ones might work with one scipy version but with an other.
Ah, well I had an older trunk installed on this machine and got the same as the OP. I updated and then it worked (ie., returned something)...
with the version I have:
stats.gamma._fitstart([25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]) (0.56777181998920023, -422.95393359547484, 1742.449869411083)
stats.gamma.fit([25.6,35.8,100.2,115.2,125.2,140.1,160.6,210.1,250.5,4500.3]) (0.37079853616614639, 25.599999999999998, 2459.7325273317128)
Charles, you could try starting parameters like (0.5, 0.0, 2000) and see if this works.
I still haven't looked at all the changes in the fit
Josef
Skipper _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
Thanks a lot for your help to you both.
I haven't found any documentation about _fitstart(). Is it a new fitting that uses the data to choose the starting parameters (as mentioned in ticket#808 http://projects.scipy.org/scipy/ticket/808)? In this case, is it available in the new trunk or in 0.9.0?
I also noticed you are both using 64bits, perhaps it could have an effect on the convergence of the fitting?
I apologize for all the questions.
Got my answer: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuou... It is in 0.9.x. Then I'll have to wait for the release as i'm not familiar with compiling. Many thanks again!
On Tue, Sep 28, 2010 at 1:12 PM, James Phillips <zunzun@zunzun.com> wrote:
Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered.
Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like:
Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract
and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line:
Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func
same here when I was running my old scripts from scipy 0.7. I opened a ticket a while ago http://projects.scipy.org/scipy/ticket/1276 I'm monkey patching while I figure out what to do stats.distributions.beta_gen._fitstart = lambda self, data : (5,5,0,1) Josef
and program flow is stopped.
In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code.
James Phillips
########################################################
import numpy, scipy, scipy.stats
# test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data)
# test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data)
print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
participants (7)
-
Alan G Isaac -
Charles -
Charles Moliere
-
James Phillips -
josef.pktd＠gmail.com -
Per.Brodtkorb＠ffi.no -
Skipper Seabold