Rosetta Code

Logistic curve fitting in epidemiology: Difference between revisions

Page Discussion

Logistic curve fitting in epidemiology (view source)

Revision as of 03:04, 6 June 2021

15,358 bytes added ,  3 years ago
m
→‎{{header|Phix}}: added syntax colouring the hard way
(Ada version)
m (→‎{{header|Phix}}: added syntax colouring the hard way)
Line 677:
=={{header|Phix}}==
Simplified, my interpretation of shift-cutting (from [[wp:Non-linear_least_squares]])
<!--<lang Phix>(phixonline)-- demo\rosetta\Curve_fit.exw>
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Curve_fit.exw</span>
constant K = 7_800_000_000, -- approx world population
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
n0 = 27, -- number of cases at day 0
<span style="color: #008080;">constant</span> <span style="color: #000000;">K</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">7_800_000_000</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- approx world population</span>
actual = {
<span style="color: #000000;">n0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">27</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- number of cases at day 0</span>
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
<span style="color: #000000;">actual</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span>
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
<span style="color: #000000;">27</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">27</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">27</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">44</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">44</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">60</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">60</span><span style="color: #0000FF;">,</span>
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
<span style="color: #000000;">61</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">61</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">66</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">83</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">219</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">239</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">392</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">534</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">631</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">897</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1350</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2023</span><span style="color: #0000FF;">,</span>
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
<span style="color: #000000;">2820</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4587</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6067</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7823</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9826</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">11946</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">14554</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">17372</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20615</span><span style="color: #0000FF;">,</span>
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
<span style="color: #000000;">24522</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">28273</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">31491</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">34933</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">37552</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">40540</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">43105</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">45177</span><span style="color: #0000FF;">,</span>
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
<span style="color: #000000;">60328</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">64543</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">67103</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">69265</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">71332</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">73327</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">75191</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">75723</span><span style="color: #0000FF;">,</span>
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
<span style="color: #000000;">76719</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">77804</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">78812</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">79339</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">80132</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">80995</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">82101</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">83365</span><span style="color: #0000FF;">,</span>
105824, 109695, 114232, 118610, 125497, 133852, 143227,
<span style="color: #000000;">85203</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">87024</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">89068</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">90664</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">93077</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">95316</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">98172</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">102133</span><span style="color: #0000FF;">,</span>
151367, 167418, 180096, 194836, 213150, 242364, 271106,
<span style="color: #000000;">105824</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">109695</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">114232</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">118610</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">125497</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">133852</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">143227</span><span style="color: #0000FF;">,</span>
305117, 338133, 377918, 416845, 468049, 527767, 591704,
<span style="color: #000000;">151367</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">167418</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">180096</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">194836</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">213150</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">242364</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">271106</span><span style="color: #0000FF;">,</span>
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
<span style="color: #000000;">305117</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">338133</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">377918</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">416845</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">468049</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">527767</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">591704</span><span style="color: #0000FF;">,</span>
1174652 }
<span style="color: #000000;">656866</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">715353</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">777796</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">851308</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">928436</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1000249</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1082054</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">1174652</span> <span style="color: #0000FF;">}</span>
function f(atom r)
atom sq = 0
<span style="color: #008080;">function</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span>
for i=1 to length(actual) do
<span style="color: #004080;">atom</span> <span style="color: #000000;">sq</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
atom eri = exp(r*(i-1)),
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">actual</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
guess = (n0*eri)/(1+n0*(eri-1)/K),
<span style="color: #004080;">atom</span> <span style="color: #000000;">eri</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)),</span>
diff = guess-actual[i]
<span style="color: #000000;">guess</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">n0</span><span style="color: #0000FF;">*</span><span style="color: #000000;">eri</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">n0</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">eri</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">K</span><span style="color: #0000FF;">),</span>
sq += diff*diff
<span style="color: #000000;">diff</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">guess</span><span style="color: #0000FF;">-</span><span style="color: #000000;">actual</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
end for
<span style="color: #000000;">sq</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">diff</span><span style="color: #0000FF;">*</span><span style="color: #000000;">diff</span>
return sq
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">sq</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function solve(integer f, atom guess=0.5, epsilon=0)
atom f0 = f(guess),
<span style="color: #008080;">function</span> <span style="color: #000000;">solve</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">guess</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">epsilon</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
delta = iff(guess?guess:1),
<span style="color: #004080;">atom</span> <span style="color: #000000;">f0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">),</span>
factor = 2 -- double until f0 best, then
<span style="color: #000000;">delta</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">?</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">:</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span>
-- halve until delta<=epsilon
<span style="color: #000000;">factor</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> <span style="color: #000080;font-style:italic;">-- double until f0 best, then
-- or we hit precision limit
while delta>epsilon -- (predefinedhalve limit)until delta&lt;=epsilon
and guess!=guess-delta do -- ( or we hit precision limit)</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">delta</span><span style="color: #0000FF;">></span><span style="color: #000000;">epsilon</span> <span style="color: #000080;font-style:italic;">-- (predefined limit)</span>
atom nf = f(guess-delta)
<span style="color: #008080;">and</span> <span style="color: #000000;">guess</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">-</span><span style="color: #000000;">delta</span> <span style="color: #008080;">do</span> <span style="color: #000080;font-style:italic;">-- (precision limit)</span>
if nf<f0 then
<span style="color: #004080;">atom</span> <span style="color: #000000;">nf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">-</span><span style="color: #000000;">delta</span><span style="color: #0000FF;">)</span>
f0 = nf
<span style="color: #008080;">if</span> <span style="color: #000000;">nf</span><span style="color: #0000FF;"><</span><span style="color: #000000;">f0</span> <span style="color: #008080;">then</span>
guess -= delta
<span style="color: #000000;">f0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nf</span>
else
<span style="color: #000000;">guess</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">delta</span>
nf = f(guess+delta)
<span style="color: if nf#008080;">else<f0 then/span>
<span style="color: #000000;">nf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">+</span><span style="color: #000000;">delta</span><span style="color: #0000FF;">)</span>
f0 = nf
<span style="color: #008080;">if</span> <span style="color: #000000;">nf</span><span style="color: #0000FF;"><</span><span style="color: #000000;">f0</span> <span style="color: #008080;">then</span>
guess += delta
<span style="color: #000000;">f0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nf</span>
else
<span style="color: #000000;">guess</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">delta</span>
factor = 0.5
end if<span style="color: #008080;">else</span>
<span style="color: #000000;">factor</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.5</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
delta *= factor
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end while
<span style="color: #000000;">delta</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">factor</span>
return guess
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">guess</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
atom r = solve(f),
R0 = exp(12 * r)
<span style="color: #004080;">atom</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">solve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">),</span>
printf(1,"r = %.10f, R0 = %.8f\n",{r,R0})</lang>
<span style="color: #000000;">R0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">12</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"r = %.10f, R0 = %.8f\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R0</span><span style="color: #0000FF;">})</span>
<!--</lang>-->
{{out}}
<pre>
Petelomax
7,806

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/299665"