Rosetta Code

Logistic curve fitting in epidemiology: Difference between revisions

Page Discussion

Logistic curve fitting in epidemiology (view source)

Revision as of 15:51, 17 December 2023

35,563 bytes added ,  5 months ago
m
→‎{{header|Wren}}: Changed to Wren S/H
m (elided the 1st task, corrected 2 misspellings, used a bigger font to show the super- and subscripts, (fine print), used a consistent font for variables in/from expressions.corrected 2 uses of n sub oh-->n sub zero, moved 2 definitions closer to their use.)
m (→‎{{header|Wren}}: Changed to Wren S/H)
 
(16 intermediate revisions by 14 users not shown)
Line 1:
{{task|Probability and statistics}}
{{draft task}}
 
The least-squares method (see references below) in statistics is used to fit data
Line 73:
:;*[[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++}}
<langsyntaxhighlight lang="c">#include <math.h>
#include <stdio.h>
 
Line 138 ⟶ 362:
printf("r = %f, R0 = %f\n", r, R0);
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>r = 0.112302, R0 = 3.848279</pre>
Line 144 ⟶ 368:
=={{header|C++}}==
{{trans|Phix}}
<langsyntaxhighlight lang="cpp">#include <cmath>
#include <functional>
#include <iostream>
Line 200 ⟶ 424:
std::cout << "r = " << r << ", R0 = " << R0 << '\n';
return 0;
}</langsyntaxhighlight>
 
{{out}}
Line 209 ⟶ 433:
=={{header|D}}==
{{trans|C++}}
<langsyntaxhighlight lang="d">import std.math;
import std.stdio;
 
Line 265 ⟶ 489:
double R0 = exp(12 * r);
writeln("r = ", r, ", R0 = ", R0);
}</langsyntaxhighlight>
{{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}}==
{{libheader|lm}}
This uses the Levenberg-Marquardt method.
<langsyntaxhighlight lang="go">package main
 
import (
Line 327 ⟶ 611:
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))
}</langsyntaxhighlight>
 
{{out}}
Line 334 ⟶ 618:
R0 is then approximately equal to 3.8482793
</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}}==
<langsyntaxhighlight lang="julia">using LsqFit
 
const K = 7_800_000_000 # approximate world population
Line 372 ⟶ 796:
confidence_interval(fit, 0.05))
println("Since R0 ≈ exp(G * r), and G ≈ 12 days, R0 ≈ ", exp(12r[1]))
</langsyntaxhighlight>{{out}}
<pre>
The logistic curve r for the world data is: [0.11230217572265622]
Line 381 ⟶ 805:
=={{header|Kotlin}}==
{{trans|D}}
<langsyntaxhighlight lang="scala">import kotlin.math.exp
 
const val K = 7.8e9
Line 442 ⟶ 866:
val r0 = exp(12.0 * r)
println("r = $r, R0 = $r0")
}</langsyntaxhighlight>
{{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}}==
{{trans|Phix}}
<langsyntaxhighlight lang="perl">use strict;
use warnings;
 
Line 496 ⟶ 982:
my $r = solve(\&logistic_func, 0.5, 0);
my $R0 = exp(12 * $r);
printf "r = %%(%.3f), R0 = %%(%.3f)\n", $r, $R0;</langsyntaxhighlight>
{{out}}
<pre>r = %(0.112), R0 = %(3.848)</pre>
Line 502 ⟶ 988:
=={{header|Phix}}==
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
while delta>epsilon -- (predefinedhalve limit)until delta&lt;=epsilon
and guess!=guess-delta do -- ( or we hit precision limit)</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">delta</span><span style="color: #0000FF;">></span><span style="color: #000000;">epsilon</span> <span style="color: #000080;font-style:italic;">-- (predefined limit)</span>
atom nf = f(guess-delta)
<span style="color: #008080;">and</span> <span style="color: #000000;">guess</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">-</span><span style="color: #000000;">delta</span> <span style="color: #008080;">do</span> <span style="color: #000080;font-style:italic;">-- (precision limit)</span>
if nf<f0 then
<span style="color: #004080;">atom</span> <span style="color: #000000;">nf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">-</span><span style="color: #000000;">delta</span><span style="color: #0000FF;">)</span>
f0 = nf
<span style="color: #008080;">if</span> <span style="color: #000000;">nf</span><span style="color: #0000FF;"><</span><span style="color: #000000;">f0</span> <span style="color: #008080;">then</span>
guess -= delta
<span style="color: #000000;">f0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nf</span>
else
<span style="color: #000000;">guess</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">delta</span>
nf = f(guess+delta)
<span style="color: if nf#008080;">else<f0 then/span>
<span style="color: #000000;">nf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">guess</span><span style="color: #0000FF;">+</span><span style="color: #000000;">delta</span><span style="color: #0000FF;">)</span>
f0 = nf
<span style="color: #008080;">if</span> <span style="color: #000000;">nf</span><span style="color: #0000FF;"><</span><span style="color: #000000;">f0</span> <span style="color: #008080;">then</span>
guess += delta
<span style="color: #000000;">f0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nf</span>
else
<span style="color: #000000;">guess</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">delta</span>
factor = 0.5
end if<span style="color: #008080;">else</span>
<span style="color: #000000;">factor</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.5</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
delta *= factor
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end while
<span style="color: #000000;">delta</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">factor</span>
return guess
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">guess</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
atom r = solve(f),
R0 = exp(12 * r)
<span style="color: #004080;">atom</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">solve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">),</span>
printf(1,"r = %.10f, R0 = %.8f\n",{r,R0})</lang>
<span style="color: #000000;">R0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">12</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"r = %.10f, R0 = %.8f\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R0</span><span style="color: #0000FF;">})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 566 ⟶ 1,055:
=={{header|Python}}==
Uses NumPy/SciPy's optimize package.
<langsyntaxhighlight lang="python">import numpy as np
import scipy.optimize as opt
 
Line 596 ⟶ 1,085:
", with covariance of", cov)
print("The calculated R0 is then", np.exp(12 * r))
</langsyntaxhighlight>{{out}}
<pre>
The r for the world Covid-19 data is: [0.11230218] , with covariance of [[2.46164331e-08]]
Line 603 ⟶ 1,092:
 
=={{header|R}}==
<syntaxhighlight lang="r">
<lang R>
library(minpack.lm)
 
Line 638 ⟶ 1,127:
 
cat("Therefore, R0 is approximately: ", exp(12 * coef(nls.out)))
</langsyntaxhighlight>{{out}}
<pre>
The r for the model is: 0.1123022
Line 650 ⟶ 1,139:
Original task numbers of cases per day as of April 05 plus updated, as of today, May 11 numbers.
 
<syntaxhighlight lang="raku" perl6line>my $K = 7_800_000_000; # population
my $n0 = 27; # cases @ day 0
 
Line 726 ⟶ 1,215:
say "Reproductive rate: R0 = ", $R0.fmt('%08f');
}
</syntaxhighlight>
</lang>
{{out}}
<pre>For period: Dec 31 - Apr 5
Line 738 ⟶ 1,227:
=={{header|Ruby}}==
{{trans|C++}}
<langsyntaxhighlight lang="ruby">K = 7.9e9
N0 = 27
ACTUAL = [
Line 796 ⟶ 1,285:
end
 
main()</langsyntaxhighlight>
{{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}}==
{{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
 
Line 826 ⟶ 1,374:
var sq = 0
for (i in 0...y.count) {
var eri = Math.exp(r*i).exp
var dst = (n0*eri)/(1+n0*(eri-1)/K) - y[i]
sq = sq + dst * dst
Line 857 ⟶ 1,405:
 
var r = (solve.call(f, 0.5, 0) * 1e10).round / 1e10
var R0 = (Math.exp(12 * r).exp * 1e8).round / 1e8
System.print("r = %(r), R0 = %(R0)")</langsyntaxhighlight>
 
{{out}}
<pre>
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>
PureFox
9,479

edits

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