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