%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

"""
Explanation: 多変量解析
圓川隆夫「多変量のデータ解析」(1988, 朝倉書店)をもとに多変量解析の演習を行います。
End of explanation
"""


data = pd.read_csv("tab22.csv") # http://shimotsu.web.fc2.com/Site/duo_bian_liang_jie_xi.html
data

"""
Explanation: 重回帰分析
以下のように変数を設定する。
$
y = 
\begin{bmatrix}
y_{1} \
y_{2} \
\vdots \
y_{n}
\end{bmatrix}
, 
X =
\begin{bmatrix}
x_{11} & \cdots & x_{1i} & \cdots & x_{1k}\
\vdots & \ddots & & & \vdots \
x_{i1} & & x_{ii} & & x_{ik} \
\vdots & & & \ddots & \vdots \
x_{n1} & \cdots & x_{ni} & \cdots & x_{nk}
\end{bmatrix}
,
\beta = 
\begin{bmatrix}
\beta_{1} \
\vdots \
\beta_{n}
\end{bmatrix}
,
\epsilon = 
\begin{bmatrix}
\epsilon_{1} \
\vdots \
\epsilon_{n}
\end{bmatrix}
$
このとき、最小自乗法により
$\hat \beta = (X'X)^{-1}X'y$となる。  
まず教科書p.28の表2.2を読み込む。
End of explanation
"""


y = np.asarray(data["y"])
X = np.c_[np.ones(len(data)), np.asarray(data[' x1'])]
beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
e = y - X.dot(beta)
R2 = 1 - e.dot(e.T)/sum((y - np.mean(y))**2)
print('y = ' + str(y))
print('X = \n' + str(X))
print('b = \n' + str(beta))
print('R2 = \n' + str(R2))

"""
Explanation: 最初に収益データ(y)に圧力(x1)を説明変数として単回帰モデルをあてはめる。
End of explanation
"""


plt.scatter(data[' x1'], data['y'])
plt.plot(data[' x1'], X.dot(beta))
plt.text(18, 31,"$R^{2}=%s$" % round(R2, 3), fontsize=17)

"""
Explanation: 故に、
$\hat y_{i} = 39.1 - 0.645 x_{1i}$
が得られ、寄与率は$R^{2} = 0.29$である。
End of explanation
"""


y = np.asarray(data["y"])
X = np.c_[np.ones(len(data)), np.asarray(data[' x2'])]
beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
e = y - X.dot(beta)
R2 = 1 - e.dot(e.T)/sum((y - np.mean(y))**2)
print('y = ' + str(y))
print('X = \n' + str(X))
print('b = \n' + str(beta))
print('R2 = \n' + str(R2))

"""
Explanation: 同様にして、温度(x2)を説明変数とした場合は、
End of explanation
"""


plt.scatter(data[' x2'], data['y'])
plt.plot(data[' x2'], X.dot(beta))
plt.text(91, 31,"$R^{2}=%s$" % round(R2, 3), fontsize=17)

"""
Explanation: 故に、
$\hat y_{i} = 11.8 +0.184x_{2i}$が得られ、寄与率は$R^{2} = 0.028$である。
End of explanation
"""


y = np.asarray(data["y"])
X = np.c_[np.ones(len(data)), np.asarray(data[[' x1', ' x2', ' x3']])]
beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
e = y - X.dot(beta)
R2 = 1 - e.dot(e.T)/sum((y - np.mean(y))**2)
print('y = ' + str(y))
print('X = \n' + str(X))
print('b = \n' + str(beta))
print('R2 = \n' + str(R2))

"""
Explanation: 収率データ(y)に圧力(x1)、温度(x2)、酸度(x3)の3個の説明変数で重回帰分析を行う。
End of explanation
"""


