こんにちは、はっこんです。
もう寒いですね、ボクも上着を着るようになりました。
コロナウィルスの感染者数は日に日に増加していく一方ですね。
明日は我が身と思いながら、行動しましょう。
とりあえず11月の都内コロナウィルス陽性者数の推移はこんな感じです。
左：日別 右：累計
今回はこのデータを元に最小二乗法を用いて定量化してみようと思います。
最小二乗法とは？
最小二乗法とは大学に入ってから最初に習うものだったと思います。
データの組(x,y)が与えられた時に、それらの関係を表す関数を求める方法です。
y=ax+bへの近似
最も一般的な線形近似です。
簡単な説明なので気になるところがあれば自分で調べてください。
目的の直線を
y=ax+b
と仮定します。
あるx地点でのy座標はax+bなので、実際のy座標との差は、
ax+b-y
符号を正にするため二乗し、総和を取ります。
E=\sum_{i=1}^n(ax_i+b-y_i)^2これを最小にするa,bを見つければ良いわけです。
a,bの関数と考えて、
\frac{\partial E}{\partial a}=\sum_{i=1}^n2x_i(ax_i+b-y_i)=0\frac{\partial E}{\partial b}=\sum_{i=1}^n2(ax_i+b-y_i)=0bの方が簡単なので、まずはそちらから、
b=\frac{1}{n}\sum_{i=1}^n (y_i-ax_i)こいつを先ほどのaで偏微分した式に代入すると、
a=\frac{\sum_{i=1}^n x_iy_i - \frac{1}{n}\sum_{i=1}^nx_i\sum_{i=1}^ny_i}{\sum_{i=1}^nx_i^2-\frac{1}{n}(\sum_{i=1}^nx_i)^2}という感じですね、
y=bxaへの近似
目的の曲線を
y=bx^a
と仮定します。
両辺に対数を取ると、
\log y=\log(bx^a)
右辺を分解すると、
\log y=a\log(x)+\log(b)
対数を新たな変数として定義し直すと、
Y=aX+B
さきほどの一次関数と同じ形になりますね。
あとは先ほどと同様に定数を求めればOKです。
結果
y=ax+bへの近似
結果はこんな感じですね、
左：日別 右：累計
import numpy as np
import matplotlib.pyplot as plt
def leastsquares(x,y):
n = len(x)
a = (np.dot(x,y)-(sum(x)*sum(y))/n) / (sum(x**2)-(sum(x)**2)/n)
b = sum(y-a*x)/n
return a,b
x = np.arange(1,29)
y = np.array([116,87,209,122,269,242,294,189,157,293,316,392,374,352,255,180,298,485,533,522,539,391,314,186,401,481,570,561])
y_sum = np.array([31198,31285,31494,31616,31885,32127,32421,32610,32767,33060,33376,33768,34142,34494,34749,34929,35227,35712,36245,36767,37306,37697,38011,38197,38598,39079,39649,40210])
a,b = leastsquares(x, y)
a_sum,b_sum = leastsquares(x, y_sum)
plt.figure(1)
plt.clf()
plt.title('Number of coronavirus positive cases in Tokyo in November')
plt.scatter(x,y,color='k')
plt.xlabel('day')
plt.ylabel('people')
plt.plot([0,x.max()],[b,a*x.max()+b])
plt.text((x.min()+x.max())/6, (y.min()+y.max())/2,"y={:.3f}x+({:.3f})".format(a,b),size=10)
plt.figure(2)
plt.clf()
plt.title('Number of coronavirus positive cases in Tokyo in November')
plt.scatter(x,y_sum,color='k')
plt.xlabel('day')
plt.ylabel('people')
plt.plot([0,x.max()],[b_sum,a_sum*x.max()+b_sum])
plt.text((x.min()+x.max())/6, (y_sum.min()+y_sum.max())/2,"y={:.3f}x+({:.3f})".format(a_sum,b_sum),size=10)
plt.show()やはり日別の方は変動が激しいので微妙ですね、
類型の方は概ね一致していると思います。
y=bxaへの近似
こちらは累計のみ行いました。
データをそのまま用いた結果がこちら、
これはひどいですね、コードはこんな感じ、
import numpy as np
import matplotlib.pyplot as plt
def leastsquares(x,y):
n = len(x)
lnx = np.log(x)
lny = np.log(y)
a = (np.dot(lnx,lny)-(sum(lnx)*sum(lny))/n) / (sum(lnx**2)-(sum(lnx)**2)/n)
B = sum(lny-a*lnx)/n
b = np.exp(B)
return b*(X**a),a,b
x = np.arange(1,29)
y = np.array([31198,31285,31494,31616,31885,32127,32421,32610,32767,33060,33376,33768,34142,34494,34749,34929,35227,35712,36245,36767,37306,37697,38011,38197,38598,39079,39649,40210])
X = np.linspace(1,x.max(),500)
Y,a,b = leastsquares(x,y)
plt.figure(1)
plt.clf()
plt.scatter(x,y,color='k')
plt.plot(X,Y)
plt.xlabel('day')
plt.ylabel('people')
plt.text((x.min()+x.max())/4, (y.min()+y.max())/2,"y={:.3f}*x^{:.3f}".format(b,a),size=10)
plt.show()おそらく切片の値が大きすぎるためだと思います。
というわけで元データの人数から10月31日の陽性者数(31082)を引いてみました。
y = np.array([31198,31285,31494,31616,31885,32127,32421,32610,32767,33060,33376,33768,34142,34494,34749,34929,35227,35712,36245,36767,37306,37697,38011,38197,38598,39079,39649,40210])
y -= 31082一行加えただけです。
その結果がこちら、
これは綺麗ですね、結果的に11月1日以降は、
y=90.465x^{1.37}+31082と表せそうです。
ちなみに日別のデータを用いなかったのは、変動が激しいし、
ちょっとめんどいなと思ったからです。
コメント