[SciPy-Dev] optimize.fsolve endless loop with nan
josef.pktd@gmai...
josef.pktd@gmai...
Mon Mar 18 10:29:17 CDT 2013
On Mon, Mar 18, 2013 at 11:14 AM, <josef.pktd@gmail.com> wrote:
> On Wed, Mar 13, 2013 at 1:24 PM, <josef.pktd@gmail.com> wrote:
>> On Wed, Mar 13, 2013 at 12:59 PM, <josef.pktd@gmail.com> wrote:
>>> On Wed, Mar 13, 2013 at 12:17 PM, <josef.pktd@gmail.com> wrote:
>>>> On Wed, Mar 13, 2013 at 10:51 AM, <josef.pktd@gmail.com> wrote:
>>>>> preliminary question, I didn't have time yet to look closely
>>>>>
>>>>>>>> scipy.__version__
>>>>> '0.9.0'
>>>>>
>>>>> I have a problem where fsolve goes into a range where the values are
>>>>> nan. After that it goes into an endless loop, as far as I can tell.
>>>>>
>>>>> Something like this has been fixed for optimize.fmin_bfgs. Was there a
>>>>> fix for this also for fsolve, since 0.9.0?
>>>>>
>>>>>
>>>>> (The weirder story: I rearranged some test, and made unfortunately
>>>>> also some other changes, and now when I run nosetests it never
>>>>> returns. Ctrl+C kills nosetests, but leaves a python process running.
>>>>> I have no clue why the test sequence should matter.)
>>>>
>>>> I had left the starting values in a module global even after I started
>>>> to adjust them in one of the cases.
>>>>
>>>> The starting value for fsolve was in a range where the curvature is
>>>> very flat, and fsolve made huge steps into the invalid range. After
>>>> getting nans, it went AWOL.
>>>>
>>>> If I return np.inf as soon as I get a nan, then fsolve seems to stop
>>>> right away. Is there a way to induce fsolve to stay out of the nan
>>>> zone, for example returning something else than inf?
>>>>
>>>> I don't want to find a very specific solution, because I'm throwing
>>>> lot's of different cases at the same generic method.
>>>
>>> same result with python 2.7, scipy version 0.11.0b1
>>>
>>>>"C:\Programs\Python27\python.exe" fsolve_endless_nan.py
>>> scipy version 0.11.0b1
>>> args 0.3 [100] 0.05
>>> args 0.3 [ 100.] 0.05
>>> args 0.3 [ 100.] 0.05
>>> args 0.3 [ 100.00000149] 0.05
>>> args 0.3 [-132.75434239] 0.05
>>> fsolve_endless_nan.py:36: RuntimeWarning: invalid value encountered in sqrt
>>> pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
>>> fsolve_endless_nan.py:39: RuntimeWarning: invalid value encountered in sqrt
>>> pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
>>> args 0.3 [ nan] 0.05
>>>
>>>
>>> standalone test case (from my power branch)
>>>
>>> Don't run in an interpreter (session) that you want to keep alive!
>>> And open TaskManager if you are on Windows :)
>>>
>>> ------------
>>> # -*- coding: utf-8 -*-
>>> """Warning: endless loop in runaway process, requires hard kill of process
>>>
>>> Created on Wed Mar 13 12:44:15 2013
>>>
>>> Author: Josef Perktold
>>> """
>>>
>>> import numpy as np
>>> from scipy import stats, optimize
>>>
>>> import scipy
>>> print "scipy version", scipy.__version__
>>>
>>>
>>> def ttest_power(effect_size, nobs, alpha, df=None, alternative='two-sided'):
>>> '''Calculate power of a ttest
>>> '''
>>> print 'args', effect_size, nobs, alpha
>>> d = effect_size
>>> if df is None:
>>> df = nobs - 1
>>>
>>> if alternative in ['two-sided', '2s']:
>>> alpha_ = alpha / 2. #no inplace changes, doesn't work
>>> elif alternative in ['smaller', 'larger']:
>>> alpha_ = alpha
>>> else:
>>> raise ValueError("alternative has to be 'two-sided', 'larger' " +
>>> "or 'smaller'")
>>>
>>> pow_ = 0
>>> if alternative in ['two-sided', '2s', 'larger']:
>>> crit_upp = stats.t.isf(alpha_, df)
>>> # use private methods, generic methods return nan with negative d
>>> pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
>>> if alternative in ['two-sided', '2s', 'smaller']:
>>> crit_low = stats.t.ppf(alpha_, df)
>>> pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
>>> return pow_
>>>
>>> func = lambda nobs, *args: ttest_power(args[0], nobs, args[1])
>>>
>>> print optimize.fsolve(func, 100, args=(0.3, 0.05))
>>> ------------
>>
>> correction to get a solution that would make sense (last two lines)
>>
>> func = lambda nobs, *args: ttest_power(args[0], nobs, args[1]) - args[2]
>>
>> print optimize.fsolve(func, 10, args=(0.76638635, 0.1, 0.8))
>> #converges to 12
>> print optimize.fsolve(func, 100, args=(0.76638635, 0.1, 0.8)) #runaway
>>
>> Josef
>> "Lost in Translation"
>
> (continuing my monologue)
>
> more problems with fsolve, looks like on Ubuntu
>
> Skipper uses scipy master, but python-xy daily testing uses
> python-scipy (0.10.1+dfsg1-4)
> TravisCI, which also runs Ubuntu doesn't have any problems (using
> python-scipy amd64 0.9.0+dfsg1-1ubuntu2 )
> https://travis-ci.org/statsmodels/statsmodels/jobs/5585028
>
> (just realized: I've tested on Windows so far only with scipy 0.9)
I should have done the testing across versions:
I'm getting the same failures on Windows with
>>> import scipy
>>> scipy.__version__
'0.11.0b1'
but not with scipy 0.9
something has changed in the recent scipy versions that makes fsolve
more fragile
Josef
>
> from Skipper's run ( https://github.com/statsmodels/statsmodels/issues/710 ):
>
> -----------
> [~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
> [67]: optimize.fsolve(func, .14, full_output=True)
> [67]:
> (array([ 0.05]),
> {'fjac': array([[-1.]]),
> 'fvec': array([ 0.]),
> 'nfev': 12,
> 'qtf': array([ 0.]),
> 'r': array([-3.408293])},
> 1,
> 'The solution converged.')
>
> [~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
> [68]: optimize.fsolve(func, .15, full_output=True)
> [68]:
> (array([ 0.15]),
> {'fjac': array([[-1.]]),
> 'fvec': array([ 0.20178728]),
> 'nfev': 4,
> 'qtf': array([-0.20178728]),
> 'r': array([ inf])},
> 1,
> 'The solution converged.')
>
> It looks like the QR factorization is failing (r == inf) and then it's
> reporting convergence still.
> ---------
>
> That it stops with a "solution converged" also doesn't trigger my
> backup root-finding.
>
>
> I can work around this since I have no idea how to debug this.
> However, there might be something fishy with fsolve.
>
> I never had problems with fsolve before, and scipy.stats.distributions
> is a heavy user of it.
>
> Josef
>
>>
>>>
>>>>
>>>> Josef
>>>>
>>>>>
>>>>> Josef
More information about the SciPy-Dev
mailing list