Jump to content

Logistic curve fitting in epidemiology

From Rosetta Code
Task
Logistic curve fitting in epidemiology
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

Translation of: C++
-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

Translation of: C++
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

Translation of: 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++

Translation of: Phix
#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

Translation of: C++
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

Translation of: Phix
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

Library: lm

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

Translation of: Kotlin
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

Translation of: Wren
Works with: 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

Translation of: D
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

Translation of: Kotlin
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

Translation of: Phix
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

Works with: Rakudo version 2020.05
Translation of: Perl

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

Translation of: C++
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)

Translation of: wren
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

Translation of: Phix
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

Translation of: C
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
Retrieved from "https://rosettacode.org/wiki/Logistic_curve_fitting_in_epidemiology?oldid=354163"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information