Pythonで単回帰分析を実装する

今週はPythonで単回帰分析をやってみましょう。
モジュールを使えば簡単にできますが、それでは面白味がないので、回帰係数の計算までスクラッチで実装します。

回帰分析とは何か？

回帰分析とは２つの変数－仮にx,yとする－の間の関係を直線で考えることです。xに対するyの値を直線の関係で捉えることとも言えるでしょう。回帰分析は、xに対するyの値の予想に使うことができます。
　
例えば、アイスの売上と気温を考えてみましょう。
横軸に気温、縦軸にアイスの売上をプロットすると、以下のような関係がありました。

f:id:HirotoTakeda:20210203204641p:plain — 気温とアイスの売上の関係

0℃のときは1日に10個程度売れ、30℃のときには約４０個が売れています。
これは明らかに、気温が上がった時にアイスの売上も増加していると言えそうです。

ここで、これらの点に可能な限り寄り添うような直線を考えてみましょう。
直線というのはつまり、1次関数なので、こんな式を考えてみたいわけです。

$アイスの売上=気温\times傾き＋切片$

気温が高いほどアイスの売上が大きいことは自明です。しかし、実際にどの程度伸びているのでしょうか。つまり、１℃気温が上がったとき、何個売上は増えるのでしょうか？このようなことを回帰分析で明らかにします。

本当は複雑な数式展開があるのですが、なんだかんだ計算すると次の式で傾きと切片を求めることができます。

$傾き＝\frac{気温とアイスの売上を-緒に考えたときのばらつき(共分散)}{気温のばらつき(分散)}$

$切片＝アイスの売上の平均-傾き\times気温の平均$

これらを実際に求めて、同じ平面にプロットすると次のような図になります。

f:id:HirotoTakeda:20210203205854p:plain — 回帰直線

回帰式は $y = 0.967x+10.58$ でした。つまり、気温が１℃の上がると約10個売上が増えると言えます。
36℃のときの売上は、回帰式から

$y = 0.967×36 + 10.58 = 45.392 ≒45個$
と予測できるわけです。

Pythonで回帰式と図を出力する

最後にPythonで上の図と回帰式を得るためのコードを考えていきます。

モジュール

図をプロットするためのモジュールです

import matplotlib.pyplot as plt
import numpy as np

標準入力からxとyのデータを得ます。

x = list(map(int,input("xを入力").split()))
y = list(map(int,input("yを入力").split()))

#これは先程の気温とアイスの関係です。
#x = [0,5,10,15,20,25,30,35]
#y = [10,19,20,21,29,36,40,45]

平均

和を要素で割るだけですね

x_ave = sum(x)/len(x)
y_ave = sum(y)/len(y)

分散

分散は平均と各要素の差を２乗し、全て足し合わせていきます。
for文を使って下のように実装できます。

#xの分散
x_var = 0
for i in range(len(x)):
    x_var += (x_ave - x[i])**2
x_var = x_var/len(x)

共分散

共分散はxの各要素と平均の差とyの各要素を掛けたものの総和を取り、要素数で割ります。

#共分散の計算
cov = 0
for i in range(len(x)):
    cov += ((x[i]-x_ave)*(y[i]-y_ave))
cov = cov/len(x)

回帰係数

回帰係数は以下の通り計算すれば出ます。
$傾き＝\frac{x,yの共分散}{x分散}$
$切片＝y平均-傾き×x平均$

a = cov/x_var
b = y_ave - a*x_ave

#回帰式
s = "y = {0}x + {1}".format(a,b)
print("回帰式 : ",s)

プロット

直線と散布図をプロットします。
scatterは散布図、plt.plotで直線を書いていきます。
legendは凡例です。

Range = np.arange(0,max(x))
reg = a*Range + b
plt.title("Regression analysis")
plt.scatter(x,y)
plt.plot(Range,reg,color="r",label="Regression line")
plt.legend()
plt.xlabel("℃")
plt.ylabel("Ice")
plt.savefig("regression_result.png")

プログラム全体像

import matplotlib.pyplot as plt
import numpy as np
x = list(map(int,input("xを入力").split()))
y = list(map(int,input("yを入力").split()))
#x = [0,5,10,15,20,25,30,35]
#y = [10,19,20,21,29,36,40,45]
#x,yの平均
x_ave = sum(x)/len(x)
y_ave = sum(y)/len(y)
#xの分散
x_var = 0
for i in range(len(x)):
    x_var += (x_ave - x[i])**2
x_var = x_var/len(x)
#共分散の計算
cov = 0
for i in range(len(x)):
    cov += ((x[i]-x_ave)*(y[i]-y_ave))
cov = cov/len(x)
#回帰係数の計算
#傾き
a = cov/x_var
#切片
b = y_ave - a*x_ave
#回帰式
s = "y = {0}x + {1}".format(a,b)
print("回帰式 : ",s)
#直線と散布図のプロット
Range = np.arange(0,max(x))
reg = a*Range + b
plt.title("Regression analysis")
plt.scatter(x,y)
plt.plot(Range,reg,color="r",label="Regression line")
plt.legend()
plt.xlabel("℃")
plt.ylabel("Ice")
plt.savefig("regression_result.png")

ここまでを実行すると、先程も提示した

この図が得られます。

宣伝

ココナラでフリーのライターをやっています。
これまでに280件以上の受注実績があります。
大学入試の小論文、学生・社会人の方のレポートの作成・添削のご相談はこちらまで！
coconala.com