[SciPy-Dev] optimize.fsolve endless loop with nan

[josef.pktd@gmai...]

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"

>
>>
>> Josef
>>
>>>
>>> Josef