[Tutor] script comparing two values

[Mic Forster]

Thanks Magnus and Jeff,

Magnus, the result I obtained by hand was 0.449. So
you were spot on!


--- Magnus Lycka <magnus@thinkware.se> wrote:
> At Mon, 24 Feb 2003 19:14:57 -0800 (PST), Mic
> Forster wrote:
> >I am attempting to fit an expected geometric series
> of
> >relative species abundance with an observed set of
> >relative species abundance. In doing so one has to
> >determine a constant, k. This is done by iterating
> the
> >following equation:
> >
> >Nmin / N = (k / (1 - k)) * (pow((1 - k),s)) / (1 -
> >(pow((1 - k),s)) = j
> >
> >where Nmin is the number of individuals in the
> least
> >abundant species, N is the total number of
> individuals
> >and s is the number of species. In my example:
> >
> >Nmin = 1
> >N = 862
> >s = 11
> 
> If this was a one time calculation, I might just
> make it
> an interactive python session. First, we need to
> define things.
> Remember as Jeff said to carefully balance
> parenthesis, that
> multiplication never happens implicitly and that **
> is the
> power operator.
> 
> Also remember that Python still thinks that 1 / 2 is
> 0, but
> that 1 / 2.0 is 0.5. In other words, dividing two
> integers will
> give an integer result. (This will change in the
> future.)
> 
>  >>> def f(k, s):
> ...     return k / (1-k) * (1-k)**s / (1-(1-k)**s)
> ...
>  >>> j = 1. / 862
>  >>> s = 11
> 
> Lets just get a rough picture to begin with:
> 
>  >>> for i in range(1,10):
> ...     k = i / 10.
> ...     print k, f(k,s) - j
> ...
> 0.1 0.04965363852
> 0.2 0.0223327647725
> 0.3 0.00748510854082
> 0.4 0.00126736096339
> 0.5 -0.000671573022373
> 0.6 -0.00109717560849
> 0.7 -0.0011559593701
> 0.8 -0.00116001088742
> 0.9 -0.00116009271742
> 
> Obviously, the correct value for k is somewhere
> between 0.4
> and 0.5. I don't know what precision you need. If
> you throw
> these values into excel of some other plotting
> capable program,
> you can see that it's a smooth curve, where the
> derivate is
> constantly negative. This means that you could use
> the bisection
> method, but let's do the manual iteration two more
> times.
> 
>  >>> for i in range(10):
> ...     k = 0.4 + i / 100.
> ...     print k, f(k,s) - j
> ...
> 0.4 0.00126736096339
> 0.41 0.000941824399188
> 0.42 0.000653817213948
> 0.43 0.000399869640732
> 0.44 0.000176721938801
> 0.45 -1.86744412629e-005
> 0.46 -0.000189157813584
> 0.47 -0.000337359101879
> 0.48 -0.000465706337899
> 0.49 -0.000576430579781
>  >>> for i in range(10):
> ...     k = 0.44 + i / 1000.
> ...     print k, f(k,s) - j
> ...
> 0.44 0.000176721938801
> 0.441 0.000155981964688
> 0.442 0.000135516397721
> 0.443 0.000115322326515
> 0.444 9.53968606932e-005
> 0.445 7.57371308685e-005
> 0.446 5.6340288613e-005
> 0.447 3.72035064319e-005
> 0.448 1.8323977733e-005
> 0.449 -3.01083204318e-007
> 
> I just do arrow-up in PythonWin to the for loop,
> press
> enter to get a copy, changes the k = ... row and
> press
> enter again.
> 
> You can continue until you reach your required
> accuracy,
> 0.4489837 or whatever.
> 
> If you want to make this into a normal program, you
> would
> need to make the code determine what interval to
> look closer
> at. In this case, you can use the bisection method.
> Note that
> the bisection method only works well if the functin
> curve is
> nice enough, that is, the derivate doesn't change
> sign in
> the tested inteval. Analysis of the function of a
> simple
> plot will often determine this.
> 
> With bisection, you first try a minimum point and a
> maximum
> point, and make sure that the solutions for these
> points have
> different signs. Then your solution is in between.
> Then you
> try a value in the middle, between min and max. The
> solution
> for this point will have the same sign as min or
> max. Look at
> values from our first run above.
> 
> 0.1 0.04965363852
> 0.5 -0.000671573022373
> 0.9 -0.00116009271742
> 
> This solution for 0.5 has the same sign as for max,
> 0.9. This
> means that the next test should be in the interval
> 0.1 to 0.5.
> So, we should make a bisection function now, right?
> 
>  >>> def bisect(min, max, delta, function):
> ...     fMax = function(max)
> ...     fMin = function(min)
> ...     assert fMax * fMin < 0
> ...     while max - min > delta:
> ...             newX = 0.5 * (min + max)
> ...             fNew = function(newX)
> ...             if fNew * fMax > 0:
> ...                     # fNew has same sign as
> fMax, so newX should 
> replace max
> ...                     max, fMax = newX, fNew
> ...             else:
> ...                     min, fMin = newX, fNew
> ...     return newX, fNew
> ...
>  >>> def fun(x):
> ...     return f(x, s) - j
> ...
> 
> First, let's test the error handling:
> 
>  >>> bisect(0.1, 0.2, delta, fun)
> Traceback (most recent call last):
>    File "<interactive input>", line 1, in ?
>    File "<interactive input>", line 4, in bisect
> AssertionError
> 
> Ok, it won't let us enter an interval where both
> ends
> lead to a solution of the same sign.
> 
> Lets make a very low resolution test.
> 
>  >>> bisect(0.01, 0.99, 0.01, fun)
> (0.3775, 0.0021564007982048661)
> 
> Ok, it's the right region... Let's use a higher
> resolution.
> 
>  >>> delta = 1e-9
>  >>> bisect(delta, 1-delta, delta, fun)
> (0.44898369918536085, 4.9493988790207111e-010)
> 
> How far can we go?
> 
>  >>> delta = 1e-17
>  >>> bisect(delta, 1-delta, delta, fun)
> Traceback (most recent call last):
>    File "<interactive input>", line 1, in ?
>    File "<interactive input>", line 2, in bisect
> 
=== message truncated ===


__________________________________________________
Do you Yahoo!?
Yahoo! Tax Center - forms, calculators, tips, more
http://taxes.yahoo.com/