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

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

```math
\ln{G(t)}=\ln{G_{max}}-\frac{A_0}{k}e^{-kt}
```
と表せます。
等しい時間間隔の 3つの ln(累積発生数) を *A*, *B*, *C* とすると、

```math
\ln{G_{max}}=\frac{AC-B^2}{A-2B+C}
```
の関係が成り立ちます。([ゴンペルツ曲線とは何か？（3）](https://qiita.com/BMJr/items/4e2782e8f1d3c1db6801))
ln *G*(*t*) = *A* の時の時間 *t* = 0 とすると

```math
\frac{A_0}{k}=\ln{G_{max}}-A
```
等しい時間間隔を *s* とすると

```math
\ k = -\frac {1}{s}\ln\frac{C-B}{B-A}\\
```
![gomp3points.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193431/f081ae49-407e-9c02-7b0b-66dd497ca343.gif)
※上図は実測値を対数化したグラフです。
##不均等な時間間隔の3点を通る Gompertz 曲線
*A*, *B*, *C* が不均等な時間間隔（*t*<sub>0</sub>, *t*<sub>1</sub>, *t*<sub>2</sub>）の時は、どうすれば良いでしょうか。
![gomp3points2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193431/d7f3d2b8-4072-492f-c7bc-a0b8ef0285dc.gif)
等しい時間間隔の仮の値（*<sub>mid</sub>B*)を設定すると、*A*, *<sub>mid</sub>B*, *C* を通る Gompertz 曲線が得られます。得られた Gompertz 曲線の時間 *t*<sub>1</sub> における値を *<sub>calc</sub>B* として実測値 *B* と比較します。*<sub>mid</sub>B*を少しずつ増減して *<sub>calc</sub>B* と *B* が一致するまで繰り返します。
数学的に解を求めるのは難しいですが、計算機学的に for ループで微調整を繰り返せば解が得られるのです。
以下が、Excel のスクリーンショットと VBA のコードになります。
セル C4 が *t*<sub>0</sub>, C5 が *t*<sub>1</sub>, C6 が *t*<sub>2</sub>, D4 が *A*, D5 が *B*, D6 が *C* です。

| day             |G(t)            |
|:-----------------:|:-----------------:|
| 4  |      6000 |
| 22  |     68300|
| 43 |   80700|
上記データ（中国のCOVID-19累計感染者数）を入力してみましょう。
※マクロ内で *A*, *B*, *C* は対数に変換されますので、対数化していない元データを入力すれば良いです。
![excelSCRN.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193431/4a5a1b9e-917e-9e77-97b8-2215935bb2bd.gif)

```swift
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
```

[ゴンペルツ曲線とは何か？（1）](https://qiita.com/BMJr/items/24fca066926d471c651e)
[ゴンペルツ曲線とは何か？（2）](https://qiita.com/BMJr/items/1f9135edfb5e68b4f7da)
[ゴンペルツ曲線とは何か？（3）](https://qiita.com/BMJr/items/4e2782e8f1d3c1db6801)
[ゴンペルツ曲線とは何か？（4）](https://qiita.com/BMJr/items/f1665fcf61fc1e182311)
[ゴンペルツ曲線とは何か？（5）](https://qiita.com/BMJr/items/6e48af3017941589bd74)
[ゴンペルツ曲線とは何か？（6）](https://qiita.com/BMJr/items/60858ef7148732b04160)
[ゴンペルツ曲線とは何か？（8）](https://qiita.com/BMJr/items/acad4ec612086d20ac79)
[ゴンペルツ曲線とは何か？（9）](https://qiita.com/BMJr/items/775ae42bb5e0cc6dc461)
[ゴンペルツ曲線とは何か？（10）](https://qiita.com/BMJr/items/5a3fa10d72d107498604)
[ゴンペルツ曲線とは何か？（11）](https://qiita.com/BMJr/items/2055c27c0cb5596816df)
[ゴンペルツ曲線とは何か？（12）](https://qiita.com/BMJr/items/0cf2f2d2801ef52cda60)
[ゴンペルツ曲線とは何か？（13）](https://qiita.com/BMJr/items/608222a2e39cb01f2f52)
[ゴンペルツ曲線とは何か？（14）](https://qiita.com/BMJr/items/4757e087752500913f1b)
[COVID-19 数理モデル 〜 SIR と Gompertz と k値の関係〜](https://qiita.com/BMJr/items/83c122fa6623d80d4509)