[Tutor] script comparing two values
Magnus Lycka
magnus@thinkware.se
Tue Feb 25 05:11:00 2003
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
File "<interactive input>", line 2, in fun
File "<interactive input>", line 2, in f
ZeroDivisionError: float division
Too far!
>>> delta = 1e-16
>>> bisect(delta, 1-delta, delta, fun)
(0.44898372593474889, -3.6212352561015848e-017)
Is this accurate enough? ;)
To determine how accurate you are, I guess you should
return both min and max. Then you can say that k must
be in that interval.
>>> def bisect(min, max, delta, function):
... fMax = function(max)
... fMin = function(min)
... assert fMax * fMin < 0
... while abs(fMax - fMin) > 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 "%.17f +/- %.17f" % (0.5*(min+max), 0.5*(max-min))
...
>>> bisect(delta, 1-delta, delta, fun)
'0.44898372593474711 +/- 0.00000000000000178'
--
Magnus Lycka, Thinkware AB
Alvans vag 99, SE-907 50 UMEA, SWEDEN
phone: int+46 70 582 80 65, fax: int+46 70 612 80 65
http://www.thinkware.se/ mailto:magnus@thinkware.se