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

12,599 bytes added ,  5 months ago
m
→‎{{header|Wren}}: Changed to Wren S/H
(Added 11l)
m (→‎{{header|Wren}}: Changed to Wren S/H)
 
(11 intermediate revisions by 10 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 77:
{{trans|C++}}
 
<langsyntaxhighlight lang="11l">-V K = 7.8e9
-V n0 = 27
-V actual = [
Line 122:
V r = solve(f)
V R0 = exp(12 * r)
print(‘r = #.6, R0 = #.6’.format(r, R0))</langsyntaxhighlight>
 
{{out}}
Line 131:
=={{header|Ada}}==
{{trans|C++}}
<langsyntaxhighlight Adalang="ada">with Ada.Text_Io;
with Ada.Numerics.Generic_Elementary_Functions;
 
Line 221:
Put (", R0 = "); Real_Io.Put (R0, Exp => 0, Aft => 7);
New_Line;
end Curve_Fitting;</langsyntaxhighlight>
{{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 289 ⟶ 362:
printf("r = %f, R0 = %f\n", r, R0);
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>r = 0.112302, R0 = 3.848279</pre>
Line 295 ⟶ 368:
=={{header|C++}}==
{{trans|Phix}}
<langsyntaxhighlight lang="cpp">#include <cmath>
#include <functional>
#include <iostream>
Line 351 ⟶ 424:
std::cout << "r = " << r << ", R0 = " << R0 << '\n';
return 0;
}</langsyntaxhighlight>
 
{{out}}
Line 360 ⟶ 433:
=={{header|D}}==
{{trans|C++}}
<langsyntaxhighlight lang="d">import std.math;
import std.stdio;
 
Line 416 ⟶ 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 478 ⟶ 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 488 ⟶ 621:
=={{header|Java}}==
{{trans|Kotlin}}
<langsyntaxhighlight lang="java">import java.util.List;
import java.util.function.Function;
 
Line 561 ⟶ 694:
System.out.printf("r = %.16f, R0 = %.16f\n", r, r0);
}
}</langsyntaxhighlight>
{{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 602 ⟶ 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 611 ⟶ 805:
=={{header|Kotlin}}==
{{trans|D}}
<langsyntaxhighlight lang="scala">import kotlin.math.exp
 
const val K = 7.8e9
Line 672 ⟶ 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 726 ⟶ 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 732 ⟶ 988:
=={{header|Phix}}==
Simplified, my interpretation of shift-cutting (from [[wp:Non-linear_least_squares]])
<!--<langsyntaxhighlight Phixlang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Curve_fit.exw</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
Line 791 ⟶ 1,047:
<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>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 799 ⟶ 1,055:
=={{header|Python}}==
Uses NumPy/SciPy's optimize package.
<langsyntaxhighlight lang="python">import numpy as np
import scipy.optimize as opt
 
Line 829 ⟶ 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 836 ⟶ 1,092:
 
=={{header|R}}==
<syntaxhighlight lang="r">
<lang R>
library(minpack.lm)
 
Line 871 ⟶ 1,127:
 
cat("Therefore, R0 is approximately: ", exp(12 * coef(nls.out)))
</langsyntaxhighlight>{{out}}
<pre>
The r for the model is: 0.1123022
Line 883 ⟶ 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 959 ⟶ 1,215:
say "Reproductive rate: R0 = ", $R0.fmt('%08f');
}
</syntaxhighlight>
</lang>
{{out}}
<pre>For period: Dec 31 - Apr 5
Line 971 ⟶ 1,227:
=={{header|Ruby}}==
{{trans|C++}}
<langsyntaxhighlight lang="ruby">K = 7.9e9
N0 = 27
ACTUAL = [
Line 1,029 ⟶ 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 1,059 ⟶ 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 1,090 ⟶ 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,476

edits

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