how to find solutions with fsolve

[Jose Guzman]

This probably has not much to do with Sage, but I need some help. I am
using the binom package of scipy.stats to solve how many trials (N) I
would need to get 40 or more successes at least in 2% of the times.
Basically what I am doing with sage is:

from scipy.stats import binom
from scipy.stats import fsolve

def mybinom(N):
return 1-binom(N,p).cdf(40)

fsolve(mybinomial(N)-0.02, 100)

however this returns an error:

TypeError: unsupported operand type(s) for *: 'numpy.ndarray' and 'numpy.bool_'

is there anything similar to the command FindRoot in mathematica ???

g[N_] := Interpolation[Table
FindRoot[g[x] == 0.02 , n, {x, 100, 101}, MaxIterations -> 1000]

I tried some of the scalar optimizations of scipy.optimize but I could
not get it working.

Thanks a lot in advance!

[Jason Grout]

On 02/26/2010 04:20 PM, Jose Guzman wrote:
> This probably has not much to do with Sage, but I need some help. I am
> using the binom package of scipy.stats to solve how many trials (N) I
> would need to get 40 or more successes at least in 2% of the times.
> Basically what I am doing with sage is:
>
> from scipy.stats import binom
> from scipy.stats import fsolve
>
> def mybinom(N):
> return 1-binom(N,p).cdf(40)

Here, I get an import error:

sage: from scipy.stats import binom
sage: from scipy.stats import fsolve
Traceback (most recent call last):
...
ImportError: cannot import name fsolve

scipy.stats does not have an fsolve function.

Also, you don't say above what p is, so mybinom doesn't work either.

Could you post a self-contained example that gives the error you get below?

>
> fsolve(mybinomial(N)-0.02, 100)
>
> however this returns an error:
>
> TypeError: unsupported operand type(s) for *: 'numpy.ndarray' and
> 'numpy.bool_'
>
>
> is there anything similar to the command FindRoot in mathematica ???
>
> g[N_] := Interpolation[Table
> FindRoot[g[x] == 0.02 , n, {x, 100, 101}, MaxIterations -> 1000]

It seems like your mathematica code got cut off in the definition of g.

Sage has a find_root command. See
http://www.sagemath.org/doc/reference/sage/numerical/optimize.html#sage.numerical.optimize.find_root

However, I think it depends on the function being continuous. I assume
that's why the mathematica code appears to find the root of an
interpolated function?

Thanks,

Jason

[Jose Guzman]

Dear Jason,

first of all, thank you very much for your kind help and time.

> Could you post a self-contained example that gives the error you get
> below?
>

here my running example. I got to solve it late last night...

>>> from scipy.stats import binom
>>> from scipy.optimize import fsolve

>>> def mybinom(N):
>>>...."""" returns the probability of finding 40 or more successes in
a population of N samples """"
>>>....return 1-binom(N, 0.15).cdf(40)

>>> mybinom(200) # I know the value should be around 200 (solved by
brute-force)
>>> 0.02199927

>>> fsolve(lambda N: mybinom(N)-0.020, 100)
>>> 199.5458489705228

a float value of N=199.54 does not make much sense for a binomial
experiment. I am new here but it seems that fsolve changes depending on
the starting estimation parameter.

>>> int(round(fsolve(lambda N: mybinom(N) - 0.02, 100)))
>>> 200

>>> int(round(fsolve(lambda N: mybinom(N) - 0.02, 90)))
>>> 199

.. is there any optimize better method to get an integer? I did not get
find_root() to work :/

> It seems like your mathematica code got cut off in the definition of g.
>
> Sage has a find_root command. See
> http://www.sagemath.org/doc/reference/sage/numerical/optimize.html#sage.numerical.optimize.find_root
>
>
> However, I think it depends on the function being continuous. I
> assume that's why the mathematica code appears to find the root of an
> interpolated function?
>

Exactly, this is how the mathematica code solve it, but I became a big
fan of Sage!

Thank you very much!

Jose.

[Harald Schilly]

On Feb 26, 11:20 pm, Jose Guzman <n...@neurohost.org> wrote:
> This probably has not much to do with Sage, but I need some help.

I do not fully get your example and I might be wrong, but I have the
feeling you are looking for the inverse survival function, i.e.
binom.isf(...) [i'm only sure that 1-binom.cdf == binom.sf and a
survival function has something to do with reliability ]

h