More than 3 years have passed since last update.

ゴンペルツ曲線とは何か？（7）

Last updated at 2021-06-20Posted at 2020-03-22

#Excel 版 Gompertz 曲線計算機
ゴンペルツ曲線とは何か？（6）で新型コロナウイルス（COVID-19）の累積患者数を任意の3点の値から数式化する方法を示しました。
本項では Excel VBA のプログラムを紹介します。
##等しい時間間隔の3点を通る Gompertz 曲線
時間を t、累積発生数を G(t)、対数化した値を ln G(t) とすると

\ln{G(t)}=\ln{G_{max}}-\frac{A_0}{k}e^{-kt}

と表せます。
等しい時間間隔の 3つの ln(累積発生数) を A, B, C とすると、

\ln{G_{max}}=\frac{AC-B^2}{A-2B+C}

の関係が成り立ちます。(ゴンペルツ曲線とは何か？（3）)
ln G(t) = A の時の時間 t = 0 とすると

\frac{A_0}{k}=\ln{G_{max}}-A

等しい時間間隔を s とすると

\ k = -\frac {1}{s}\ln\frac{C-B}{B-A}\\

※上図は実測値を対数化したグラフです。
##不均等な時間間隔の3点を通る Gompertz 曲線
A, B, C が不均等な時間間隔（t₀, t₁, t₂）の時は、どうすれば良いでしょうか。

等しい時間間隔の仮の値（_midB)を設定すると、A, _midB, C を通る Gompertz 曲線が得られます。得られた Gompertz 曲線の時間 t₁ における値を _calcB として実測値 B と比較します。_midBを少しずつ増減して _calcB と B が一致するまで繰り返します。
数学的に解を求めるのは難しいですが、計算機学的に for ループで微調整を繰り返せば解が得られるのです。
以下が、Excel のスクリーンショットと VBA のコードになります。
セル C4 が t₀, C5 が t₁, C6 が t₂, D4 が A, D5 が B, D6 が C です。

day	G(t)
4	6000
22	68300
43	80700
上記データ（中国のCOVID-19累計感染者数）を入力してみましょう。
※マクロ内で A, B, C は対数に変換されますので、対数化していない元データを入力すれば良いです。

Sub Gompertz()
If ActiveSheet.ChartObjects.Count = 1 Then ActiveSheet.ChartObjects(1).Delete
   t0 = Range("C4").Value
   t1 = Range("C5").Value
   t2 = Range("C6").Value
   actualA = Range("D4").Value
   actualB = Range("D5").Value
   actualC = Range("D6").Value
    
   A = Log(actualA)
   B = Log(actualB)
   C = Log(actualC)

   Dim midBupper As Double
   Dim midBlower As Double
   Dim midB As Double

   midBupper = (A + C) / 2
   midBlower = C
   midB = (midBupper + midBlower) / 2

    If t0 >= t1 Or t1 >= t2 Then
        MsgBox "t0 < t1 <t2としてください。": Exit Sub
    End If
    
    abslope = (B - A) / (t1 - t0)
    bcslope = (C - B) / (t2 - t1)
    acslope = (C - A) / (t2 - t0)

    If abslope * bcslope <= 0 Then
        MsgBox "正->負や負->正 の傾きにならないようにして下さい。": Exit Sub
    End If
    
    If acslope > 0 And abslope > bcslope Then
        midBupper = C
        midBlower = (A + C) / 2
    End If
     
    If acslope > 0 And abslope < bcslope Then
        midBupper = (A + C) / 2
        midBlower = A
    End If
     
    If acslope < 0 And abslope < bcslope Then
        midBupper = (A + C) / 2
        midBlower = C
    End If
     
    If acslope < 0 And abslope > bcslope Then
        midBupper = A
        midBlower = (A + C) / 2
    End If
     
      midB = (midBupper + midBlower) / 2

    For i = 0 To 100
           lnGmax = (A * C - midB ^ 2) / (A - 2 * midB + C)
           k = -1 / ((t2 - t0) / 2) * Log((midB - C) / (A - midB))
           calcB = lnGmax - (lnGmax - A) * Exp(-k * (t1 - t0))
           If calcB = B Then Exit For
           If B > calcB Then midBlower = midB
           If B < calcB Then midBupper = midB
           midB = (midBupper + midBlower) / 2
    Next i
    
            Range("C11").Value = "lnG(t) = lnGmax -  (A0/k)*EXP(-k*(t-t0)) "
            Range("D12").Value = k
            Range("D13").Value = Exp(lnGmax)
            Range("C12").Value = "k"
            Range("C13").Value = "Gmax"
            Range("C14").Value = "lnGmax"
            Range("C15").Value = "A0/k"
            Range("C16").Value = " t0"
            Range("D14").Value = lnGmax
            Range("D15").Value = lnGmax - A
            Range("D16").Value = t0
            
            Cells(18, 3) = " time"
            Cells(18, 4) = "G(t)"
            Cells(18, 5) = "lnG(t)"
            
          For i = t0 To t2
             Cells(i - t0 + 19, 3) = i
             Cells(i - t0 + 19, 4) = Exp(lnGmax - (lnGmax - A) * Exp(-k * (i - t0)))
             Cells(i - t0 + 19, 5) = lnGmax - (lnGmax - A) * Exp(-k * (i - t0))
          Next i
             
     Set trange = Worksheets("Sheet1").Range(Cells(19, 3), Cells(t2 - t0 + 19, 4))
      Charts.Add
      With ActiveChart
            .HasLegend = False
            .ChartType = xlXYScatterLines
            .SetSourceData Source:=trange, PlotBy:=xlColumns
            .Location Where:=xlLocationAsObject, Name:="Sheet1"
      End With
      With ActiveChart
            .Axes(xlValue, xlPrimary).HasTitle = True
            .Axes(xlValue, xlPrimary).AxisTitle.Characters.Text = "G(t)"
            .Axes(xlValue, xlPrimary).AxisTitle.AutoScaleFont = True
            .Axes(xlCategory, xlPrimary).HasTitle = True
            .Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = "time"
            .Axes(xlCategory, xlPrimary).AxisTitle.AutoScaleFont = True
            .PlotArea.Interior.ColorIndex = 2
            .Axes(xlValue).MajorGridlines.Delete
      End With
      With ActiveSheet
            .ChartObjects(1).Top = 50
            .ChartObjects(1).Width = 350
            .ChartObjects(1).Height = 350
            .ChartObjects(1).Left = 600
      End With
             
End Sub

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up