🔰
数理モデルで慢性上咽頭炎の治療期間を予測する

2025/08/02に公開
初めてzennに投稿します。分析の記事を書くならnoteよりzennの方が良くないかと思い、今に至ります。

とはいえたくさんネタが出てくるわけでもないので、練習がてら以前noteに書いた記事をそのまま書き直しました。せっかくなのでちょっとだけグラフを追加しました。

 治療予測を行う背景僕は結構前から慢性上咽頭炎を患っている（た）。

概ね治療は終わったのだが、こいつは結構治療に時間がかかる。というのも、そもそも有効とされている上咽頭擦過療法：EATと呼ばれる治療を的確にできるお医者さんが少ないのである。
いわゆる最短治療に近い先生は各地に一人くらいしかいないため、通うのも手間で、遠隔地から治療に通うという人もいるらしい。

とあるコミュニティによると、100回EATをして半分程度治ったという声もあった。なかなか根気の必要な治療である。
上咽頭炎はこの画像の部位にあたり、治療の時は長い綿棒で強烈に擦られる。

ちなみに僕が初めて治療した日、痛くないだろと舐め過ぎた結果、治療後1時間待合室であまりの痛みに泣きまくった。
そんなわけで、闇雲に病院に通うよりも、ある程度見通しが立った状態で通院した方が希望が持てるのでは？と思い、得意の数理モデルで治療期間を予測してみることにした。
※念のためだが、医療知識は皆無なので、専門的な観点から予測を行ったものではない点はご留意いただきたい。

あくまで汎用的フレームワークに則った解釈である。

 方法患部がどのようにして治っていくのかを数理的にシミュレーションしていくには、一定程度仮定を置く必要がある。（「1回の治療では、n%の治療効果があるとする」など。）

 患部の表現患部は nマス×nマス のグリッドで表現
各グリッド (i,j) の回復度 R_{i,j} は 0〜1の間の値
初期状態はすべて0（未回復）

 薬の塗布確率1回の治療でグリッド上の1セグメントに薬の塗布ができる
1セグメントの選択は確率的に選択され、この確率分布は多変量正規分布に従うものとする
中央が最も薬が当たりやすいとする

 治療効果1回の治療で選択されたセグメントの回復度は治療回数とともに増加するが、効果は小さくなる
各セグメントについて、具体的に1回の治療の回復量は、以下の数式で表現する。治療のたびにα分回復するとする。2回目の治療では、残りの部分に効果があるとする。これを繰り返すと効果は低減していく

ΔR_{i,j} = α \cdot (1 - R_{i,j})

 全体の回復度全体に占める治療された部位の割合である。

R_{total} = \frac{1}{n \times n} \sum_{i=1}^{n} \sum_{j=1}^{n} R_{i,j}



 シミュレーション
 手順n×n 個のグリッドのうち 1つ選択
選択したセグメントが α 治る
全体の治癒度を計算
これを 500回繰り返す（治療500回）
1〜4のシミュレーションを 30回実施
ちなみに、一部の人から聞いた「治療回数100回で50％程度回復した」という話を参考にして、1回の回復度 α=0.3 と仮定してシミュレーションしています。

 結果シミュレーション結果は以下の通り。

現時点で100回程度治療した方は、さらに100回通院すると 60〜70％程度まで回復する見込みがあることがわかる。



※「治療回数100回で50％程度回復」の目安ラインもグラフに表示
次に、100回治療した時点でどれくらいの全体回復度になることが多いかを把握するために、上述のシミュレーショングラフを100回時点で切り取って回復度をヒストグラムにしてみた。

すると大体52~54%くらいの回復度であることがわかった。


 示唆今回の数理シミュレーションから、ある程度現実的に通院回数を推定することができそうだった。

さらに、回復の速度を上げたい場合は、以下のような制御が現実的なレバーになるとわかった。
1回の治療率 α を上げる
1回で治療できる範囲（セグメント）を広くする



 コードPython によるシミュレーションコードはこちらです。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
# 日本語のmatplotlibプロットをサポート
import japanize_matplotlib 

# グリッドサイズ
GRID_SIZE = 5
NUM_CELLS = GRID_SIZE * GRID_SIZE

# 初期回復度 (0: 未回復, 1: 完全回復)
recovery_grid = np.zeros((GRID_SIZE, GRID_SIZE))

# ガウス分布のパラメータ設定（中心を (2,2) にして、広がりを設定）
mean = [2, 2]  # 中心
cov = [[1, 0], [0, 1]]  # 分散共分散行列（中心付近が高確率）
x, y = np.meshgrid(np.arange(GRID_SIZE), np.arange(GRID_SIZE))
pos = np.dstack((x, y))
probability_distribution = multivariate_normal(mean, cov).pdf(pos)
probability_distribution /= probability_distribution.sum()  # 正規化

# 治療パラメータ
# 試行回数
num_trials = 30
alpha = 0.3  # 固定値
num_treatments = 500  # 治療回数

plt.figure(figsize=(8, 5))

for trial in range(num_trials):
    # 初期回復度リセット
    recovery_grid = np.zeros((GRID_SIZE, GRID_SIZE))
    recovery_rates = []

    for t in range(num_treatments):
        # どのセグメントに治療を行うかを確率的に選択
        indices = np.random.choice(NUM_CELLS, p=probability_distribution.ravel())
        print(indices)
        i, j = divmod(indices, GRID_SIZE)

        # 治療による回復度の増加（逓減効果）
        recovery_grid[i, j] += alpha * (1 - recovery_grid[i, j])
        recovery_grid[i, j] = min(recovery_grid[i, j], 1)  # 1を超えないように制限

        # 全体の回復率を計算
        recovery_rate = recovery_grid.mean()
        recovery_rates.append(recovery_rate)

    # プロット（各試行ごとに異なる色で描画）
    plt.plot(recovery_rates, 
            #  label=f"Trial {trial+1}"
             )

# 50%回復ライン
plt.axhline(0.5, color='r', linestyle='--', label="50% Recovery")

plt.xlabel("治療回数")
plt.ylabel("全体回復度")
plt.title(f"治療回数と回復度 α={alpha}")
plt.legend()
plt.grid()
plt.show()