BLog

ImprintImpressum
PrivacyDatenschutz
DisclaimerHaftung
Downloads 

Zeitreihenanalyse zur Ausbreitung von Covid-19 ausserhalb Chinas

Im vorigen BLog-Artikel hatte ich die Ausbreitung des n-Coronavirus-2019 in China vom 8. Februar bis zum 24. Februar durch Fortschreiben der Kurvenanpassung einer modifizierten logistischen Funktion an die von der WHO berichteten Infektionsfälle begleitet. Dabei hatte bereits die erste Kurvenanpassung die Ausbreitung, nämlich Wendepunkt und Grenzwert des logistischen Wachstums erstaunlich gut vorhergesagt, und ich beendete die betreffende Analyse beim Erreichen des Grenzwertes. Seitdem hat sich die Situation in China stabilisiert, aber nun gibt es mehrere andere heiße Herde, darunter Italien und Deutschland, und in Brasilien, dort wo ich als Deutscher lebe, geht es nun auch los. Ich liebe Italien und die Italiener, mein Herz schlägt immer noch für Deutschland und in Brasilien sowieso. Können also Zeitreihenanalysen mittels fortgeschriebener Kurvenanpassungen auch für die Covid-19-Ausbrüche in den genannten Ländern valide Vorhersagen liefern, so wie für China?

Die Daten für die China-Zeitreihe hatte ich noch mühsam händisch aus den WHO-Situationsberichten zusammengetragen, und in mein Program CVA zur Auswertung von Meßkurven zwecks Darstellung und Kurvenanpassung eingespeist. Mittlerweile unterhält das Center for Systems Science and Engineering (CSSE) an der Johns Hopkins University in ihrem GitHub-Repository tagesaktualisierte Daten zur Ausbreitung von Covid-19 in vielen Ländern und Regionen, und man kann sich die Daten als CSV-Tabelle herunterladen. Die Zeitreihen für die jeweiligen Länder/Regionen sind zeilenweise angeordnet.

Ich habe ein Kommandozeilen-Programm für FreeBSD, macOS und Linux geschrieben, mit dem man die Zeitreihe für spezifizierte Länder bzw. Regionen automatisch extrahieren und in eine spaltenorientierte TSV-Datei exportieren kann, die man dann direkt mit einem Analyse-Programm seiner Wahl (ich benutze CVA) zwecks weiterer Auswertung und graphischer Darstellung öffnen kann - Excel geht auch. Im selben Atemzug wird auch eine Kurvenanpassung an die logistische Funktion versucht, und der Verlauf der angepaßten Kurve wird in einer separaten Spalte selbiger TSV-Datei gespeichert. Der Quellcode und die Xcode-Projekt-Datei befinden sich in meinem GitHub-Repository.

Für die logistische Funktion bieten sich zwei Darstellungsformen an:

  1. die generische Form:  a/(1 + exp(-b·(x - c)))
  2. die alternative Form:  a·0.5/(0.5 + (a - 0.5)·exp(-b·a·(x - c)))

Die logistische Funktion hat eine hochinteressante Eigenschaft, sie ist nämlich um den Wendepunkt symmetrisch, und sobald man den Wendepunkt sicher eingefangen hat, weiß man mit an Sicherheit grenzender Wahrscheinlichkeit bereits wie es ausgehen wird.

Für die Corona-Zeitreihenanalyse zur Ausbreitung in China hatte ich die um ein lineares Glied modifizierte alternative From verwendet, weil sie sich bei der Kurvenanpassung insbesondere in der Anfangsphase als stabiler erwies. Die generische Form hat allerdings den Vorteil, daß die Lage des Wendepunktes direkt aus dem Parameter c hervorgeht. Deshalb versucht man eine Anpassung zunächst mit der Form 1, und falls diese absurde Ergebnisse liefern sollte, kann man es noch einmal mit der alternativen Form 2 versuchen.

//  xcssecovid.c
//  xcssecovid
//
//  Created by Dr. Rolf Jansen on 2020-03-22.
//  Copyright © 2020 Dr. Rolf Jansen. All rights reserved.
//
//  Redistribution and use in source and binary forms, with or without modification,
//  are permitted provided that the following conditions are met:
//
//  1. Redistributions of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//  2. Redistributions in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
//  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
//  INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
//  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
//  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
//  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
//  OF THE POSSIBILITY OF SUCH DAMAGE.
//
//  Usage:
//
//  1. Compile this file on either of FreeBSD or macOS:
//
//     clang -g0 -O3 -march=native xcssecovid.c -Wno-parentheses -lm -o xcssecovid
//
//  2. Download the daily updated time series of confirmed Covid-19 cases
//     from CSSE's (at Johns Hopkins University) GitHub site CSSE COVID-19 Dataset
//     https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data
//
//     fetch https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
//
//  3. Extract the time series row of a given country starting on 2020-01-21 as day 0
//     and write it out together with the country's cases and a simulated curve-fitted
//     logistic function into a three-column TSV output file:
//
//     ./xcssecovid Germany time_series_covid19_confirmed_global.csv Germany.tsv
//
//  4. Open the TSV file with your favorite graphing and/or data analysis application.


#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <sys/stat.h>


typedef long double ldouble;
typedef ldouble    *Vector;
typedef ldouble   **Matrix;
typedef int        *Index;


static ldouble sqr(ldouble x)
{
   return x*x;
}

ldouble LogisticFunctionGen(ldouble t, Vector A)
{  // generic form: https://en.wikipedia.org/wiki/Logistic_function - parameters A[3] to A[5] = 0
   return A[0]/(1 + exp(-A[1]*(t - A[2]))) - (A[3] + A[4]*t + A[5]*sqr(t));
}


ldouble LogisticFunctionAlt(ldouble t, Vector A)
{  // alternate form: https://de.wikipedia.org/wiki/Logistische_Funktion#Weitere_Darstellungen - parameters A[3] to A[5] = 0
   return A[0]*0.5/(0.5 + (A[0] - 0.5)*exp(-A[1]*A[0]*(t - A[2]))) - (A[3] + A[4]*t + A[5]*sqr(t));
}

ldouble curveFit(int n, Vector T, Vector C, Vector S,
                 int k, Vector A, Vector dA,
                 ldouble (*model)(ldouble t, Vector A));

static inline size_t collen(const char *col)
{
   if (!col || !*col)
      return 0;

   size_t l;
   for (l = 0; col[l] && col[l] != ','; l++)
      ;
   return l;
}

int main(int argc, char *const argv[])
{
   FILE  *csv, *tsv;
   char  *country = argv[1];
   size_t cl = strlen(country);

   if (argc == 4 && *country)
      if (csv = (*(uint16_t *)argv[2] == *(uint16_t *)"-")
                ? stdin
                : fopen(argv[2], "r"))
      {
         if (tsv = (*(uint16_t *)argv[3] == *(uint16_t *)"-")
                   ? stdout
                   : fopen(argv[3], "w"))
         {
            int     i, m = 0, n = 0;
            char   *line;
            char    d[65536];
            ldouble t[366], c[366], l[366];  // 2020 is a leap year - t[0] = 2020-01-01; t[365] = 2020-12-31

            // day 1 of the CSSE time series is 2020-01-22, hence day 0 is 2020-01-21, i.e. t[20]
            for (i = 0; i < 366; i++)
               t[i] = i - 20, c[i] = l[i] = NAN;

            // find the lines with the series of the specified country
            while (line = fgets(d, 65536, csv))
            {
               size_t len = strlen(line);
               size_t col = collen(line);
               if (len != col && (line = strstr(line+col+1, country)) && line[cl] == ',')
               {
                  // skip 3 more fields (Country/Region,Lat,Long)
                  for (i = 0; i < 3; i++)
                     line += collen(line) + 1;

                  // read the case numbers and sum it up into the c array, starting at day 1 = index 20+1;
                  int p, q = p = 21;
                  char *chk;
                  while (*line && *line != '\n' && *line != '\r')
                  {
                     ldouble v = strtod(line, &chk);
                     if (chk > line && isfinite(v))
                     {
                        if (isnan(c[q]))
                           c[q]  = v;
                        else
                           c[q] += v;
                        line = chk + 1;
                     }
                     else
                        line += collen(line) + 1;
                     q++;
                  }

                  if (q > p)
                  {
                     if (p > m) m = p;
                     if (q > n) n = q;
                  }
               }
            }

            if (n > m)
            {
               ldouble (*model)(ldouble t, Vector A) = LogisticFunctionGen;

               int j, k;
               ldouble ChiSqr,
                       A[6] = {2.0L*c[n-1], (model == LogisticFunctionGen) ? 0.2L : 5.0e-6L, 1.0L, 0.0L, 0.0L, 0.0L},
                      dA[6] = {};

               char *funcStr = (model == LogisticFunctionGen)
                             ? "a/(1 + exp(-b·(x - c)))"
                             : "a·0.5/(0.5 + (a - 0.5)·exp(-b·a·(x - c)))";

               if (isfinite(ChiSqr = curveFit(n - m, &t[m], &c[m], &l[m], k = 1, A, dA, model)))
               {
                  dA[0] = 0.0;
                  if (isfinite(ChiSqr = curveFit(n - m, &t[m], &c[m], &l[m], k = 3, A, dA, model)))
                  {
                     if (k <= 3)
                        fprintf(tsv, "# Model: %s\n", funcStr);
                     else if (k == 5)
                        fprintf(tsv, "# Model: %s\n", funcStr);
                     else if (k == 6)
                        fprintf(tsv, "# Model: %s\n", funcStr);

                     for (j = 0; j < k; j++)
                        fprintf(tsv, "#        %c = %9.6Lg ± %.5Lg %%\n", 'a'+j, A[j], dA[j]);
                     fprintf(tsv, "#   ChiSqr = %9.7Lg\n", ChiSqr);
                  }
                  else
                     fprintf(tsv, "# Curve fit failed\n");
               }

               // write the column header formular symbols and units.
               // - the formular symbol of time is 't', the unit symbol of day is 'd'
               // - the formular symbol of number of cases is C without a unit
               // - the formular symbol of the siumulated loogistic function is L without a unit
               fprintf(tsv, "t/d\tC\tL\n");
               for (i = 0; i < 366; i++)
                  if (isfinite(c[i]))
                     fprintf(tsv, "%.0Lf\t%.0Lf\t%.6Lf\n", t[i], c[i], l[i] = model(t[i], A));
                  else
                     fprintf(tsv, "%.0Lf\t*\t%.6Lf\n", t[i], l[i] = model(t[i], A));
            }
            else
               fprintf(tsv, "# No values for country %s encountered.\n", country);

            if (tsv != stdout)
               fclose(tsv);
         }

         if (csv != stdin)
            fclose(csv);
      }

   return 0;
}


// Curve fitting using the Levenberg–Marquardt least squares minimization algorithm
// https://en.wikipedia.org/wiki/Levenberg–Marquardt_algorithm
//
// GNU Octave's leasqr() function should give (almost) the same reseuls
// https://www.gnu.org/software/octave/


ldouble LUdecomposition(int m, Matrix A, Matrix LU, Index idx);
void LUbacksubstitution(int m, Matrix LU, Index idx, Vector B, Vector X);
void LUrefinment(int m, Matrix A, Matrix LU, Index idx, Vector B, Vector X);
void LUinversion(int m, Matrix A, Matrix LU, Index idx);


ldouble calcChiSqr(int n, Vector T, Vector C, Vector S,
                   int k, Vector A,
                   ldouble (*model)(ldouble t, Vector A))
{
   int     i, cnt = -1;
   ldouble chiSqr = 0.0L;

   for (i = 0; i < n; i++)
   {
      S[i] = model(T[i], A);
      chiSqr += sqr(S[i] - C[i]);
      cnt++;
   }

   if (cnt - k > 0)
      return chiSqr/(cnt - k);
   else if (cnt > 0)
      return chiSqr/cnt;
   else
      return NAN;
}


ldouble calcGradientCurvature(int n, Vector T, Vector  C,
                              int k, Vector A, Vector dA,
                              Vector beta, Matrix alpha, ldouble lambda,
                              ldouble (*model)(ldouble t, Vector A))
{
   int       i, p, q, cnt = 0;
   ldouble   ap, app, dYdA, yi;
   ldouble **V = malloc(k*sizeof(ldouble*));
   for (p = 0; p < k; p++)
      V[p] = calloc(n, sizeof(ldouble));

   for (p = 0; p < k; p++)
   {
      beta[p] = 0.0L;
      for (q = 0; q < k; q++)
         alpha[p][q] = 0.0L;

      ap = A[p];
      A[p] += dA[p];

      for (i = 0; i < n; i++)
         V[p][i] = model(T[i], A);

      A[p] = ap;
   }

   for (i = 0; i < n; i++)
   {
      cnt++;
      yi = model(T[i], A);
      for (p = 0; p < k; p++)
      {
         dYdA = (V[p][i] - yi)/dA[p];
         beta[p] += dYdA*(yi - C[i]);
         for (q = 0; q <= p; q++)
            alpha[p][q] += dYdA*(V[q][i] - yi)/dA[q];
      }
   }

   for (p = 1; p < k; p++)
      for (q = 0; q < p; q++)
         alpha[q][p] = alpha[p][q];

   for (p = 0; p < k; p++)
      for (q = 0; q < k; q++)
         alpha[p][q] /= beta[p];

   if (lambda == 0.0L)
   {
      for (p = 0; p < k; p++)
         if ((app = fabsl(alpha[p][p])) > lambda)
            lambda = app;
      lambda = sqrt(lambda)/cnt;
   }

   for (p = 0; p < k; p++)
   {
      alpha[p][p] *= (1.0L + lambda);
      free(V[p]);
   }
   free(V);
   return lambda;
}


ldouble ERelNorm(int m, Vector B, Vector R)
{
   int     i;
   ldouble sqSum = 0.0L;

   for (i = 0; i < m; i++)
      if (R[i] != 0.0L && isfinite(R[i]))
         sqSum += sqr(B[i]/R[i]);
      else
         sqSum += sqr(B[i]);
   return sqrt(sqSum);
}

ldouble curveFit(int n, Vector T, Vector C, Vector S,
                int k, Vector A, Vector dA,
                ldouble (*model)(ldouble t, Vector A))
{
   #define rtol 1.0e-9L
   #define atol 1.0e-12L

   static ldouble unitv[10] = {1.0L, 1.0L, 1.0L, 1.0L, 1.0L, 1.0L, 1.0L, 1.0L, 1.0L, 1.0L};

   bool      revoke;
   int       i, j, iter = 0;
   int       badIter = 0;
   int       idx[k];
   ldouble   ChiSqr, ChiSqr0, dChiSqr, dRelA = NAN, lambda0 = 0.0L, lambda = 0.0L;
   ldouble   beta[k], delta[k], B[k];
   ldouble **alpha, **alphaLU;

   alpha   = malloc(k*sizeof(ldouble *));
   alphaLU = malloc(k*sizeof(ldouble *));
   for (j = 0; j < k; j++)
   {
      dA[j]      = (A[j] != 0.0L) ? fabsl(A[j])*rtol : rtol;
      alpha[j]   = calloc(k, sizeof(ldouble));
      alphaLU[j] = calloc(k, sizeof(ldouble));
   }

   ChiSqr = calcChiSqr(n, T, C, S, k, A, model);
   do
   {
      iter++;
      lambda = calcGradientCurvature(n, T, C,
                                     k, A, dA, beta, alpha,
                                     lambda, model);
      if (!isfinite(dRelA = LUdecomposition(k, alpha, alphaLU, idx)))
         break;
      LUbacksubstitution(k, alphaLU, idx, unitv, delta);
      LUrefinment(k, alpha, alphaLU, idx, unitv, delta);

      dRelA = ERelNorm(k, delta, A);

      for (j = 0; j < k; j++)
         B[j] = A[j], A[j] -= delta[j];

      ChiSqr0 = ChiSqr;
      ChiSqr  = calcChiSqr(n, T, C, S, k, A, model);
      dChiSqr = ChiSqr - ChiSqr0;
      revoke  = dChiSqr >= 0;

      lambda0 = lambda;
      if (!revoke)
         lambda *= 0.2L;
      else
      {
         lambda = (lambda < 100.0L) ? lambda*3.0L : 0.0L;
         for (j = 0; j < k; j++)
            A[j] = B[j];
         ChiSqr = ChiSqr0;

         if (dChiSqr < atol)
            badIter++;
         else
            badIter = 0;
      }
   }
   while (iter < 1000 && badIter < 10 && isfinite(dRelA) && (dRelA > atol && fabsl(dChiSqr) > atol || dChiSqr > 0));

   if (isfinite(dRelA))
   {
      // extract lambda for preparation of the covarince matrix
      for (i = 0; i < k; i++)
      {
         for (j = 0; j < k; j++)
            alpha[i][j] *= beta[i];
         alpha[i][i] /= (1.0 + lambda0);
      }

      ldouble det = LUdecomposition(k, alpha, alphaLU, idx);
      if (det != 0.0L && isfinite(det))
      {
         LUinversion(k, alpha, alphaLU, idx);
         for (j = 0; j < k; j++)
            dA[j] = sqrtl(ChiSqr*fabsl(alpha[j][j]))*100.0L/fabsl(A[j]);
      }
   }
   else
      ChiSqr = NAN;

   for (j = 0; j < k; j++)
   {
      free(alpha[j]);
      free(alphaLU[j]);
   }
   free(alpha);
   free(alphaLU);

   #undef rtol
   #undef atol

   return ChiSqr;
}


// LU decomposition
// https://en.wikipedia.org/wiki/LU_decomposition

ldouble LUdecomposition(int m, Matrix A, Matrix LU, Index idx)
{
   int     i, j, k, imax = 0;
   ldouble max, sum, dum, d;
   Vector  V = calloc(m, sizeof(ldouble));

   if (LU != A)
      for (i = 0; i < m; i++)
         for (j = 0; j < m; j++)
            LU[i][j] = A[i][j];

   for (i = 0; i < m; i++)
   {
      max = 0.0L;
      for (j = 0; j < m; j++)
         if ((dum = fabsl(LU[i][j])) > max)
            max = dum;

      if (max != 0.0L)
         V[i] = 1.0L/max;
      else
      {
         free(V);
         return NAN;
      }
   }

   d = 1.0L;
   for (j = 0; j < m; j++)
   {
      for (i = 0; i < j; i++)
      {
         sum = LU[i][j];
         for (k = 0; k < i; k++)
            sum -= LU[i][k]*LU[k][j];
         LU[i][j] = sum;
      }

      max = 0.0L;
      for (; i < m; i++)
      {
         sum = LU[i][j];
         for (k = 0; k < j; k++)
            sum -= LU[i][k]*LU[k][j];
         LU[i][j] = sum;

         if ((dum = V[i]*fabsl(sum)) >= max)
         {
            max = dum;
            imax = i;
         }
      }

      if (j != imax)
      {
         for (k = 0; k < m; k++)
         {
            dum = LU[imax][k];
            LU[imax][k] = LU[j][k];
            LU[j][k] = dum;
         }
         V[imax] = V[j];
         d = -d;
      }
      idx[j] = imax;

      if (LU[j][j] == 0.0L)
         LU[j][j] = __LDBL_EPSILON__;

      if (j < m-1)
      {
         dum = 1.0L/LU[j][j];
         for (i = j+1; i < m; i++)
            LU[i][j] *= dum;
      }
   }

   free(V);
   return d;
}

void LUbacksubstitution(int m, Matrix LU, Index idx, Vector B, Vector X)
{
   int     i, j, k, l = -1;
   ldouble sum;

   if (X != B)
      for (i = 0; i < m; i++)
         X[i] = B[i];

   for (i = 0; i < m; i++)
   {
      k = idx[i];
      sum = X[k];
      X[k] = X[i];
      if (l >= 0)
         for (j = l; j <= i-1; j++)
            sum -= LU[i][j]*X[j];
      else if (sum != 0.0L)
         l = i;
      X[i] = sum;
   }

   for (i = m-1; i >= 0; i--)
   {
      sum = X[i];
      for (j = i+1; j < m; j++)
         sum -= LU[i][j]*X[j];
      X[i] = sum/LU[i][i];
   }
}

void LUrefinment(int m, Matrix A, Matrix LU, Index idx, Vector B, Vector X)
{
   int     i, j;
   ldouble sdp;
   Vector  R = calloc(m, sizeof(ldouble));

   for (i = 0; i < m; i++)
   {
      sdp = -B[i];
      for (j = 0; j < m; j++)
         sdp += A[i][j]*X[j];
      R[i] = sdp;
   }
   LUbacksubstitution(m, LU, idx, R, R);

   for (i = 0; i < m; i++)
      X[i] -= R[i];

   free(R);
}

void LUinversion(int m, Matrix A, Matrix LU, Index idx)
{
   int    i, j;
   Matrix Ai = malloc(m*sizeof(ldouble*));
   for (j = 0; j < m; j++)
      Ai[j] = calloc(m, sizeof(ldouble));

   for (i = 0; i < m; i++)
      Ai[i][i] = 1.0L;

   for (j = 0; j < m; j++)
      LUbacksubstitution(m, LU, idx, Ai[j], Ai[j]);

   for (i = 0; i < m; i++)
      for (j = 0; j < m; j++)
         A[i][j] = Ai[i][j];

   for (j = 0; j < m; j++)
      free(Ai[j]);
   free(Ai);
}
  1. Man öffne das Terminal und kompiliere unter FreeBSD oder macOS Linux obiges Programm mit dem folgenden Befehl:

        clang -g0 -O3 -march=native xcssecovid.c -Wno-parentheses -lm -o xcssecovid
  2. Man benutze fetch(1) oder curl(1), um die aktuelle Zeitreihen-CSV-Datei vom CSSE/JHU-GitHub-Repository herunterzuladen:

    curl -O https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv

    bzw.:

    fetch https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
  3. Man extrahiert den Datensatz für z.B. Deutschland (Germany) und exportiert die Zeitreihe, und den Verlauf der angepaßten logistischen Funktion in eine TSV-Datei:

        ./xcssecovid Germany time_series_covid19_confirmed_global.csv Germany.tsv

Stand 2020-03-22

Deutschland

# Model: a/(1 + exp(-b·(x - c)))
#        a =   41332.2 ± 4.5872 %
#        b =   0.32398 ± 2.6622 %
#        c =   59.5766 ± 0.48832 %
#   ChiSqr =  49811.09
t/d	C	L
-20	*	0.000000
-19	*	0.000000
-18	*	0.000001
-17	*	0.000001
-16	*	0.000001
-15	*	0.000001
-14	*	0.000002
-13	*	0.000003
-12	*	0.000004
-11	*	0.000005
-10	*	0.000007
-9	*	0.000009
-8	*	0.000013
-7	*	0.000018
-6	*	0.000025
-5	*	0.000034
-4	*	0.000047
-3	*	0.000065
-2	*	0.000090
-1	*	0.000124
0	*	0.000171
1	0	0.000237
2	0	0.000327
3	0	0.000453
4	0	0.000626
5	0	0.000865
6	1	0.001196
7	4	0.001654
8	4	0.002287
9	4	0.003162
10	5	0.004372
11	8	0.006045
12	10	0.008359
13	12	0.011557
14	12	0.015978
15	12	0.022092
16	12	0.030545
17	13	0.042232
18	13	0.058391
19	14	0.080732
20	14	0.111622
21	16	0.154330
22	16	0.213380
23	16	0.295022
24	16	0.407902
25	16	0.563972
26	16	0.779754
27	16	1.078096
28	16	1.490581
29	16	2.060878
30	16	2.849356
31	16	3.939471
32	16	5.446590
33	16	7.530182
34	16	10.410651
35	17	14.392582
36	27	19.896813
37	46	27.504655
38	48	38.018792
39	79	52.547028
40	130	72.617221
41	159	100.334536
42	196	138.595754
43	262	191.379648
44	482	264.137274
45	670	364.310540
46	799	502.009827
47	1040	690.877921
48	1176	949.150901
49	1457	1300.884357
50	1908	1777.225779
51	2078	2417.454787
52	3675	3269.257983
53	4585	4387.360258
54	5795	5829.295423
55	7272	7647.041481
56	9257	9873.970012
57	12327	12508.594135
58	15320	15499.833968
59	19848	18741.284670
60	22213	22081.172078
61	24873	25348.407929
62	*	28386.238438
63	*	31080.217599
64	*	33370.824992
65	*	35249.798128
66	*	36746.254413
67	*	37910.277774
68	*	38799.210377
69	*	39468.571690
70	*	39967.271905
71	*	40335.890401
72	*	40606.765036
73	*	40804.957415
74	*	40949.512917
75	*	41054.704735
76	*	41131.123695
77	*	41186.572370
78	*	41226.769702
79	*	41255.891950
80	*	41276.980700
81	*	41292.246891
82	*	41303.295429
83	*	41311.290133
84	*	41317.074350
85	*	41321.258880
86	*	41324.285931
87	*	41326.475567
88	*	41328.059399
89	*	41329.205004
90	*	41330.033620

Mit CVA (oder wenn es nicht anders geht eben mit Excel) läßt man sich den Kurvenverlauf anzeigen.

Solange man den Wendepunkt noch nicht sicher im Kasten hat, kann es sein, daß es noch zu groben Schwankungen bei den vorhergesagten Parametern der Kurvenanpassung kommt. So hatte z.B. der Ausreißer am vergangenen Freitag (20.3. = Tag 59) bewirkt, daß bei der Kurvenanpassung ein Grenzwert von ca. 140000 Infektionsfällen in Deutschland herauskam. Mit den zusätzlichen Daten von Samstag und Sonntag hat sich das dann wieder auf weniger als ⅓ reduziert. Man muß die Kurvenanpassungen also mit jeweils neuen Daten fortschreiben. Es hatte sich auch bei der Zeitreihenanalyse für China gezeigt, daß die Vorhersage mit jedem zusätzlichen Punkt genauer wird. Wenn die Daten nicht zu stark streuen, dann sollte man im Laufe dieser Woche für Deutschland zu einigermaßen validen Vorhersagen kommen können.

Italien

fetch -qo - https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv \
| ./xcssecovid Italy - Italy.tsv

# Model: a/(1 + exp(-b·(x - c)))
#        a =    119874 ± 5.2686 %
#        b =  0.195765 ± 2.0158 %
#        c =   61.1827 ± 0.80438 %
#   ChiSqr =  199699.2
t/d	C	L
-20	*	0.015017
-19	*	0.018265
-18	*	0.022214
-17	*	0.027018
-16	*	0.032860
-15	*	0.039966
-14	*	0.048608
-13	*	0.059119
-12	*	0.071903
-11	*	0.087452
-10	*	0.106362
-9	*	0.129362
-8	*	0.157335
-7	*	0.191358
-6	*	0.232737
-5	*	0.283064
-4	*	0.344274
-3	*	0.418721
-2	*	0.509265
-1	*	0.619388
0	*	0.753325
1	0	0.916224
2	0	1.114347
3	0	1.355313
4	0	1.648383
5	0	2.004826
6	0	2.438343
7	0	2.965601
8	0	3.606867
9	0	4.386792
10	2	5.335355
11	2	6.489017
12	2	7.892118
13	2	9.598584
14	2	11.673992
15	2	14.198093
16	2	17.267865
17	3	21.001237
18	3	25.541605
19	3	31.063324
20	3	37.778382
21	3	45.944498
22	3	55.874966
23	3	67.950594
24	3	82.634192
25	3	100.488140
26	3	122.195678
27	3	148.586681
28	3	180.668836
29	3	219.665315
30	3	267.060210
31	20	324.653263
32	62	394.625606
33	155	479.618546
34	229	582.827644
35	322	708.114582
36	453	860.139505
37	655	1044.516535
38	888	1267.995017
39	1128	1538.668490
40	1694	1866.212341
41	2036	2262.149140
42	2502	2740.137619
43	3089	3316.276480
44	3858	4009.407309
45	4636	4841.391173
46	5883	5837.320463
47	7375	7025.611119
48	9172	8437.900873
49	10149	10108.658581
50	12462	12074.391722
51	12462	14372.330601
52	17660	17038.478710
53	21157	20104.961756
54	24747	23596.696104
55	27980	27527.537811
56	31506	31896.259230
57	35713	36682.900685
58	41035	41846.201543
59	47021	47322.850462
60	53578	53029.138272
61	59138	58865.226360
62	*	64721.718317
63	*	70487.686377
64	*	76058.936901
65	*	81345.236876
66	*	86275.490172
67	*	90800.347149
68	*	94892.280475
69	*	98543.598845
70	*	101763.106164
71	*	104572.141774
72	*	107000.615502
73	*	109083.458570
74	*	110857.715215
75	*	112360.340989
76	*	113626.666687
77	*	114689.428656
78	*	115578.245094
79	*	116319.420779
80	*	116935.978106
81	*	117447.832556
82	*	117872.050870
83	*	118223.147912
84	*	118513.392415
85	*	118753.102714
86	*	118950.921402
87	*	119114.063348
88	*	119248.535118
89	*	119359.326169
90	*	119450.573489
91	*	119525.702090
92	*	119587.544051
93	*	119638.438795
94	*	119680.317183
95	*	119714.771783
96	*	119743.115400
97	*	119766.429727
98	*	119785.605691
99	*	119801.376869
100	*	119814.347127

Italien hat voraussichtlich noch einen schlimmen Monat vor sich.

Brasilien

fetch -qo - https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv \
| ./xcssecovid Brazil - Brazil.tsv

# Model: a/(1 + exp(-b·(x - c)))
#        a = 2.47301e+07 ± 2.0495e+05 %
#        b =  0.327684 ± 5.2709 %
#        c =   90.5531 ± 6976.6 %
#   ChiSqr =  385.0287
t/d	C	L
-20	*	0.000000
-19	*	0.000000
-18	*	0.000000
-17	*	0.000000
-16	*	0.000000
-15	*	0.000000
-14	*	0.000000
-13	*	0.000000
-12	*	0.000000
-11	*	0.000000
-10	*	0.000000
-9	*	0.000000
-8	*	0.000000
-7	*	0.000000
-6	*	0.000000
-5	*	0.000001
-4	*	0.000001
-3	*	0.000001
-2	*	0.000002
-1	*	0.000002
0	*	0.000003
1	0	0.000004
2	0	0.000006
3	0	0.000009
4	0	0.000012
5	0	0.000017
6	0	0.000023
7	0	0.000032
8	0	0.000044
9	0	0.000061
10	0	0.000085
11	0	0.000118
12	0	0.000164
13	0	0.000227
14	0	0.000315
15	0	0.000438
16	0	0.000607
17	0	0.000843
18	0	0.001170
19	0	0.001623
20	0	0.002253
21	0	0.003126
22	0	0.004339
23	0	0.006021
24	0	0.008356
25	0	0.011595
26	0	0.016091
27	0	0.022331
28	0	0.030990
29	0	0.043006
30	0	0.059682
31	0	0.082823
32	0	0.114938
33	0	0.159505
34	0	0.221353
35	0	0.307183
36	1	0.426294
37	1	0.591589
38	1	0.820977
39	2	1.139311
40	2	1.581079
41	2	2.194143
42	2	3.044921
43	4	4.225589
44	4	5.864061
45	13	8.137849
46	13	11.293299
47	20	15.672272
48	25	21.749190
49	31	30.182428
50	38	41.885643
51	52	58.126762
52	151	80.665338
53	151	111.943178
54	162	155.348869
55	200	215.584863
56	321	299.176876
57	372	415.180752
58	621	576.163324
59	793	799.563410
60	1021	1109.580186
61	1593	1539.793078
62	*	2136.796204
63	*	2965.239516
64	*	4114.819730
65	*	5709.972552
66	*	7923.304859
67	*	10994.200227
68	*	15254.571399
69	*	21164.468818
70	*	29361.247552
71	*	40727.320945
72	*	56483.269460
73	*	78315.297916
74	*	108548.761988
75	*	150382.594546
76	*	208202.521442
77	*	287992.898651
78	*	397865.563864
79	*	548714.902328
80	*	754983.615654
81	*	1035470.988032
82	*	1414017.615460
83	*	1919740.723057
84	*	2586272.862203
85	*	3449230.586999
86	*	4541081.192354
87	*	5883009.501097
88	*	7474671.897518
89	*	9284817.721708
90	*	11247588.711370
91	*	13268827.207360
92	*	15242645.959610
93	*	17072708.500855
94	*	18689654.910001
95	*	20058588.378917
96	*	21176273.231950
97	*	22062112.372848
98	*	22747811.715419
99	*	23268949.045561
100	*	23659526.693240

Für Brasilien ist es für eine sinnvolle Kurvenanpassung an die logistische Funktion noch zu früh, denn der Effekt der logistischen Begrenzung auf den exponentiellen Anstieg ist noch zu gering und geht in der Streuung der erhobenen Daten unter. Die Anpassung liefert für einen geraumen Zeitraum im Grunde die Kurve einer reinen Exponentialfunktion und ergibt dann zum Ende hin einen absurd hohen Grenzwert. Hier zeige ich den Verlauf in logarithmischer Darstellung, und zwar nur damit man sieht, wie eine fehlgeschlagene Kurvenanpassung aussieht. Mit diesen Werten arbeiten wir natürlich nicht. Auf der anderen Seite zeigt dieses Ergebnis im Vergleich zu den Kurven für Deutschland bzw. Italien, daß eine Kurvenanpassung die logistische Begrenzung herausarbeiten kann, wenn deren Effekt jenseits der Streuung durch die Datenerhebung bemerkbar wird.

Brasilien hat recht früh einschneidende Eindämmungsmaßnahmen ergriffen, praktisch zeitgleich wie Deutschland, aber auf dem Stand von wesentlich geringeren Fallzahlen. Die Chancen stehen also gut, daß der Infektionsverlauf weniger dramatisch ausfällt als anderswo.

Update 2020-03-23

Die neuen Anpassungen für Deutschland und Italien zeigen nur geringfügige Abweichungen zu denjenigen von gestern. Die Wendepunkte wurden offenbar überschritten. In DE vorgestern und in IT gestern. In Deutschland werden wir in ca. 2 Wochen den Grenzwert von vielleicht 45000 bestätigten Infektionsfällen erreichen. In Italien wird es wohl noch 3 Wochen dauern, und es wird wohl zu mehr als 115000 bestätigten Infektionsfällen kommen.

Für Brasilien gelang heute erstmals eine sinnvolle Kurvenpassung. Es ging hier später los, und der Wendepunkt wird voraussichtlich erst in etwa 1 Woche erreicht werden. In den Grenzwert von etwa 15000 ± 52 % bestätigten Infektionsfällen kämen wir demnach in etwa 3 Wochen. Die Unsicherheit der Anpassung ist allerdings noch recht groß, es können über 20000 Fälle oder auch unter 10000 Fälle werden. Sobald der Wendepunkt tatsächlich erreicht wurde, wissen wir mehr.

Deutschland

Parameter Wert Erläuterung
a 43857,6 ± 3,4 % Grenzwert der Infektionsfälle
b 0,316036 Steigung der e-Funktion
c 59,9318 Wendepunkt = 21.3.2020

Italien

Parameter Wert Erläuterung
a 115429 ± 3,8 % Grenzwert der Infektionsfälle
b 0,197855 Steigung der e-Funktion
c 60,8509 Wendepunkt = 22.3.2020

Brasilien

Parameter Wert Erläuterung
a 14416,8 ± 52 % Grenzwert der Infektionsfälle
b 0,327183 Steigung der e-Funktion
c 67,668 Wendepunkt = 28|29.3.2020

Update 2020-03-25

Deutschland

Die neuen Anpassungen für Deutschland befinden sich im Fluß. Der Grenzwert steigt von Mal zu Mal, und der Wendepunkt hat sich auch um einen Tag verschoben. Das ist den Streuungen bei den gemeldeten Fallzahlen zu verdanken.

Parameter Wert Erläuterung
a 52852,0 ± 2,7 % Grenzwert der Infektionsfälle
b 0,289108 Steigung der e-Funktion
c 61,1637 Wendepunkt = 22.3.2020

Italien

Die Anpassungen für Italien konvergieren, d.h. sie zeigen nur geringfügige Abweichungen zu den vorhergehenden. Der Wendepunkt und der vorhergesagte Grenzwert sind bei gleichzeitig abnehmender Standardabweichung stabil. Die halblogarithmische Auftragung zeigt, daß wir hier bereits das exponentielle Regime verlassen und uns schon sehr deutlich im logistischen Regime befinden.

Parameter Wert Erläuterung
a 112172 ± 2,1 % Grenzwert der Infektionsfälle
b 0,199858 Steigung der e-Funktion
c 60,5864 Wendepunkt = 22.3.2020

Brasilien

Die Anpassungen für Brasilien ergeben seit 2 Tagen einen deutlich rückläufigen Grenzwert und der Wendepunkt war demnach bereits vorgestern durchlaufen worden. Womöglich haben die einschneidenden Eindämmungsmaßnahmen, die Brasilien bereits vor über einer Woche bei der noch relativ niedrigen Fallzahl von 200 ergriffen hatte, eine entscheidende Wirkung gezeitigt. Wir sollten uns aber nicht zu früh freuen, d.h. wir warten die Ergebnisse der nächsten Tage ab, und zwar zu Hause.

Parameter Wert Erläuterung
a 3758,8 ± 4,3 % Grenzwert der Infektionsfälle
b 0,409505 Steigung der e-Funktion
c 62,0552 Wendepunkt = 23.3.2020

Update 2020-03-28

Deutschland

Die Rate, mit der neue Infektionsfälle gemeldet wurden, beschleunigte sich wieder im Laufe der jetzt vergangenen Woche, und zwar nach dem in der Vorwoche eigentlich schon ein Abklingen sichtbar wurde. Könnte es sein, daß in D von Woche zu Woche einfach mehr getestet wird, und deshalb einfach mehr Fälle erkannt werden, die in den Vorwochen in der Statistik nicht auftauchten? Die Logistische Funktion bildet sowas natürlich nicht vernünftig ab, und die Vorhersagen befinden sich im Fluß. So wie die Sache derzeit aussieht, könnten wir in D bei über 100k bestätigten Infektionsfällen landen.

Parameter Wert Erläuterung
a 104688 ± 6,5 % Grenzwert der Infektionsfälle
b 0,222636 Steigung der e-Funktion
c 66,2892 Wendepunkt = 27.3.2020

Italien

Die logistische Funktion läßt sich immer noch gut an den Verlauf der gemeldeten Infektionen in Italien anpassen. Aber auch hier steigt der Grenzwert von Mal zu Mal, wenn auch im Rahmen der Fehlerbandbreite der Erhebungen.

Parameter Wert Erläuterung
a 123623 ± 1,5 % Grenzwert der Infektionsfälle
b 0,190612 Steigung der e-Funktion
c 61,59 Wendepunkt = 22.3.2020

Brasilien

Die Anpassungen in Brasilien befinden sich nun auch im Fluß. Der weitere Verlauf muß verfolgt werden, und vielleicht gelingen ja dann in absehbarer Zukunft einigermaßen zuverlässige Aussagen.

Parameter Wert Erläuterung
a 5050,9 ± 2,8 % Grenzwert der Infektionsfälle
b 0,343718 Steigung der e-Funktion
c 63,7258 Wendepunkt = 25.3.2020

Diskussion auf Twitter: 1242279525955842048

Copyright © Dr. Rolf Jansen - 2020-03-22 23:45:39

PROMOTION