Logistic curve fitting in epidemiology: Difference between revisions
(Added C++ solution) |
m (→{{header|Wren}}: Changed to Wren S/H) |
||
(24 intermediate revisions by 16 users not shown) | |||
Line 1: | Line 1: | ||
{{task|Probability and statistics}} |
|||
{{draft task}} |
|||
The least-squares method (see references below) in statistics is used to fit data |
The least-squares method (see references below) in statistics is used to fit data |
||
Line 12: | Line 12: | ||
including the growth of epidemics, are curves akin to the logistic curve: |
including the growth of epidemics, are curves akin to the logistic curve: |
||
<big><big><big><b> ƒ(x) = L / (1 + e<sup>-k(x-x<sub>0</sub>)</sup>) </b></big></big></big> |
|||
Though predictions based on fitting to such curves may |
Though predictions based on fitting to such curves may error, especially if used to |
||
extrapolate from incomplete data, curves similar to the logistic curve have had |
extrapolate from incomplete data, curves similar to the logistic curve have had |
||
good fits in population studies, including modeling the growth of past epidemics. |
good fits in population studies, including modeling the growth of past epidemics. |
||
The task: |
|||
;Task |
|||
* Given the following daily world totals since December 31, 2019 for persons |
|||
who have become infected with the novel coronavirus Covid-19: |
|||
Given the following daily world totals since <tt> December 31, 2019 </tt> for persons who have become infected with the novel coronavirus Covid-19: |
|||
<pre> |
<pre> |
||
Daily totals: |
Daily totals: |
||
Line 38: | Line 34: | ||
</pre> |
</pre> |
||
* Use the following variant of the logistic curve as a formula: |
|||
Use the following variant of the logistic curve as a formula: |
|||
'''f(t) = n<sub>0</sub> e<sup>(r t)</sup> / ((1 + n<sub>0</sub> (e<sup>(r t)</sup> - 1) / K)''' |
|||
<big><big><big><b> ƒ(t) = n<sub>0</sub> e<sup>(r t)</sup> / ((1 + n<sub>0</sub> (e<sup>(r t)</sup> - 1) / K) </b></big></big></big> |
|||
Where: |
|||
Where r is the rate of growth of the infection in the population. |
|||
::* <big><b> r </b></big> is the rate of growth of the infection in the population. |
|||
::* <big><b> K </b></big> is the world population, about 7.8 billion. |
|||
::* <big><b> n<sub>0</sub> </b></big> is 27, the number of cases found in China at the start of the pandemic. |
|||
The R0 of an infection (different from r above) is a measure of how many |
|||
The <big><b> R0 </b></big> of an infection (different from <big><b>r</b></big> above) is a measure of how many |
|||
new individuals will become infected for every individual currently infected. |
new individuals will become infected for every individual currently infected. |
||
It is an important measure of how quickly an infectious disease may spread. |
It is an important measure of how quickly an infectious disease may spread. |
||
R0 is related to the logistic curve's r parameter by the formula |
<big><b>R0</b></big> is related to the logistic curve's <big><b>r</b></big> parameter by the formula: |
||
<big><big><big><b> r ≈ ln(R0) / G </b></big></big></big> |
|||
where G |
where <big><b> G </b></big> the generation time, is roughly the sum of the incubation time, perhaps 5 days, and |
||
the mean contagion period, perhaps 7 days, so, for covid-19, roughly we have |
the mean contagion period, perhaps 7 days, so, for covid-19, roughly we have: |
||
<big><big><big><b> R0 ≈ e<sup>12r</sup> </b></big></big></big> |
|||
'''R0 ≈ e<sup>12r</sup>''' |
|||
Assume the following constants hold in the formula above: |
|||
* K is the world population, about 7.8 billion |
|||
* n0 is 27, the number of cases found in China at the start of the pandemic. |
|||
;Task: |
;Task: |
||
* Demonstrate code that finds a least-squares fits of the curve to the data. |
* Demonstrate code that finds a least-squares fits of the curve to the data. |
||
* Show the calculated r for the logistic curve. |
* Show the calculated <big><b> r </b></big> for the logistic curve. |
||
* Show the final R0 parameter you calculate from the logistic curve r value parameter. |
* Show the final <big><b> R0 </b></big> parameter you calculate from the logistic curve <big><b> r </b></big> value parameter. |
||
;See also |
|||
;See also |
|||
:;*[[https://en.wikipedia.org/wiki/Basic_reproduction_number Basic reproduction number]] |
:;*[[https://en.wikipedia.org/wiki/Basic_reproduction_number Basic reproduction number]] |
||
:;*[[https://en.wikipedia.org/wiki/Logistic_function Logistic functions]] |
:;*[[https://en.wikipedia.org/wiki/Logistic_function Logistic functions]] |
||
Line 73: | Line 72: | ||
:;*[[https://ourworldindata.org/coronavirus#all-charts-preview World covid-19 case tallies]] |
:;*[[https://ourworldindata.org/coronavirus#all-charts-preview World covid-19 case tallies]] |
||
:;*[[https://www.zoology.ubc.ca/~bio310/Blank%20Page%204_files/DL%20R%20and%20r.htm Calculating r and R0]] |
:;*[[https://www.zoology.ubc.ca/~bio310/Blank%20Page%204_files/DL%20R%20and%20r.htm Calculating r and R0]] |
||
<br><br> |
|||
=={{header|11l}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="11l">-V K = 7.8e9 |
|||
-V n0 = 27 |
|||
-V actual = [ |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, |
|||
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, |
|||
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, |
|||
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, |
|||
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, |
|||
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, |
|||
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, |
|||
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 |
|||
] |
|||
F f(Float r) -> Float |
|||
V sq = 0.0 |
|||
L(act) :actual |
|||
V eri = exp(r * L.index) |
|||
V guess = (:n0 * eri) / (1 + :n0 * (eri - 1) / :K) |
|||
V diff = guess - act |
|||
sq += diff * diff |
|||
R sq |
|||
F solve(func, =guess = 0.5, epsilon = 0.0) |
|||
V delta = I guess != 0 {guess} E 1 |
|||
V f0 = func(guess) |
|||
V factor = 2.0 |
|||
L delta > epsilon & guess != guess - delta |
|||
V nf = func(guess - delta) |
|||
I nf < f0 |
|||
f0 = nf |
|||
guess -= delta |
|||
E |
|||
nf = func(guess + delta) |
|||
I nf < f0 |
|||
f0 = nf |
|||
guess += delta |
|||
E |
|||
factor = 0.5 |
|||
delta *= factor |
|||
R guess |
|||
V r = solve(f) |
|||
V R0 = exp(12 * r) |
|||
print(‘r = #.6, R0 = #.6’.format(r, R0))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
r = 0.112302, R0 = 3.848279 |
|||
</pre> |
|||
=={{header|Ada}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="ada">with Ada.Text_Io; |
|||
with Ada.Numerics.Generic_Elementary_Functions; |
|||
procedure Curve_Fitting is |
|||
type Real is new Long_Float; |
|||
type Fuction_Access is access function (X : Real) return Real; |
|||
type Time_Type is new Natural; -- Days |
|||
package Real_Io is new Ada.Text_Io.Float_Io (Real); |
|||
package Math is new Ada.Numerics.Generic_Elementary_Functions (Real); |
|||
Actual : constant array (Time_Type range <>) of Integer := |
|||
(27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, |
|||
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, |
|||
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, |
|||
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, |
|||
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, |
|||
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, |
|||
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, |
|||
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652); |
|||
N0 : constant Real := Real (Actual (Actual'First)); -- Initially infected |
|||
K : constant Real := 7.8e9; -- World population apx. |
|||
function Logistic_Curve_Function (R : Real) return Real is |
|||
sq : Real := 0.0; |
|||
begin |
|||
for I in Actual'range loop |
|||
declare |
|||
Eri : constant Real := Math.Exp (R * Real (I)); |
|||
Guess : constant Real := (N0 * Eri) / (1.0 + N0 * (Eri - 1.0) / K); |
|||
Diff : constant Real := Guess - Real (Actual (I)); |
|||
begin |
|||
Sq := Sq + Diff * Diff; |
|||
end; |
|||
end loop; |
|||
return sq; |
|||
end Logistic_Curve_Function; |
|||
procedure Solve (Func : Fuction_Access; |
|||
Guess : Real := 0.5; |
|||
Epsilon : Real := 0.0; |
|||
New_Guess : out Real) |
|||
is |
|||
Delt : Real := (if guess /= 0.0 then guess else 1.0); |
|||
f0 : Real := Func (Guess); |
|||
factor : Real := 2.0; |
|||
begin |
|||
New_Guess := Guess; |
|||
while |
|||
Delt > Epsilon and New_Guess /= New_Guess - Delt |
|||
loop |
|||
declare |
|||
nf : real := func (New_guess - delt); |
|||
begin |
|||
if nf < f0 then |
|||
f0 := nf; |
|||
New_guess := New_guess - delt; |
|||
else |
|||
nf := func (New_guess + delt); |
|||
if nf < f0 then |
|||
f0 := nf; |
|||
New_guess := New_guess + delt; |
|||
else |
|||
Factor := 0.5; |
|||
end if; |
|||
end if; |
|||
end; |
|||
Delt := Delt * factor; |
|||
end loop; |
|||
end Solve; |
|||
use Ada.Text_Io; |
|||
Incubation_Time : constant Real := 5.0; -- Days |
|||
Contagion_Period : constant Real := 7.0; -- Days |
|||
Generation_Time : constant Real := Incubation_Time + Contagion_Period; |
|||
R : Real; -- Rate |
|||
R0 : Real; -- Infection rate |
|||
begin |
|||
Solve (Logistic_Curve_Function'Access, New_Guess => R); |
|||
R0 := Math.Exp (Generation_Time * R); |
|||
Put ("r = "); Real_IO.Put (R, Exp => 0, Aft => 7); |
|||
Put (", R0 = "); Real_Io.Put (R0, Exp => 0, Aft => 7); |
|||
New_Line; |
|||
end Curve_Fitting;</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.1123022, R0 = 3.8482793</pre> |
|||
=={{header|Arturo}}== |
|||
<syntaxhighlight lang="arturo">K: 7.8e9 |
|||
N0: 27 |
|||
Actual: [ |
|||
27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0, |
|||
59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0, |
|||
60.0, 61.0, 61.0, 66.0, 83.0, 219.0, 239.0, |
|||
392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0, |
|||
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0, |
|||
20615.0, 24522.0, 28273.0, 31491.0, 34933.0, 37552.0, 40540.0, |
|||
43105.0, 45177.0, 60328.0, 64543.0, 67103.0, 69265.0, 71332.0, |
|||
73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0, |
|||
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0, |
|||
90664.0, 93077.0, 95316.0, 98172.0, 102133.0, 105824.0, 109695.0, |
|||
114232.0, 118610.0, 125497.0, 133852.0, 143227.0, 151367.0, 167418.0, |
|||
180096.0, 194836.0, 213150.0, 242364.0, 271106.0, 305117.0, 338133.0, |
|||
377918.0, 416845.0, 468049.0, 527767.0, 591704.0, 656866.0, 715353.0, |
|||
777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0 |
|||
] |
|||
f: function [r][ |
|||
result: 0 |
|||
loop 0..dec size Actual 'i [ |
|||
eri: exp r * to :floating i |
|||
guess: (N0 * eri) // (1 + N0 * (eri - 1) // K) |
|||
diff: guess - Actual\[i] |
|||
result: result + diff * diff |
|||
] |
|||
return result |
|||
] |
|||
solve: function [fn][ |
|||
guess: 0.5 |
|||
eps: 0.0 |
|||
result: guess |
|||
delta: guess |
|||
f0: call fn @[result] |
|||
factor: 2.0 |
|||
while [and? -> delta > eps |
|||
-> result <> result - delta][ |
|||
nf: call fn @[result-delta] |
|||
if? nf < f0 [ |
|||
f0: nf |
|||
result: result - delta |
|||
] |
|||
else [ |
|||
nf: call fn @[result+delta] |
|||
if? nf < f0 [ |
|||
f0: nf |
|||
result: result + delta |
|||
] |
|||
else [ |
|||
factor: 0.5 |
|||
] |
|||
] |
|||
delta: delta * factor |
|||
] |
|||
return result |
|||
] |
|||
r: solve 'f |
|||
r0: exp 12 * r |
|||
print ["r =" r] |
|||
print ["R0 =" r0]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.11230217569755 |
|||
R0 = 3.848279280916719</pre> |
|||
=={{header|C}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="c">#include <math.h> |
|||
#include <stdio.h> |
|||
const double K = 7.8e9; |
|||
const int n0 = 27; |
|||
const double actual[] = { |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, |
|||
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, |
|||
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, |
|||
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, |
|||
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, |
|||
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, |
|||
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, |
|||
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 |
|||
}; |
|||
const size_t actual_size = sizeof(actual) / sizeof(double); |
|||
double f(double r) { |
|||
double sq = 0; |
|||
size_t i; |
|||
for (i = 0; i < actual_size; ++i) { |
|||
double eri = exp(r * i); |
|||
double guess = (n0 * eri) / (1 + n0 * (eri - 1) / K); |
|||
double diff = guess - actual[i]; |
|||
sq += diff * diff; |
|||
} |
|||
return sq; |
|||
} |
|||
double solve(double (*fn)(double), double guess, double epsilon) { |
|||
double delta, f0, factor; |
|||
for (delta = guess ? guess : 1, f0 = fn(guess), factor = 2; |
|||
delta > epsilon && guess != guess - delta; |
|||
delta *= factor) { |
|||
double nf = (*fn)(guess - delta); |
|||
if (nf < f0) { |
|||
f0 = nf; |
|||
guess -= delta; |
|||
} else { |
|||
nf = fn(guess + delta); |
|||
if (nf < f0) { |
|||
f0 = nf; |
|||
guess += delta; |
|||
} else { |
|||
factor = 0.5; |
|||
} |
|||
} |
|||
} |
|||
return guess; |
|||
} |
|||
double solve_default(double (*fn)(double)) { |
|||
return solve(fn, 0.5, 0); |
|||
} |
|||
int main() { |
|||
double r = solve_default(f); |
|||
double R0 = exp(12 * r); |
|||
printf("r = %f, R0 = %f\n", r, R0); |
|||
return 0; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.112302, R0 = 3.848279</pre> |
|||
=={{header|C++}}== |
=={{header|C++}}== |
||
{{trans|Phix}} |
{{trans|Phix}} |
||
< |
<syntaxhighlight lang="cpp">#include <cmath> |
||
#include <functional> |
#include <functional> |
||
#include <iostream> |
#include <iostream> |
||
Line 98: | Line 389: | ||
double f(double r) { |
double f(double r) { |
||
double sq = 0; |
double sq = 0; |
||
size_t len = |
constexpr size_t len = std::size(actual); |
||
for (size_t i = 0; i < len; ++i) { |
for (size_t i = 0; i < len; ++i) { |
||
double eri = std::exp(r * i); |
double eri = std::exp(r * i); |
||
Line 109: | Line 400: | ||
double solve(std::function<double(double)> fn, double guess=0.5, double epsilon=0) { |
double solve(std::function<double(double)> fn, double guess=0.5, double epsilon=0) { |
||
for (double delta = guess, f0 = fn(guess), factor = 2; |
for (double delta = guess ? guess : 1, f0 = fn(guess), factor = 2; |
||
delta > epsilon && guess != guess - delta; |
delta > epsilon && guess != guess - delta; |
||
delta *= factor) { |
delta *= factor) { |
||
Line 133: | Line 424: | ||
std::cout << "r = " << r << ", R0 = " << R0 << '\n'; |
std::cout << "r = " << r << ", R0 = " << R0 << '\n'; |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 139: | Line 430: | ||
r = 0.112302, R0 = 3.84828 |
r = 0.112302, R0 = 3.84828 |
||
</pre> |
</pre> |
||
=={{header|D}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="d">import std.math; |
|||
import std.stdio; |
|||
immutable K = 7.8e9; |
|||
immutable n0 = 27; |
|||
immutable actual = [ |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, |
|||
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, |
|||
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, |
|||
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, |
|||
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, |
|||
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, |
|||
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, |
|||
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 |
|||
]; |
|||
double f(double r) { |
|||
double sq = 0.0; |
|||
auto len = actual.length; |
|||
for (size_t i = 0; i < len; ++i) { |
|||
auto eri = exp(r * i); |
|||
auto guess = (n0 * eri) / (1 + n0 * (eri - 1) / K); |
|||
auto diff = guess - actual[i]; |
|||
sq += diff * diff; |
|||
} |
|||
return sq; |
|||
} |
|||
double solve(double function(double) fn, double guess = 0.5, double epsilon = 0) { |
|||
for (double delta = guess ? guess : 1, f0 = fn(guess), factor = 2; |
|||
delta > epsilon && guess != guess - delta; |
|||
delta *= factor |
|||
) { |
|||
double nf = fn(guess - delta); |
|||
if (nf < f0) { |
|||
f0 = nf; |
|||
guess -= delta; |
|||
} else { |
|||
nf = fn(guess + delta); |
|||
if (nf < f0) { |
|||
f0 = nf; |
|||
guess += delta; |
|||
} else { |
|||
factor = 0.5; |
|||
} |
|||
} |
|||
} |
|||
return guess; |
|||
} |
|||
void main() { |
|||
double r = solve(&f); |
|||
double R0 = exp(12 * r); |
|||
writeln("r = ", r, ", R0 = ", R0); |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.112302, R0 = 3.84828</pre> |
|||
=={{header|FreeBASIC}}== |
|||
{{trans|Phix}} |
|||
<syntaxhighlight lang="freebasic">Const K = 7.8e9 ' approx world population |
|||
Const n0 = 27 ' number of cases at day 0 |
|||
Dim Shared As Integer actual(1 To ...) => { _ |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, _ |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, _ |
|||
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, _ |
|||
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177, _ |
|||
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723, _ |
|||
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, _ |
|||
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133, _ |
|||
105824, 109695, 114232, 118610, 125497, 133852, 143227, _ |
|||
151367, 167418, 180096, 194836, 213150, 242364, 271106, _ |
|||
305117, 338133, 377918, 416845, 468049, 527767, 591704, _ |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 } |
|||
Function f1(r As Double) As Double |
|||
Dim As Double sq = 0 |
|||
For i As Integer = 1 To Ubound(actual) |
|||
Dim As Double eri = Exp(r*(i-1)) |
|||
Dim As Double guess = (n0*eri) / (1+n0*(eri-1) / K) |
|||
Dim As Double diff = guess-actual(i) |
|||
sq += diff*diff |
|||
Next |
|||
Return sq |
|||
End Function |
|||
Function solve(f As Function (As Double) As Double, guess As Double = 0.5, epsilon As Double = 0.0) As Double |
|||
Dim As Double f0 = f(guess) |
|||
Dim As Double delta = Iif(guess, guess, 1) |
|||
Dim As Double factor = 2 ' double until f0 best then halve until delta <= epsilon |
|||
While delta > epsilon And guess <> guess-delta |
|||
Dim As Double nf = f(guess-delta) |
|||
If nf < f0 Then |
|||
f0 = nf |
|||
guess -= delta |
|||
Else |
|||
nf = f(guess+delta) |
|||
If nf < f0 Then |
|||
f0 = nf |
|||
guess += delta |
|||
Else |
|||
factor = 0.5 |
|||
End If |
|||
End If |
|||
delta *= factor |
|||
Wend |
|||
Return guess |
|||
End Function |
|||
Dim As Double r = solve(@f1) |
|||
Dim As Double R0 = Exp(12 * r) |
|||
Print Using "r = #.##########, R0 = #.########"; r; R0 |
|||
Sleep</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.1123021757, R0 = 3.84827928</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{libheader|lm}} |
{{libheader|lm}} |
||
This uses the Levenberg-Marquardt method. |
This uses the Levenberg-Marquardt method. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 198: | Line 611: | ||
fmt.Printf("The logistic curve r for the world data is %.8f\n", r) |
fmt.Printf("The logistic curve r for the world data is %.8f\n", r) |
||
fmt.Printf("R0 is then approximately equal to %.7f\n", math.Exp(12*r)) |
fmt.Printf("R0 is then approximately equal to %.7f\n", math.Exp(12*r)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 205: | Line 618: | ||
R0 is then approximately equal to 3.8482793 |
R0 is then approximately equal to 3.8482793 |
||
</pre> |
</pre> |
||
=={{header|Java}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="java">import java.util.List; |
|||
import java.util.function.Function; |
|||
public class LogisticCurveFitting { |
|||
private static final double K = 7.8e9; |
|||
private static final int N0 = 27; |
|||
private static final List<Double> ACTUAL = List.of( |
|||
27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0, 60.0, |
|||
61.0, 61.0, 66.0, 83.0, 219.0, 239.0, 392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0, |
|||
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0, 20615.0, 24522.0, 28273.0, |
|||
31491.0, 34933.0, 37552.0, 40540.0, 43105.0, 45177.0, 60328.0, 64543.0, 67103.0, |
|||
69265.0, 71332.0, 73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0, |
|||
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0, 90664.0, 93077.0, |
|||
95316.0, 98172.0, 102133.0, 105824.0, 109695.0, 114232.0, 118610.0, 125497.0, |
|||
133852.0, 143227.0, 151367.0, 167418.0, 180096.0, 194836.0, 213150.0, 242364.0, |
|||
271106.0, 305117.0, 338133.0, 377918.0, 416845.0, 468049.0, 527767.0, 591704.0, |
|||
656866.0, 715353.0, 777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0 |
|||
); |
|||
private static double f(double r) { |
|||
var sq = 0.0; |
|||
var len = ACTUAL.size(); |
|||
for (int i = 0; i < len; i++) { |
|||
var eri = Math.exp(r * i); |
|||
var guess = (N0 * eri) / (1.0 + N0 * (eri - 1.0) / K); |
|||
var diff = guess - ACTUAL.get(i); |
|||
sq += diff * diff; |
|||
} |
|||
return sq; |
|||
} |
|||
private static double solve(Function<Double, Double> fn) { |
|||
return solve(fn, 0.5, 0.0); |
|||
} |
|||
private static double solve(Function<Double, Double> fn, double guess, double epsilon) { |
|||
double delta; |
|||
if (guess != 0.0) { |
|||
delta = guess; |
|||
} else { |
|||
delta = 1.0; |
|||
} |
|||
var f0 = fn.apply(guess); |
|||
var factor = 2.0; |
|||
while (delta > epsilon && guess != guess - delta) { |
|||
var nf = fn.apply(guess - delta); |
|||
if (nf < f0) { |
|||
f0 = nf; |
|||
guess -= delta; |
|||
} else { |
|||
nf = fn.apply(guess + delta); |
|||
if (nf < f0) { |
|||
f0 = nf; |
|||
guess += delta; |
|||
} else { |
|||
factor = 0.5; |
|||
} |
|||
} |
|||
delta *= factor; |
|||
} |
|||
return guess; |
|||
} |
|||
public static void main(String[] args) { |
|||
var r = solve(LogisticCurveFitting::f); |
|||
var r0 = Math.exp(12.0 * r); |
|||
System.out.printf("r = %.16f, R0 = %.16f\n", r, r0); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.1123021756975501, R0 = 3.8482792809167194</pre> |
|||
=={{header|jq}}== |
|||
{{trans|Wren}} |
|||
{{works with|jq}} |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
<syntaxhighlight lang="jq">def K: 7800000000; # approx world population |
|||
def n0: 27; # number of cases at day 0 |
|||
def y: [ |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, |
|||
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, |
|||
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177, |
|||
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723, |
|||
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, |
|||
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133, |
|||
105824, 109695, 114232, 118610, 125497, 133852, 143227, |
|||
151367, 167418, 180096, 194836, 213150, 242364, 271106, |
|||
305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, |
|||
1174652 |
|||
]; |
|||
def f: |
|||
. as $r |
|||
| reduce range(0; y|length) as $i (0; |
|||
(($r*$i)|exp) as $eri |
|||
| ((n0*$eri)/(1+n0*($eri-1)/K) - y[$i]) as $dst |
|||
| . + $dst * $dst); |
|||
def solve(f; $guess; $epsilon): |
|||
{ f0: ($guess|f), |
|||
delta: $guess, |
|||
factor: 2, |
|||
$guess} |
|||
# double until f0 best, then halve until delta <= epsilon |
|||
| until(.delta <= $epsilon; |
|||
.nf = ((.guess - .delta)|f) |
|||
| if .nf < .f0 |
|||
then .f0 = .nf |
|||
| .guess = .guess - .delta |
|||
else .nf = ((.guess + .delta)|f) |
|||
| if .nf < .f0 |
|||
then .f0 = .nf |
|||
| .guess = .guess + .delta |
|||
else .factor = 0.5 |
|||
end |
|||
end |
|||
| .delta = .delta * .factor ) |
|||
| .guess; |
|||
((solve(f; 0.5; 0) * 1e10)|round / 1e10 ) as $r |
|||
| ((((12 * $r)|exp) * 1e8)|round / 1e8) as $R0 |
|||
| "r = \($r), R0 = \($R0)"</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
r = 0.1123021757, R0 = 3.84827928 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
< |
<syntaxhighlight lang="julia">using LsqFit |
||
const K = 7_800_000_000 # approximate world population |
const K = 7_800_000_000 # approximate world population |
||
Line 243: | Line 796: | ||
confidence_interval(fit, 0.05)) |
confidence_interval(fit, 0.05)) |
||
println("Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ ", exp(12r[1])) |
println("Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ ", exp(12r[1])) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
The logistic curve r for the world data is: [0.11230217572265622] |
The logistic curve r for the world data is: [0.11230217572265622] |
||
Line 249: | Line 802: | ||
Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ 3.8482792820761063 |
Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ 3.8482792820761063 |
||
</pre> |
</pre> |
||
=={{header|Kotlin}}== |
|||
{{trans|D}} |
|||
<syntaxhighlight lang="scala">import kotlin.math.exp |
|||
const val K = 7.8e9 |
|||
const val n0 = 27 |
|||
val actual = arrayOf( |
|||
27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0, 60.0, |
|||
61.0, 61.0, 66.0, 83.0, 219.0, 239.0, 392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0, |
|||
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0, 20615.0, 24522.0, 28273.0, |
|||
31491.0, 34933.0, 37552.0, 40540.0, 43105.0, 45177.0, 60328.0, 64543.0, 67103.0, |
|||
69265.0, 71332.0, 73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0, |
|||
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0, 90664.0, 93077.0, |
|||
95316.0, 98172.0, 102133.0, 105824.0, 109695.0, 114232.0, 118610.0, 125497.0, |
|||
133852.0, 143227.0, 151367.0, 167418.0, 180096.0, 194836.0, 213150.0, 242364.0, |
|||
271106.0, 305117.0, 338133.0, 377918.0, 416845.0, 468049.0, 527767.0, 591704.0, |
|||
656866.0, 715353.0, 777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0 |
|||
) |
|||
fun f(r: Double): Double { |
|||
var sq = 0.0 |
|||
val len = actual.size |
|||
for (i in 0 until len) { |
|||
val eri = exp(r * i) |
|||
val guess = (n0 * eri) / (1.0 + n0 * (eri - 1.0) / K) |
|||
val diff = guess - actual[i] |
|||
sq += diff * diff |
|||
} |
|||
return sq |
|||
} |
|||
fun solve(fn: (Double) -> Double, guess: Double = 0.5, epsilon: Double = 0.0): Double { |
|||
var guess2 = guess |
|||
var delta = if (guess2 != 0.0) guess2 else 1.0 |
|||
var f0 = fn(guess2) |
|||
var factor = 2.0 |
|||
while (delta > epsilon && guess2 != guess2 - delta) { |
|||
var nf = fn(guess2 - delta) |
|||
if (nf < f0) { |
|||
f0 = nf |
|||
guess2 -= delta |
|||
} else { |
|||
nf = fn(guess2 + delta) |
|||
if (nf < f0) { |
|||
f0 = nf |
|||
guess2 += delta |
|||
} else { |
|||
factor = 0.5 |
|||
} |
|||
} |
|||
delta *= factor |
|||
} |
|||
return guess2 |
|||
} |
|||
fun main() { |
|||
val r = solve(::f) |
|||
val r0 = exp(12.0 * r) |
|||
println("r = $r, R0 = $r0") |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.11230217569755005, R0 = 3.8482792809167194</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="nim">import math, sugar |
|||
const |
|||
K = 7.8e9 |
|||
N0 = 27 |
|||
Actual = [ 27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0, |
|||
59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0, |
|||
60.0, 61.0, 61.0, 66.0, 83.0, 219.0, 239.0, |
|||
392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0, |
|||
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0, |
|||
20615.0, 24522.0, 28273.0, 31491.0, 34933.0, 37552.0, 40540.0, |
|||
43105.0, 45177.0, 60328.0, 64543.0, 67103.0, 69265.0, 71332.0, |
|||
73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0, |
|||
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0, |
|||
90664.0, 93077.0, 95316.0, 98172.0, 102133.0, 105824.0, 109695.0, |
|||
114232.0, 118610.0, 125497.0, 133852.0, 143227.0, 151367.0, 167418.0, |
|||
180096.0, 194836.0, 213150.0, 242364.0, 271106.0, 305117.0, 338133.0, |
|||
377918.0, 416845.0, 468049.0, 527767.0, 591704.0, 656866.0, 715353.0, |
|||
777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0] |
|||
type Func = (float) -> float |
|||
func f(r: float): float = |
|||
for i in 0..Actual.high: |
|||
let |
|||
eri = exp(r * i.toFloat) |
|||
guess = (N0 * eri) / (1 + N0 * (eri - 1) / K) |
|||
diff = guess - Actual[i] |
|||
result += diff * diff |
|||
func solve(fn: Func; guess = 0.5, epsilon = 0.0): float = |
|||
result = guess |
|||
var delta = if result != 0: result else: 1.0 |
|||
var f0 = fn(result) |
|||
var factor = 2.0 |
|||
while delta > epsilon and result != result - delta: |
|||
var nf = fn(result - delta) |
|||
if nf < f0: |
|||
f0 = nf |
|||
result -= delta |
|||
else: |
|||
nf = fn(result + delta) |
|||
if nf < f0: |
|||
f0 = nf |
|||
result += delta |
|||
else: |
|||
factor = 0.5 |
|||
delta *= factor |
|||
let r = f.solve() |
|||
let r0 = exp(12 * r) |
|||
echo "r = ", r, ", R0 = ", r0</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.11230217569755, R0 = 3.848279280916719</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{trans|Phix}} |
{{trans|Phix}} |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
Line 300: | Line 982: | ||
my $r = solve(\&logistic_func, 0.5, 0); |
my $r = solve(\&logistic_func, 0.5, 0); |
||
my $R0 = exp(12 * $r); |
my $R0 = exp(12 * $r); |
||
printf "r = %%(%.3f), R0 = %%(%.3f)\n", $r, $R0;</ |
printf "r = %%(%.3f), R0 = %%(%.3f)\n", $r, $R0;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>r = %(0.112), R0 = %(3.848)</pre> |
<pre>r = %(0.112), R0 = %(3.848)</pre> |
||
Line 306: | Line 988: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Simplified, my interpretation of shift-cutting (from [[wp:Non-linear_least_squares]]) |
Simplified, my interpretation of shift-cutting (from [[wp:Non-linear_least_squares]]) |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>-- 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 |
|||
-- halve until delta<=epsilon |
|||
-- 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: #008080;">else</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 |
|||
<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> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 370: | Line 1,055: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
Uses NumPy/SciPy's optimize package. |
Uses NumPy/SciPy's optimize package. |
||
< |
<syntaxhighlight lang="python">import numpy as np |
||
import scipy.optimize as opt |
import scipy.optimize as opt |
||
Line 400: | Line 1,085: | ||
", with covariance of", cov) |
", with covariance of", cov) |
||
print("The calculated R0 is then", np.exp(12 * r)) |
print("The calculated R0 is then", np.exp(12 * r)) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
The r for the world Covid-19 data is: [0.11230218] , with covariance of [[2.46164331e-08]] |
The r for the world Covid-19 data is: [0.11230218] , with covariance of [[2.46164331e-08]] |
||
Line 407: | Line 1,092: | ||
=={{header|R}}== |
=={{header|R}}== |
||
<syntaxhighlight lang="r"> |
|||
<lang R> |
|||
library(minpack.lm) |
library(minpack.lm) |
||
Line 442: | Line 1,127: | ||
cat("Therefore, R0 is approximately: ", exp(12 * coef(nls.out))) |
cat("Therefore, R0 is approximately: ", exp(12 * coef(nls.out))) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
The r for the model is: 0.1123022 |
The r for the model is: 0.1123022 |
||
Line 454: | Line 1,139: | ||
Original task numbers of cases per day as of April 05 plus updated, as of today, May 11 numbers. |
Original task numbers of cases per day as of April 05 plus updated, as of today, May 11 numbers. |
||
<lang |
<syntaxhighlight lang="raku" line>my $K = 7_800_000_000; # population |
||
my $n0 = 27; # cases @ day 0 |
my $n0 = 27; # cases @ day 0 |
||
Line 530: | Line 1,215: | ||
say "Reproductive rate: R0 = ", $R0.fmt('%08f'); |
say "Reproductive rate: R0 = ", $R0.fmt('%08f'); |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre>For period: Dec 31 - Apr 5 |
<pre>For period: Dec 31 - Apr 5 |
||
Line 539: | Line 1,224: | ||
Instantaneous rate of growth: r = 0.093328 |
Instantaneous rate of growth: r = 0.093328 |
||
Reproductive rate: R0 = 3.064672</pre> |
Reproductive rate: R0 = 3.064672</pre> |
||
=={{header|Ruby}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="ruby">K = 7.9e9 |
|||
N0 = 27 |
|||
ACTUAL = [ |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, |
|||
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, |
|||
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, |
|||
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, |
|||
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, |
|||
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, |
|||
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, |
|||
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 |
|||
] |
|||
def f(r) |
|||
sq = 0.0 |
|||
len = ACTUAL.length |
|||
for i in 1 .. len |
|||
j = i - 1 |
|||
eri = Math.exp(r * j) |
|||
guess = (N0 * eri) / (1 + N0 * (eri - 1.0) / K) |
|||
diff = guess - ACTUAL[j] |
|||
sq += diff * diff |
|||
end |
|||
return sq |
|||
end |
|||
def solve(fn, guess=0.5, epsilon=0.0) |
|||
delta = guess ? guess : 1.0 |
|||
f0 = send(fn, guess) |
|||
factor = 2.0 |
|||
while delta > epsilon and guess != guess - delta |
|||
nf = send(fn, guess - delta) |
|||
if nf < f0 then |
|||
f0 = nf |
|||
guess -= delta |
|||
else |
|||
nf = send(fn, guess + delta) |
|||
if nf < f0 then |
|||
f0 = nf |
|||
guess += delta |
|||
else |
|||
factor = 0.5 |
|||
end |
|||
end |
|||
delta *= factor |
|||
end |
|||
return guess |
|||
end |
|||
def main |
|||
r = solve(:f) |
|||
r0 = Math.exp(12.0 * r) |
|||
print "r = ", r, ", R0 = ", r0, "\n" |
|||
end |
|||
main()</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.11230215850810021, R0 = 3.8482784871191575</pre> |
|||
=={{header|V (Vlang)}}== |
|||
{{trans|wren}} |
|||
<syntaxhighlight lang="v (vlang)">import math |
|||
const ( |
|||
k = 7800000000 // approx world population |
|||
n0 = 27 // number of cases at day 0 |
|||
y = [ |
|||
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, |
|||
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, |
|||
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, |
|||
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177, |
|||
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723, |
|||
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, |
|||
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133, |
|||
105824, 109695, 114232, 118610, 125497, 133852, 143227, |
|||
151367, 167418, 180096, 194836, 213150, 242364, 271106, |
|||
305117, 338133, 377918, 416845, 468049, 527767, 591704, |
|||
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, |
|||
1174652 |
|||
] |
|||
) |
|||
fn f(r f64) f64 { |
|||
mut sq := 0.0 |
|||
for i in 0..y.len { |
|||
eri := math.exp(r * i) |
|||
dst := (n0 * eri) / (1 + n0 * (eri - 1) / k) -y[i] |
|||
sq = sq + dst * dst |
|||
} |
|||
return sq |
|||
} |
|||
fn solve(fu fn(f64)f64, g f64, epsilon int) f64 { |
|||
mut guess := g |
|||
mut f0 := fu(guess) |
|||
mut delta := guess |
|||
mut factor := 2.0 |
|||
for delta > epsilon { |
|||
mut nf := fu(guess - delta) |
|||
if nf < f0 { |
|||
f0 = nf |
|||
guess -= delta |
|||
} else { |
|||
nf = fu(guess + delta) |
|||
if nf < f0 { |
|||
f0 = nf |
|||
guess += delta |
|||
} else { |
|||
factor = 0.5 |
|||
} |
|||
} |
|||
delta *= factor |
|||
} |
|||
return guess |
|||
} |
|||
fn main() { |
|||
r := math.round(solve(f, 0.5, 0) * 1e10) / 1e10 |
|||
r0 := math.round(math.exp(12 * r) * 1e8) / 1e8 |
|||
println("r = $r, R0 = $r0") |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>r = 0.1123021757, R0 = 3.84827928</pre> |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
{{trans|Phix}} |
{{trans|Phix}} |
||
<syntaxhighlight lang="wren">var K = 7800000000 // approx world population |
|||
{{libheader|Wren-math}} |
|||
<lang ecmascript>import "/math" for Math |
|||
var K = 7800000000 // approx world population |
|||
var n0 = 27 // number of cases at day 0 |
var n0 = 27 // number of cases at day 0 |
||
Line 566: | Line 1,374: | ||
var sq = 0 |
var sq = 0 |
||
for (i in 0...y.count) { |
for (i in 0...y.count) { |
||
var eri = |
var eri = (r*i).exp |
||
var dst = (n0*eri)/(1+n0*(eri-1)/K) - y[i] |
var dst = (n0*eri)/(1+n0*(eri-1)/K) - y[i] |
||
sq = sq + dst * dst |
sq = sq + dst * dst |
||
Line 597: | Line 1,405: | ||
var r = (solve.call(f, 0.5, 0) * 1e10).round / 1e10 |
var r = (solve.call(f, 0.5, 0) * 1e10).round / 1e10 |
||
var R0 = ( |
var R0 = ((12 * r).exp * 1e8).round / 1e8 |
||
System.print("r = %(r), R0 = %(R0)")</ |
System.print("r = %(r), R0 = %(R0)")</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
r = 0.1123021757, R0 = 3.84827928 |
r = 0.1123021757, R0 = 3.84827928 |
||
</pre> |
|||
=={{header|XPL0}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang "XPL0">def K = 7.8e9; |
|||
def N0 = 27.; |
|||
real Actual; |
|||
int ActualSize; |
|||
func real Func(R); |
|||
real R; |
|||
real Sq, Eri, Guess, Diff; |
|||
int I; |
|||
[Sq:= 0.; |
|||
for I:= 0 to ActualSize-1 do |
|||
[Eri:= Exp(R * float(I)); |
|||
Guess:= N0 * Eri / (1. + N0 * (Eri-1.) / K); |
|||
Diff:= Guess - Actual(I); |
|||
Sq:= Sq + Diff*Diff; |
|||
]; |
|||
return Sq; |
|||
]; |
|||
func real Solve(Guess, Epsilon); |
|||
real Guess, Epsilon; |
|||
real Delta, F0, Factor, NF; |
|||
[Delta:= if Guess # 0. then Guess else 1.; |
|||
F0:= Func(Guess); |
|||
Factor:= 2.; |
|||
while Delta > Epsilon and Guess # Guess-Delta do |
|||
[NF:= Func(Guess - Delta); |
|||
if NF < F0 then |
|||
[F0:= NF; |
|||
Guess:= Guess - Delta; |
|||
] |
|||
else |
|||
[NF:= Func(Guess + Delta); |
|||
if NF < F0 then |
|||
[F0:= NF; |
|||
Guess:= Guess + Delta; |
|||
] |
|||
else Factor:= 0.5; |
|||
]; |
|||
Delta:= Delta * Factor; |
|||
]; |
|||
return Guess; |
|||
]; |
|||
real R, R0; |
|||
[Actual:= [ |
|||
27., 27., 27., 44., 44., 59., 59., 59., 59., 59., 59., 59., 59., 60., 60., |
|||
61., 61., 66., 83., 219., 239., 392., 534., 631., 897., 1350., 2023., 2820., |
|||
4587., 6067., 7823., 9826., 11946., 14554., 17372., 20615., 24522., 28273., |
|||
31491., 34933., 37552., 40540., 43105., 45177., 60328., 64543., 67103., |
|||
69265., 71332., 73327., 75191., 75723., 76719., 77804., 78812., 79339., |
|||
80132., 80995., 82101., 83365., 85203., 87024., 89068., 90664., 93077., |
|||
95316., 98172., 102133., 105824., 109695., 114232., 118610., 125497., |
|||
133852., 143227., 151367., 167418., 180096., 194836., 213150., 242364., |
|||
271106., 305117., 338133., 377918., 416845., 468049., 527767., 591704., |
|||
656866., 715353., 777796., 851308., 928436., 1000249., 1082054., 1174652.]; |
|||
ActualSize:= 97; |
|||
R:= Solve(0.5, 0.0); |
|||
R0:= Exp(12.*R); |
|||
Text(0, "R = "); RlOut(0, R); CrLf(0); |
|||
Text(0, "R0 = "); RlOut(0, R0); CrLf(0); |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
R = 0.11230 |
|||
R0 = 3.84828 |
|||
</pre> |
</pre> |
Latest revision as of 15:51, 17 December 2023
You are encouraged to solve this task according to the task description, using any language you may know.
The least-squares method (see references below) in statistics is used to fit data to the best of a family of similar curves by finding the parameters for a curve which minimizes the total of the distances from each data point to the curve.
Often, the curve used is a straight line, in which case the method is also called linear regression. If a curve which uses logarithmic growth is fit, the method can be called logistic regression.
A commonly used family of functions used in statistical studies of populations, including the growth of epidemics, are curves akin to the logistic curve:
ƒ(x) = L / (1 + e-k(x-x0))
Though predictions based on fitting to such curves may error, especially if used to extrapolate from incomplete data, curves similar to the logistic curve have had good fits in population studies, including modeling the growth of past epidemics.
Given the following daily world totals since December 31, 2019 for persons who have become infected with the novel coronavirus Covid-19:
Daily totals: 27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, 61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, 133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, 271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, 656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652
Use the following variant of the logistic curve as a formula:
ƒ(t) = n0 e(r t) / ((1 + n0 (e(r t) - 1) / K)
Where:
- r is the rate of growth of the infection in the population.
- K is the world population, about 7.8 billion.
- n0 is 27, the number of cases found in China at the start of the pandemic.
The R0 of an infection (different from r above) is a measure of how many new individuals will become infected for every individual currently infected. It is an important measure of how quickly an infectious disease may spread.
R0 is related to the logistic curve's r parameter by the formula:
r ≈ ln(R0) / G
where G the generation time, is roughly the sum of the incubation time, perhaps 5 days, and the mean contagion period, perhaps 7 days, so, for covid-19, roughly we have:
R0 ≈ e12r
- Task
- Demonstrate code that finds a least-squares fits of the curve to the data.
- Show the calculated r for the logistic curve.
- Show the final R0 parameter you calculate from the logistic curve r value parameter.
- See also
11l
-V K = 7.8e9
-V n0 = 27
-V actual = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820,
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273,
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103,
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339,
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497,
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364,
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652
]
F f(Float r) -> Float
V sq = 0.0
L(act) :actual
V eri = exp(r * L.index)
V guess = (:n0 * eri) / (1 + :n0 * (eri - 1) / :K)
V diff = guess - act
sq += diff * diff
R sq
F solve(func, =guess = 0.5, epsilon = 0.0)
V delta = I guess != 0 {guess} E 1
V f0 = func(guess)
V factor = 2.0
L delta > epsilon & guess != guess - delta
V nf = func(guess - delta)
I nf < f0
f0 = nf
guess -= delta
E
nf = func(guess + delta)
I nf < f0
f0 = nf
guess += delta
E
factor = 0.5
delta *= factor
R guess
V r = solve(f)
V R0 = exp(12 * r)
print(‘r = #.6, R0 = #.6’.format(r, R0))
- Output:
r = 0.112302, R0 = 3.848279
Ada
with Ada.Text_Io;
with Ada.Numerics.Generic_Elementary_Functions;
procedure Curve_Fitting is
type Real is new Long_Float;
type Fuction_Access is access function (X : Real) return Real;
type Time_Type is new Natural; -- Days
package Real_Io is new Ada.Text_Io.Float_Io (Real);
package Math is new Ada.Numerics.Generic_Elementary_Functions (Real);
Actual : constant array (Time_Type range <>) of Integer :=
(27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820,
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273,
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103,
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339,
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497,
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364,
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652);
N0 : constant Real := Real (Actual (Actual'First)); -- Initially infected
K : constant Real := 7.8e9; -- World population apx.
function Logistic_Curve_Function (R : Real) return Real is
sq : Real := 0.0;
begin
for I in Actual'range loop
declare
Eri : constant Real := Math.Exp (R * Real (I));
Guess : constant Real := (N0 * Eri) / (1.0 + N0 * (Eri - 1.0) / K);
Diff : constant Real := Guess - Real (Actual (I));
begin
Sq := Sq + Diff * Diff;
end;
end loop;
return sq;
end Logistic_Curve_Function;
procedure Solve (Func : Fuction_Access;
Guess : Real := 0.5;
Epsilon : Real := 0.0;
New_Guess : out Real)
is
Delt : Real := (if guess /= 0.0 then guess else 1.0);
f0 : Real := Func (Guess);
factor : Real := 2.0;
begin
New_Guess := Guess;
while
Delt > Epsilon and New_Guess /= New_Guess - Delt
loop
declare
nf : real := func (New_guess - delt);
begin
if nf < f0 then
f0 := nf;
New_guess := New_guess - delt;
else
nf := func (New_guess + delt);
if nf < f0 then
f0 := nf;
New_guess := New_guess + delt;
else
Factor := 0.5;
end if;
end if;
end;
Delt := Delt * factor;
end loop;
end Solve;
use Ada.Text_Io;
Incubation_Time : constant Real := 5.0; -- Days
Contagion_Period : constant Real := 7.0; -- Days
Generation_Time : constant Real := Incubation_Time + Contagion_Period;
R : Real; -- Rate
R0 : Real; -- Infection rate
begin
Solve (Logistic_Curve_Function'Access, New_Guess => R);
R0 := Math.Exp (Generation_Time * R);
Put ("r = "); Real_IO.Put (R, Exp => 0, Aft => 7);
Put (", R0 = "); Real_Io.Put (R0, Exp => 0, Aft => 7);
New_Line;
end Curve_Fitting;
- Output:
r = 0.1123022, R0 = 3.8482793
Arturo
K: 7.8e9
N0: 27
Actual: [
27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0,
59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0,
60.0, 61.0, 61.0, 66.0, 83.0, 219.0, 239.0,
392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0,
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0,
20615.0, 24522.0, 28273.0, 31491.0, 34933.0, 37552.0, 40540.0,
43105.0, 45177.0, 60328.0, 64543.0, 67103.0, 69265.0, 71332.0,
73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0,
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0,
90664.0, 93077.0, 95316.0, 98172.0, 102133.0, 105824.0, 109695.0,
114232.0, 118610.0, 125497.0, 133852.0, 143227.0, 151367.0, 167418.0,
180096.0, 194836.0, 213150.0, 242364.0, 271106.0, 305117.0, 338133.0,
377918.0, 416845.0, 468049.0, 527767.0, 591704.0, 656866.0, 715353.0,
777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0
]
f: function [r][
result: 0
loop 0..dec size Actual 'i [
eri: exp r * to :floating i
guess: (N0 * eri) // (1 + N0 * (eri - 1) // K)
diff: guess - Actual\[i]
result: result + diff * diff
]
return result
]
solve: function [fn][
guess: 0.5
eps: 0.0
result: guess
delta: guess
f0: call fn @[result]
factor: 2.0
while [and? -> delta > eps
-> result <> result - delta][
nf: call fn @[result-delta]
if? nf < f0 [
f0: nf
result: result - delta
]
else [
nf: call fn @[result+delta]
if? nf < f0 [
f0: nf
result: result + delta
]
else [
factor: 0.5
]
]
delta: delta * factor
]
return result
]
r: solve 'f
r0: exp 12 * r
print ["r =" r]
print ["R0 =" r0]
- Output:
r = 0.11230217569755 R0 = 3.848279280916719
C
#include <math.h>
#include <stdio.h>
const double K = 7.8e9;
const int n0 = 27;
const double actual[] = {
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820,
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273,
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103,
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339,
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497,
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364,
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652
};
const size_t actual_size = sizeof(actual) / sizeof(double);
double f(double r) {
double sq = 0;
size_t i;
for (i = 0; i < actual_size; ++i) {
double eri = exp(r * i);
double guess = (n0 * eri) / (1 + n0 * (eri - 1) / K);
double diff = guess - actual[i];
sq += diff * diff;
}
return sq;
}
double solve(double (*fn)(double), double guess, double epsilon) {
double delta, f0, factor;
for (delta = guess ? guess : 1, f0 = fn(guess), factor = 2;
delta > epsilon && guess != guess - delta;
delta *= factor) {
double nf = (*fn)(guess - delta);
if (nf < f0) {
f0 = nf;
guess -= delta;
} else {
nf = fn(guess + delta);
if (nf < f0) {
f0 = nf;
guess += delta;
} else {
factor = 0.5;
}
}
}
return guess;
}
double solve_default(double (*fn)(double)) {
return solve(fn, 0.5, 0);
}
int main() {
double r = solve_default(f);
double R0 = exp(12 * r);
printf("r = %f, R0 = %f\n", r, R0);
return 0;
}
- Output:
r = 0.112302, R0 = 3.848279
C++
#include <cmath>
#include <functional>
#include <iostream>
constexpr double K = 7.8e9;
constexpr int n0 = 27;
constexpr double actual[] = {
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820,
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273,
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103,
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339,
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497,
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364,
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652
};
double f(double r) {
double sq = 0;
constexpr size_t len = std::size(actual);
for (size_t i = 0; i < len; ++i) {
double eri = std::exp(r * i);
double guess = (n0 * eri)/(1 + n0 * (eri - 1)/K);
double diff = guess - actual[i];
sq += diff * diff;
}
return sq;
}
double solve(std::function<double(double)> fn, double guess=0.5, double epsilon=0) {
for (double delta = guess ? guess : 1, f0 = fn(guess), factor = 2;
delta > epsilon && guess != guess - delta;
delta *= factor) {
double nf = fn(guess - delta);
if (nf < f0) {
f0 = nf;
guess -= delta;
} else {
nf = fn(guess + delta);
if (nf < f0) {
f0 = nf;
guess += delta;
} else
factor = 0.5;
}
}
return guess;
}
int main() {
double r = solve(f);
double R0 = std::exp(12 * r);
std::cout << "r = " << r << ", R0 = " << R0 << '\n';
return 0;
}
- Output:
r = 0.112302, R0 = 3.84828
D
import std.math;
import std.stdio;
immutable K = 7.8e9;
immutable n0 = 27;
immutable actual = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820,
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273,
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103,
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339,
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497,
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364,
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652
];
double f(double r) {
double sq = 0.0;
auto len = actual.length;
for (size_t i = 0; i < len; ++i) {
auto eri = exp(r * i);
auto guess = (n0 * eri) / (1 + n0 * (eri - 1) / K);
auto diff = guess - actual[i];
sq += diff * diff;
}
return sq;
}
double solve(double function(double) fn, double guess = 0.5, double epsilon = 0) {
for (double delta = guess ? guess : 1, f0 = fn(guess), factor = 2;
delta > epsilon && guess != guess - delta;
delta *= factor
) {
double nf = fn(guess - delta);
if (nf < f0) {
f0 = nf;
guess -= delta;
} else {
nf = fn(guess + delta);
if (nf < f0) {
f0 = nf;
guess += delta;
} else {
factor = 0.5;
}
}
}
return guess;
}
void main() {
double r = solve(&f);
double R0 = exp(12 * r);
writeln("r = ", r, ", R0 = ", R0);
}
- Output:
r = 0.112302, R0 = 3.84828
FreeBASIC
Const K = 7.8e9 ' approx world population
Const n0 = 27 ' number of cases at day 0
Dim Shared As Integer actual(1 To ...) => { _
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, _
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, _
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, _
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177, _
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723, _
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, _
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133, _
105824, 109695, 114232, 118610, 125497, 133852, 143227, _
151367, 167418, 180096, 194836, 213150, 242364, 271106, _
305117, 338133, 377918, 416845, 468049, 527767, 591704, _
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 }
Function f1(r As Double) As Double
Dim As Double sq = 0
For i As Integer = 1 To Ubound(actual)
Dim As Double eri = Exp(r*(i-1))
Dim As Double guess = (n0*eri) / (1+n0*(eri-1) / K)
Dim As Double diff = guess-actual(i)
sq += diff*diff
Next
Return sq
End Function
Function solve(f As Function (As Double) As Double, guess As Double = 0.5, epsilon As Double = 0.0) As Double
Dim As Double f0 = f(guess)
Dim As Double delta = Iif(guess, guess, 1)
Dim As Double factor = 2 ' double until f0 best then halve until delta <= epsilon
While delta > epsilon And guess <> guess-delta
Dim As Double nf = f(guess-delta)
If nf < f0 Then
f0 = nf
guess -= delta
Else
nf = f(guess+delta)
If nf < f0 Then
f0 = nf
guess += delta
Else
factor = 0.5
End If
End If
delta *= factor
Wend
Return guess
End Function
Dim As Double r = solve(@f1)
Dim As Double R0 = Exp(12 * r)
Print Using "r = #.##########, R0 = #.########"; r; R0
Sleep
- Output:
r = 0.1123021757, R0 = 3.84827928
Go
This uses the Levenberg-Marquardt method.
package main
import (
"fmt"
"github.com/maorshutman/lm"
"log"
"math"
)
const (
K = 7_800_000_000 // approx world population
n0 = 27 // number of cases at day 0
)
var y = []float64{
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652,
}
func f(dst, p []float64) {
for i := 0; i < len(y); i++ {
t := float64(i)
dst[i] = (n0*math.Exp(p[0]*t))/(1+n0*(math.Exp(p[0]*t)-1)/K) - y[i]
}
}
func main() {
j := lm.NumJac{Func: f}
prob := lm.LMProblem{
Dim: 1,
Size: len(y),
Func: f,
Jac: j.Jac,
InitParams: []float64{0.5},
Tau: 1e-6,
Eps1: 1e-8,
Eps2: 1e-8,
}
res, err := lm.LM(prob, &lm.Settings{Iterations: 100, ObjectiveTol: 1e-16})
if err != nil {
log.Fatal(err)
}
r := res.X[0]
fmt.Printf("The logistic curve r for the world data is %.8f\n", r)
fmt.Printf("R0 is then approximately equal to %.7f\n", math.Exp(12*r))
}
- Output:
The logistic curve r for the world data is 0.11230218 R0 is then approximately equal to 3.8482793
Java
import java.util.List;
import java.util.function.Function;
public class LogisticCurveFitting {
private static final double K = 7.8e9;
private static final int N0 = 27;
private static final List<Double> ACTUAL = List.of(
27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0, 60.0,
61.0, 61.0, 66.0, 83.0, 219.0, 239.0, 392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0,
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0, 20615.0, 24522.0, 28273.0,
31491.0, 34933.0, 37552.0, 40540.0, 43105.0, 45177.0, 60328.0, 64543.0, 67103.0,
69265.0, 71332.0, 73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0,
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0, 90664.0, 93077.0,
95316.0, 98172.0, 102133.0, 105824.0, 109695.0, 114232.0, 118610.0, 125497.0,
133852.0, 143227.0, 151367.0, 167418.0, 180096.0, 194836.0, 213150.0, 242364.0,
271106.0, 305117.0, 338133.0, 377918.0, 416845.0, 468049.0, 527767.0, 591704.0,
656866.0, 715353.0, 777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0
);
private static double f(double r) {
var sq = 0.0;
var len = ACTUAL.size();
for (int i = 0; i < len; i++) {
var eri = Math.exp(r * i);
var guess = (N0 * eri) / (1.0 + N0 * (eri - 1.0) / K);
var diff = guess - ACTUAL.get(i);
sq += diff * diff;
}
return sq;
}
private static double solve(Function<Double, Double> fn) {
return solve(fn, 0.5, 0.0);
}
private static double solve(Function<Double, Double> fn, double guess, double epsilon) {
double delta;
if (guess != 0.0) {
delta = guess;
} else {
delta = 1.0;
}
var f0 = fn.apply(guess);
var factor = 2.0;
while (delta > epsilon && guess != guess - delta) {
var nf = fn.apply(guess - delta);
if (nf < f0) {
f0 = nf;
guess -= delta;
} else {
nf = fn.apply(guess + delta);
if (nf < f0) {
f0 = nf;
guess += delta;
} else {
factor = 0.5;
}
}
delta *= factor;
}
return guess;
}
public static void main(String[] args) {
var r = solve(LogisticCurveFitting::f);
var r0 = Math.exp(12.0 * r);
System.out.printf("r = %.16f, R0 = %.16f\n", r, r0);
}
}
- Output:
r = 0.1123021756975501, R0 = 3.8482792809167194
jq
Works with gojq, the Go implementation of jq
def K: 7800000000; # approx world population
def n0: 27; # number of cases at day 0
def y: [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652
];
def f:
. as $r
| reduce range(0; y|length) as $i (0;
(($r*$i)|exp) as $eri
| ((n0*$eri)/(1+n0*($eri-1)/K) - y[$i]) as $dst
| . + $dst * $dst);
def solve(f; $guess; $epsilon):
{ f0: ($guess|f),
delta: $guess,
factor: 2,
$guess}
# double until f0 best, then halve until delta <= epsilon
| until(.delta <= $epsilon;
.nf = ((.guess - .delta)|f)
| if .nf < .f0
then .f0 = .nf
| .guess = .guess - .delta
else .nf = ((.guess + .delta)|f)
| if .nf < .f0
then .f0 = .nf
| .guess = .guess + .delta
else .factor = 0.5
end
end
| .delta = .delta * .factor )
| .guess;
((solve(f; 0.5; 0) * 1e10)|round / 1e10 ) as $r
| ((((12 * $r)|exp) * 1e8)|round / 1e8) as $R0
| "r = \($r), R0 = \($R0)"
- Output:
r = 0.1123021757, R0 = 3.84827928
Julia
using LsqFit
const K = 7_800_000_000 # approximate world population
const n0 = 27 # starting at day 0 with 27 Chinese cases
""" The model for logistic regression with a given r """
@. model(t, r) = (n0 * exp(r * t)) / (( 1 + n0 * (exp(r * t) - 1) / K))
# Daily world totals of covid cases, all countries
ydata = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652,
]
tdata = collect(LinRange(0.0, 96, 97))
# starting approximation for r of 1/2
rparam = [0.5]
fit = curve_fit(model, tdata, ydata, rparam)
# Our answer for r given the world data and simplistic model
r = fit.param
println("The logistic curve r for the world data is: ", r)
println("The confidence interval at 5% significance is: ",
confidence_interval(fit, 0.05))
println("Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ ", exp(12r[1]))
- Output:
The logistic curve r for the world data is: [0.11230217572265622] The confidence interval at 5% significance is: [(0.11199074156706985, 0.11261360987824258)] Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ 3.8482792820761063
Kotlin
import kotlin.math.exp
const val K = 7.8e9
const val n0 = 27
val actual = arrayOf(
27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0, 60.0,
61.0, 61.0, 66.0, 83.0, 219.0, 239.0, 392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0,
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0, 20615.0, 24522.0, 28273.0,
31491.0, 34933.0, 37552.0, 40540.0, 43105.0, 45177.0, 60328.0, 64543.0, 67103.0,
69265.0, 71332.0, 73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0,
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0, 90664.0, 93077.0,
95316.0, 98172.0, 102133.0, 105824.0, 109695.0, 114232.0, 118610.0, 125497.0,
133852.0, 143227.0, 151367.0, 167418.0, 180096.0, 194836.0, 213150.0, 242364.0,
271106.0, 305117.0, 338133.0, 377918.0, 416845.0, 468049.0, 527767.0, 591704.0,
656866.0, 715353.0, 777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0
)
fun f(r: Double): Double {
var sq = 0.0
val len = actual.size
for (i in 0 until len) {
val eri = exp(r * i)
val guess = (n0 * eri) / (1.0 + n0 * (eri - 1.0) / K)
val diff = guess - actual[i]
sq += diff * diff
}
return sq
}
fun solve(fn: (Double) -> Double, guess: Double = 0.5, epsilon: Double = 0.0): Double {
var guess2 = guess
var delta = if (guess2 != 0.0) guess2 else 1.0
var f0 = fn(guess2)
var factor = 2.0
while (delta > epsilon && guess2 != guess2 - delta) {
var nf = fn(guess2 - delta)
if (nf < f0) {
f0 = nf
guess2 -= delta
} else {
nf = fn(guess2 + delta)
if (nf < f0) {
f0 = nf
guess2 += delta
} else {
factor = 0.5
}
}
delta *= factor
}
return guess2
}
fun main() {
val r = solve(::f)
val r0 = exp(12.0 * r)
println("r = $r, R0 = $r0")
}
- Output:
r = 0.11230217569755005, R0 = 3.8482792809167194
Nim
import math, sugar
const
K = 7.8e9
N0 = 27
Actual = [ 27.0, 27.0, 27.0, 44.0, 44.0, 59.0, 59.0,
59.0, 59.0, 59.0, 59.0, 59.0, 59.0, 60.0,
60.0, 61.0, 61.0, 66.0, 83.0, 219.0, 239.0,
392.0, 534.0, 631.0, 897.0, 1350.0, 2023.0, 2820.0,
4587.0, 6067.0, 7823.0, 9826.0, 11946.0, 14554.0, 17372.0,
20615.0, 24522.0, 28273.0, 31491.0, 34933.0, 37552.0, 40540.0,
43105.0, 45177.0, 60328.0, 64543.0, 67103.0, 69265.0, 71332.0,
73327.0, 75191.0, 75723.0, 76719.0, 77804.0, 78812.0, 79339.0,
80132.0, 80995.0, 82101.0, 83365.0, 85203.0, 87024.0, 89068.0,
90664.0, 93077.0, 95316.0, 98172.0, 102133.0, 105824.0, 109695.0,
114232.0, 118610.0, 125497.0, 133852.0, 143227.0, 151367.0, 167418.0,
180096.0, 194836.0, 213150.0, 242364.0, 271106.0, 305117.0, 338133.0,
377918.0, 416845.0, 468049.0, 527767.0, 591704.0, 656866.0, 715353.0,
777796.0, 851308.0, 928436.0, 1000249.0, 1082054.0, 1174652.0]
type Func = (float) -> float
func f(r: float): float =
for i in 0..Actual.high:
let
eri = exp(r * i.toFloat)
guess = (N0 * eri) / (1 + N0 * (eri - 1) / K)
diff = guess - Actual[i]
result += diff * diff
func solve(fn: Func; guess = 0.5, epsilon = 0.0): float =
result = guess
var delta = if result != 0: result else: 1.0
var f0 = fn(result)
var factor = 2.0
while delta > epsilon and result != result - delta:
var nf = fn(result - delta)
if nf < f0:
f0 = nf
result -= delta
else:
nf = fn(result + delta)
if nf < f0:
f0 = nf
result += delta
else:
factor = 0.5
delta *= factor
let r = f.solve()
let r0 = exp(12 * r)
echo "r = ", r, ", R0 = ", r0
- Output:
r = 0.11230217569755, R0 = 3.848279280916719
Perl
use strict;
use warnings;
my $K = 7_800_000_000; # population
my $n0 = 27; # cases @ day 0
my @y = (
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59,
60, 60, 61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350,
2023, 2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, 31491,
34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, 133852, 143227, 151367, 167418, 180096,
194836, 213150, 242364, 271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, 656866, 715353,
777796, 851308, 928436,1000249,1082054,1174652
);
sub logistic_func {
my($r) = @_;
my $sq = 0;
for my $i (0 .. @y-1) {
my $eri = exp($r * $i);
my $dst = ($n0 * $eri) / (1 + $n0 * ($eri-1) / $K) - $y[$i];
$sq = $sq + $dst**2;
}
$sq
}
sub solve {
my($fn, $guess, $epsilon) = @_;
my($nfm,$nfp);
my $f0 = &$fn($guess);
my $delta = $guess;
my $factor = 2;
while ($delta > $epsilon) {
($nfm = &$fn($guess - $delta)) < $f0 ?
($f0 = $nfm, $guess -= $delta, $delta *= $factor)
:
($nfp = &$fn($guess + $delta)) < $f0 ?
($f0 = $nfp, $guess += $delta, $delta *= $factor)
:
$delta /= $factor
}
$guess
}
my $r = solve(\&logistic_func, 0.5, 0);
my $R0 = exp(12 * $r);
printf "r = %%(%.3f), R0 = %%(%.3f)\n", $r, $R0;
- Output:
r = %(0.112), R0 = %(3.848)
Phix
Simplified, my interpretation of shift-cutting (from wp:Non-linear_least_squares)
-- demo\rosetta\Curve_fit.exw with javascript_semantics constant K = 7_800_000_000, -- approx world population n0 = 27, -- number of cases at day 0 actual = { 27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60, 61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497, 133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364, 271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704, 656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652 } function f(atom r) atom sq = 0 for i=1 to length(actual) do atom eri = exp(r*(i-1)), guess = (n0*eri)/(1+n0*(eri-1)/K), diff = guess-actual[i] sq += diff*diff end for return sq end function function solve(integer f, atom guess=0.5, epsilon=0) atom f0 = f(guess), delta = iff(guess?guess:1), factor = 2 -- double until f0 best, then -- halve until delta<=epsilon -- or we hit precision limit while delta>epsilon -- (predefined limit) and guess!=guess-delta do -- (precision limit) atom nf = f(guess-delta) if nf<f0 then f0 = nf guess -= delta else nf = f(guess+delta) if nf<f0 then f0 = nf guess += delta else factor = 0.5 end if end if delta *= factor end while return guess end function atom r = solve(f), R0 = exp(12 * r) printf(1,"r = %.10f, R0 = %.8f\n",{r,R0})
- Output:
r = 0.1123021757, R0 = 3.84827928
Python
Uses NumPy/SciPy's optimize package.
import numpy as np
import scipy.optimize as opt
n0, K = 27, 7_800_000_000
def f(t, r):
return (n0 * np.exp(r * t)) / (( 1 + n0 * (np.exp(r * t) - 1) / K))
y = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652,
]
x = np.linspace(0.0, 96, 97)
r, cov = opt.curve_fit(f, x, y, [0.5])
# Our answer for r given the world data and simplistic model
print("The r for the world Covid-19 data is:", r,
", with covariance of", cov)
print("The calculated R0 is then", np.exp(12 * r))
- Output:
The r for the world Covid-19 data is: [0.11230218] , with covariance of [[2.46164331e-08]] The calculated R0 is then [3.8482793]
R
library(minpack.lm)
K <- 7800000000 # approximate world population
n0 <- 27 # starting at day 0 with 27 Chinese cases
x <- seq(from=0.0, 96.0, length=97)
# The model for logistic regression with a given r0
getPred <- function(par, xx) (n0 * exp(par$r * xx)) / (( 1 + n0 * (exp(par$r * xx) - 1) / K))
residFun <- function(p, observed, xx) observed - getPred(p, xx)
parStart <- list(r=0.5)
# Daily world totals of covid cases, all countries
observed <- c(
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652)
nls.out <- nls.lm(par=parStart, fn=residFun, observed=observed, xx = x)
cat("The r for the model is: ", coef(nls.out), "\n")
cat("Therefore, R0 is approximately: ", exp(12 * coef(nls.out)))
- Output:
The r for the model is: 0.1123022 Therefore, R0 is approximately: 3.848279
Raku
Original task numbers of cases per day as of April 05 plus updated, as of today, May 11 numbers.
my $K = 7_800_000_000; # population
my $n0 = 27; # cases @ day 0
my @Apr05 = <
27 27 27 44 44 59 59 59 59 59
59 59 59 60 60 61 61 66 83 219
239 392 534 631 897 1350 2023 2820 4587 6067
7823 9826 11946 14554 17372 20615 24522 28273 31491 34933
37552 40540 43105 45177 60328 64543 67103 69265 71332 73327
75191 75723 76719 77804 78812 79339 80132 80995 82101 83365
85203 87024 89068 90664 93077 95316 98172 102133 105824 109695
114232 118610 125497 133852 143227 151367 167418 180096 194836 213150
242364 271106 305117 338133 377918 416845 468049 527767 591704 656866
715353 777796 851308 928436 1000249 1082054 1174652
>;
my @May11 = <
27 27 27 44 44 59 59 59 59 59
59 59 59 60 60 61 61 66 83 219
239 392 534 631 897 1350 2023 2820 4587 6067
7823 9826 11946 14554 17372 20615 24522 28273 31491 34933
37552 40540 43105 45177 60328 64543 67103 69265 71332 73327
75191 75723 76719 77804 78812 79339 80132 80995 82101 83365
85203 87026 89068 90865 93077 95316 98172 102133 105824 109695
114235 118613 125497 133853 143259 154774 167418 180094 194843 213149
242374 271116 305235 338235 377968 416881 468092 527839 591796 656834
715377 777187 851587 928491 1006063 1087822 1174306 1245601 1316988 1391881
1476792 1563819 1653160 1734868 1807256 1873639 1953786 2033745 2117297 2199460
2278484 2350993 2427353 2513399 2579823 2657910 2731217 2832750 2915977 2981542
3054404 3131487 3216467 3308341 3389459 3468047 3545486 3624789 3714816 3809262
3899379 3986931 4063525
>;
sub logistic-func ($rate, @y) {
my $sq = 0;
for ^@y -> $time {
my $ert = exp($rate * $time);
my $δt = ($n0 * $ert) / (1 + $n0 * ($ert-1) / $K) - @y[$time];
$sq += $δt²;
}
$sq
}
sub solve (&f, $guess, \ε, @y) {
my $fₙ-minus;
my $fₙ-plus;
my $rate = $guess;
my $fₙ = f $rate, @y;
my $Δ = $rate;
my $factor = 2;
while $Δ > ε {
($fₙ-minus = f $rate - $Δ, @y) < $fₙ ??
do {
$fₙ = $fₙ-minus;
$rate -= $Δ;
$Δ *= $factor;
} !!
($fₙ-plus = f $rate + $Δ, @y) < $fₙ ??
do {
$fₙ = $fₙ-plus;
$rate += $Δ;
$Δ *= $factor;
} !!
$Δ /= $factor
}
$rate
}
for @Apr05, 'Dec 31 - Apr 5',
@May11, 'Dec 31 - May 11' -> @y, $period {
my $rate = solve(&logistic-func, 0.5, 0, @y);
my $R0 = exp(12 * $rate);
say "\nFor period: $period";
say "Instantaneous rate of growth: r = " , $rate.fmt('%08f');
say "Reproductive rate: R0 = ", $R0.fmt('%08f');
}
- Output:
For period: Dec 31 - Apr 5 Instantaneous rate of growth: r = 0.112302 Reproductive rate: R0 = 3.848279 For period: Dec 31 - May 11 Instantaneous rate of growth: r = 0.093328 Reproductive rate: R0 = 3.064672
Ruby
K = 7.9e9
N0 = 27
ACTUAL = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023, 2820,
4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615, 24522, 28273,
31491, 34933, 37552, 40540, 43105, 45177, 60328, 64543, 67103,
69265, 71332, 73327, 75191, 75723, 76719, 77804, 78812, 79339,
80132, 80995, 82101, 83365, 85203, 87024, 89068, 90664, 93077,
95316, 98172, 102133, 105824, 109695, 114232, 118610, 125497,
133852, 143227, 151367, 167418, 180096, 194836, 213150, 242364,
271106, 305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054, 1174652
]
def f(r)
sq = 0.0
len = ACTUAL.length
for i in 1 .. len
j = i - 1
eri = Math.exp(r * j)
guess = (N0 * eri) / (1 + N0 * (eri - 1.0) / K)
diff = guess - ACTUAL[j]
sq += diff * diff
end
return sq
end
def solve(fn, guess=0.5, epsilon=0.0)
delta = guess ? guess : 1.0
f0 = send(fn, guess)
factor = 2.0
while delta > epsilon and guess != guess - delta
nf = send(fn, guess - delta)
if nf < f0 then
f0 = nf
guess -= delta
else
nf = send(fn, guess + delta)
if nf < f0 then
f0 = nf
guess += delta
else
factor = 0.5
end
end
delta *= factor
end
return guess
end
def main
r = solve(:f)
r0 = Math.exp(12.0 * r)
print "r = ", r, ", R0 = ", r0, "\n"
end
main()
- Output:
r = 0.11230215850810021, R0 = 3.8482784871191575
V (Vlang)
import math
const (
k = 7800000000 // approx world population
n0 = 27 // number of cases at day 0
y = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652
]
)
fn f(r f64) f64 {
mut sq := 0.0
for i in 0..y.len {
eri := math.exp(r * i)
dst := (n0 * eri) / (1 + n0 * (eri - 1) / k) -y[i]
sq = sq + dst * dst
}
return sq
}
fn solve(fu fn(f64)f64, g f64, epsilon int) f64 {
mut guess := g
mut f0 := fu(guess)
mut delta := guess
mut factor := 2.0
for delta > epsilon {
mut nf := fu(guess - delta)
if nf < f0 {
f0 = nf
guess -= delta
} else {
nf = fu(guess + delta)
if nf < f0 {
f0 = nf
guess += delta
} else {
factor = 0.5
}
}
delta *= factor
}
return guess
}
fn main() {
r := math.round(solve(f, 0.5, 0) * 1e10) / 1e10
r0 := math.round(math.exp(12 * r) * 1e8) / 1e8
println("r = $r, R0 = $r0")
}
- Output:
r = 0.1123021757, R0 = 3.84827928
Wren
var K = 7800000000 // approx world population
var n0 = 27 // number of cases at day 0
var y = [
27, 27, 27, 44, 44, 59, 59, 59, 59, 59, 59, 59, 59, 60, 60,
61, 61, 66, 83, 219, 239, 392, 534, 631, 897, 1350, 2023,
2820, 4587, 6067, 7823, 9826, 11946, 14554, 17372, 20615,
24522, 28273, 31491, 34933, 37552, 40540, 43105, 45177,
60328, 64543, 67103, 69265, 71332, 73327, 75191, 75723,
76719, 77804, 78812, 79339, 80132, 80995, 82101, 83365,
85203, 87024, 89068, 90664, 93077, 95316, 98172, 102133,
105824, 109695, 114232, 118610, 125497, 133852, 143227,
151367, 167418, 180096, 194836, 213150, 242364, 271106,
305117, 338133, 377918, 416845, 468049, 527767, 591704,
656866, 715353, 777796, 851308, 928436, 1000249, 1082054,
1174652
]
var f = Fn.new { |r|
var sq = 0
for (i in 0...y.count) {
var eri = (r*i).exp
var dst = (n0*eri)/(1+n0*(eri-1)/K) - y[i]
sq = sq + dst * dst
}
return sq
}
var solve = Fn.new { |f, guess, epsilon|
var f0 = f.call(guess)
var delta = guess
var factor = 2 // double until f0 best then halve until delta <= epsilon
while (delta > epsilon) {
var nf = f.call(guess - delta)
if (nf < f0) {
f0 = nf
guess = guess - delta
} else {
nf = f.call(guess + delta)
if (nf < f0) {
f0 = nf
guess = guess + delta
} else {
factor = 0.5
}
}
delta = delta * factor
}
return guess
}
var r = (solve.call(f, 0.5, 0) * 1e10).round / 1e10
var R0 = ((12 * r).exp * 1e8).round / 1e8
System.print("r = %(r), R0 = %(R0)")
- Output:
r = 0.1123021757, R0 = 3.84827928
XPL0
def K = 7.8e9;
def N0 = 27.;
real Actual;
int ActualSize;
func real Func(R);
real R;
real Sq, Eri, Guess, Diff;
int I;
[Sq:= 0.;
for I:= 0 to ActualSize-1 do
[Eri:= Exp(R * float(I));
Guess:= N0 * Eri / (1. + N0 * (Eri-1.) / K);
Diff:= Guess - Actual(I);
Sq:= Sq + Diff*Diff;
];
return Sq;
];
func real Solve(Guess, Epsilon);
real Guess, Epsilon;
real Delta, F0, Factor, NF;
[Delta:= if Guess # 0. then Guess else 1.;
F0:= Func(Guess);
Factor:= 2.;
while Delta > Epsilon and Guess # Guess-Delta do
[NF:= Func(Guess - Delta);
if NF < F0 then
[F0:= NF;
Guess:= Guess - Delta;
]
else
[NF:= Func(Guess + Delta);
if NF < F0 then
[F0:= NF;
Guess:= Guess + Delta;
]
else Factor:= 0.5;
];
Delta:= Delta * Factor;
];
return Guess;
];
real R, R0;
[Actual:= [
27., 27., 27., 44., 44., 59., 59., 59., 59., 59., 59., 59., 59., 60., 60.,
61., 61., 66., 83., 219., 239., 392., 534., 631., 897., 1350., 2023., 2820.,
4587., 6067., 7823., 9826., 11946., 14554., 17372., 20615., 24522., 28273.,
31491., 34933., 37552., 40540., 43105., 45177., 60328., 64543., 67103.,
69265., 71332., 73327., 75191., 75723., 76719., 77804., 78812., 79339.,
80132., 80995., 82101., 83365., 85203., 87024., 89068., 90664., 93077.,
95316., 98172., 102133., 105824., 109695., 114232., 118610., 125497.,
133852., 143227., 151367., 167418., 180096., 194836., 213150., 242364.,
271106., 305117., 338133., 377918., 416845., 468049., 527767., 591704.,
656866., 715353., 777796., 851308., 928436., 1000249., 1082054., 1174652.];
ActualSize:= 97;
R:= Solve(0.5, 0.0);
R0:= Exp(12.*R);
Text(0, "R = "); RlOut(0, R); CrLf(0);
Text(0, "R0 = "); RlOut(0, R0); CrLf(0);
]
- Output:
R = 0.11230 R0 = 3.84828