import numpy as np
A = np.array([[1, 1], [1, 2]])
b = np.array([100, 180])
x = np.linalg.solve(A, b)
x
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
# NumPyを用いてデータを作成してプロットしてみる．
n = 5
x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([15.0, 20.0, 34.0, 42.0, 58.0])
plt.plot(x, y, ".")
import seaborn as sns
sns.jointplot(x=x, y=y, kind="reg")
print(np.dot(1 / np.dot(x, x) * x, y)) # NumPyの内積は np.dot()で計算できる．
print(1 / (x @ x) * x @ y) # 同じプログラム；Python 3.5からはベクトル（行列）の積は@ でできる．x@xはスカラであることに注意．
w = (x @ y) / (x @ x)
SSE = ((y - w * x) ** 2).sum()
SST = ((y - np.mean(y)) ** 2).sum()
R2 = 1.0 - SSE / SST
print(R2, SSE, SST)
import statsmodels.api as sm
model1 = sm.OLS(y, x)
r1 = model1.fit()
print(r1.summary2())
model2 = sm.GLM(y, x)
r2 = model2.fit()
print(r2.summary2())
import statsmodels.formula.api as smf
model3 = smf.glm("y ~ x - 1", {"y": y, "x": x})
r3 = model3.fit()
print(r3.summary2())
for i, r in enumerate([r1, r2, r3]):
plt.subplot(1, 3, i + 1)
plt.plot(x, y, ".")
if i < 2:
plt.plot(x, r.predict(x))
else:
plt.plot(x, r.predict({"x": x}))
plt.title(["OLS", "GLM", "glm"][i])
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0, x.max() + 1))
plt.ylim((0, y.max() + 5))
中心化
ベクトル $( y^{(1)},y^{(2)}, \ldots, y^{(n)} )$ の平均 $\bar{y}$ は， $$ \bar{y} = \left(\sum_i y^{(i)} \right)/n $$ と計算できる．
すべての要素が $1$ の長さ $n$ のベクトルとの内積を $1/n$ にすることによって計算できるので，ベクトル表記すると以下のように書ける． $$ \bar{y} = \frac{1}{n} (1,\ldots,1) \begin{pmatrix} y^{(1)}\\ \vdots\\ y^{(n)} \end{pmatrix} $$
もちろんNumPyのmean関数を使っても結果は同じである．
import numpy as np
x = np.array([1, 2, 3, 4, 5])
y = np.array([15, 20, 34, 42, 58])
ybar = (1 / 5) * np.ones(5) @ y
print(ybar)
print(np.mean(y))
平均が $0$ になるようにスケーリングすることを中心化(centering)とよぶ．
ベクトル $( y^{(1)},y^{(2)}, \ldots, y^{(n)} )$ の中心化 $\tilde{y}$ は， $$ \tilde{y} = y - \bar{y} \begin{pmatrix} 1 \\ \vdots\\ 1 \end{pmatrix} $$ と計算できる．
ytilde = y - ybar * np.ones(5)
print(ytilde)
$\bar{y}$ の定義を代入すると， $$ \tilde{y} = y - \bar{y} = \left( \begin{array}{ccc} 1&\cdots &0\\ &\ddots& \\ 0&\cdots &1 \end{array} \right) y - \frac{1}{n} \left( \begin{array}{ccc} 1&\cdots &1\\ &\ddots& \\ 1 &\cdots &1 \end{array} \right) y $$ となる．$n$行$n$列の単位行列を $I_n$，$n$行$n$列のすべての要素が $1$ の行列を $1_n$ とすると，中心化は $$ \left(I_n -\frac{1}{n} 1_n \right) y $$ と書くことができる．行列 $ \left(I_n -\frac{1}{n} 1_n \right)$ を $H_n$ と書くことにすると，中心化は $$ \tilde{y} = H_n y $$ という線形変換とみなすことができる．
H = np.identity(5) - (1 / 5) * np.ones((5, 5)) # H makes every column centered
print("H y = \n", H @ y)
行列 $H_n$ は正規化を行うための行列であるが，べき等性とよばれる重要な性質をもつ．
$$H_n H_n = H_n$$
これは2回正規化しても，1回した場合と同じ結果であることを表す．
print("H=\n", H)
print("H H=\n", H @ H)
行列 $X$ に対しても同様に，(列ごとの）平均をとるためには， $$ \bar{X} = \frac{1}{n} (1,\ldots,1) X $$ とすればよく，(列ごとに）中心化するには， $$ \tilde{X} = H_x X $$ とすればよい．
X = np.array([[1, 2, 3, 4, 5], [1, 0, 0, 0, 1]]).T
print("X=\n", X)
print("columnwise average=", (1 / 5) * np.ones(5) @ X)
print("H X =\n", H @ X)
分散
ベクトル $( y^{(1)},y^{(2)}, \ldots, y^{(n)} )$ の分散 $Var(y)$ は， $$ Var(y) = \frac{1}{n} \sum_i \left( y^{(i)} - \bar{y} \right)^2 $$ と定義される．これは $\frac{1}{n} \tilde{y}^T \tilde{y}$ であるので，$\tilde{y}= H_n y$ を代入すると， $$ Var(y) = \frac{1}{n} (H_n y)^T (H_n y) = \frac{1}{n} y^T H_n y $$ と書くことができる．ここで $H_n H_n = H_n$ という性質（べき等性）を利用している．
分散はNumPyのvar関数を用いても計算できる．
print("y^T H y =", y.T @ H @ y / 5)
print("variance=", np.var(y))
共分散
ベクトル $( x^{(1)},x^{(2)}, \ldots, x^{(n)} )$ と$( y^{(1)},y^{(2)}, \ldots, y^{(n)} )$ の共分散 $Cov(x,y)$ は， $$ Cov(y) = \frac{1}{n} \sum_i \left( x^{(i)} - \bar{x} \right) \left( y^{(i)} - \bar{y} \right) $$ と定義される．これは $\frac{1}{n} \tilde{x}^T \tilde{y}$ であるので，$\tilde{x}= H_n x$ と $\tilde{y}= H_n y$ を代入すると， $$ Cov(x,y) = \frac{1}{n} (H_n x)^T (H_n y) = \frac{1}{n} x^T H_n y $$ と書くことができる．ここで $H_n H_n = H_n$ という性質（べき等性）を利用している．
ベクトル $y$ にではなく，行列 $X$ に対して $$ \frac{1}{n} X^T H_n X $$ を計算すると分散・共分散行列を得ることができる．
print("X^T H X =\n", X.T @ H @ X / 5)
print(
"covariance= \n", np.cov(X, rowvar=False, bias=True)
) # rowbar=False：列(column)が変数, bias=True（Falseのとき不偏分散：n-1で割る）
原点を通過する1変数の回帰分析
平均が $0$ になるように中心化した後で，原点を通過する誤差最小の直線を求める． 直線の傾きを $w$ とすると，誤差は以下の関数になる． $$ f(w) = || \tilde{y} - w \tilde{x} ||^2 = \tilde{y}^T \tilde{y} - 2 w \tilde{y}^T \tilde{x} + w^2 \tilde{x}^T \tilde{x} $$ これを最小にする $w$ を求めればよい．上式を微分すると $$ f'(w)= 2 w \tilde{x}^T \tilde{x} - 2 \tilde{y}^T \tilde{x} $$ となるので，これを $0$ と置いて $w$ について解くと $$ w = \frac{\tilde{y}^T \tilde{x}}{ \tilde{x}^T \tilde{x} } $$ を得る．(2次微分の条件から関数$f(w)$が凸であることは各自確認されたい．）
$$ \tilde{y}^T \tilde{x} = x^T H_n y =n Cov(x,y) $$と $$ \tilde{x}^T \tilde{x} = x^T H_n x = n Var(x) $$ を用いると， $$ w =\frac{w \tilde{y}^T \tilde{x}}{ \tilde{x}^T \tilde{x} } = \frac{Cov(x,y)}{Var(x)} $$ を得る．
確率変数 $x,y$ の相関係数 $\rho_{xy}$ の定義は，$x,y$ の標準偏差を $\sigma_x,\sigma_y$ としたとき， $$ \rho_{xy} =\frac{Cov(x,y)}{ \sigma_x \sigma_y } $$ であったことを思い起こすと， $$ w = \frac{ \rho_{xy} \sigma_x \sigma_y }{\sigma_x^2 } = \rho_{xy} \frac{ \sigma_y }{ \sigma_x } $$ となる．
相関係数はNumPyのcorrcoef関数で得ることができる．
同じ結果が出るかどうか，確認してみよう．
xbar = (1.0 / n) * np.ones(n) @ x
xtilde = x - xbar * np.ones(n)
ybar = (1.0 / n) * np.ones(n) @ y
ytilde = y - ybar * np.ones(n)
print("正規化されたベクトルx,yに対するy^T x/x^T x=", (ytilde @ xtilde) / (xtilde @ xtilde))
print(
"正規化されたベクトルに対する 相関×yの標準偏差/xの標準偏差=",
np.corrcoef(xtilde, ytilde)[0, 1] * np.std(ytilde) / np.std(xtilde),
)
print("もとのベクトルに対する 相関×yの標準偏差/xの標準偏差=", np.corrcoef(x, y)[0, 1] * np.std(y) / np.std(x))
中心化した後で線形回帰を行って得た傾きは，上で求めた $w$ と一致するはずである．
確認してみよう．
model1 = sm.OLS(ytilde, xtilde)
r1 = model1.fit()
print(r1.summary2())
model = smf.glm("y ~ x -1", {"y": ytilde, "x": xtilde})
r = model.fit()
print(r.summary2())
model3 = smf.glm("y ~ x ", {"y": y, "x": x})
r3 = model3.fit()
print(r3.summary2())
$1.4$人の固定客がいることが推定されているが，$p$-値をみると大分大きい値であるので，固定客が $0$ であるという帰無仮説は棄却できない． つまり，友人は固定客になっていないということだ．
回帰分析の結果を相関と標準偏差から計算した結果と較べてみる．結果は同じになる．
w1 = np.corrcoef(xtilde, ytilde)[0, 1] * np.std(ytilde) / np.std(xtilde)
w0 = ybar - xbar * w1
print(w0, w1)
X = np.array([[1, 1, 1, 1, 1], [1, 2, 3, 4, 5]]).T
y = np.array([15, 20, 34, 42, 58])
XINV = np.linalg.inv(X.T @ X)
print((XINV @ X.T) @ y.T)
from sklearn.linear_model import LinearRegression
model = LinearRegression()
# 入力は縦ベクトル（行列）でなければならない．
# x=np.array([[1,2,3,4,5]]).T
x = np.array([1, 2, 3, 4, 5]).reshape(-1, 1)
y = np.array([[15, 20, 34, 42, 58]]).T
print(x)
model.fit(x, y)
print(model.coef_)
print(model.intercept_)
X = np.array([[1, 2, 3, 4, 5], [1, 0, 0, 0, 1]]).T # .Tで行列の転置を求めることができる．
y = np.array([15, 20, 34, 42, 58])
XINV = np.linalg.inv(X.T @ X)
(XINV @ X.T) @ y.T
model3 = smf.glm("y ~ x - 1", {"y": y, "x": x})
r3 = model3.fit()
print(r3.summary2())
print(r3.resid_pearson)
print(y - r3.predict())
XINV = np.linalg.inv(X.T @ X)
Hx = np.identity(5) - X @ XINV @ X.T
Hx @ y.T
from sklearn.linear_model import LinearRegression
model = LinearRegression(fit_intercept=False) # y-切片を0にする
X = np.array([[1, 2, 3, 4, 5], [1, 0, 0, 0, 1]]).T # .Tで行列の転置を求めることができる．
y = np.array([[15, 20, 34, 42, 58]]).T
model.fit(X, y)
print("係数=", model.coef_)
print("y-切片=", model.intercept_)
print("残差", model.predict(X))
import seaborn as sns
tips = sns.load_dataset("tips")
tips.head()
model = smf.glm("tip ~ total_bill + sex +smoker + day + size", tips)
# model = smf.glm('tip ~ total_bill + size ', tips)
# model = smf.glm('tip ~ total_bill', tips)
r = model.fit()
print(r.summary2())
sns.jointplot(x="tip", y="total_bill", data=tips, kind="reg")
tips["predict"] = r.predict()
tips["residual"] = tips["predict"] - tips["tip"]
tips.head()
sns.jointplot(x="total_bill", y="predict", data=tips, kind="reg")
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import OneHotEncoder
import pandas as pd
model = LinearRegression()
# カテゴリカルデータをダミー列を用いて数値に変換
tips = pd.get_dummies(
tips, columns=["sex", "smoker", "day", "time"], drop_first=True
) # drop_firstは男性，喫煙，木曜は従属なので入れないことを意味する．
X = tips[
[
"total_bill",
"size",
"sex_Female",
"smoker_No",
"day_Fri",
"day_Sat",
"day_Sun",
"time_Dinner",
]
]
y = tips["tip"]
model.fit(X, y)
print("係数=", model.coef_)
print("y-切片=", model.intercept_)
print("残差自乗和=", sum((y - model.predict(X)) ** 2))
import seaborn as sns
anscombe = sns.load_dataset("anscombe")
anscombe.head()
test = {}
test[0] = anscombe[anscombe.dataset == "I"]
test[1] = anscombe[anscombe.dataset == "II"]
test[2] = anscombe[anscombe.dataset == "III"]
test[3] = anscombe[anscombe.dataset == "IV"]
# grouped = anscombe.groupby('dataset') #とグループ化してもよい
for i in range(4):
print(test[i].describe())
for i in range(4):
# model1 = sm.OLS(test[i].y, test[i].x)
model1 = smf.glm("y ~ x ", {"y": test[i].y, "x": test[i].x})
r1 = model1.fit()
print(r1.summary2())
sns.lmplot(x="x", y="y", col="dataset", data=anscombe)
import seaborn as sns
titanic = sns.load_dataset("titanic")
titanic.head()
print(len(titanic))
mod = smf.glm(
formula="survived ~ sex +pclass ", data=titanic, family=sm.families.Binomial()
)
res = mod.fit()
print(res.summary2())
titanic["predict"] = res.predict()
titanic.head()
print("正解率", sum((titanic.predict > 0.5) == titanic.survived) / len(titanic))
sns.lmplot(x="fare", y="predict", col="sex", data=titanic)
titanic.dropna(subset=["age"], inplace=True)
mod = smf.glm(
formula="survived ~ pclass + sex + age", data=titanic, family=sm.families.Binomial()
)
res = mod.fit()
print(res.summary2())
titanic["predict"] = res.predict()
print("正解率", sum((titanic.predict > 0.5) == titanic.survived) / len(titanic))
一般化線形モデル（ポアソン回帰）
米国の州による死刑執行数
Number of Observations - 17
Number of Variables - 7
Variable name definitions::
EXECUTIONS - Executions in 1996
INCOME - Median per capita income in 1996 dollars
PERPOVERTY - Percent of the population classified as living in poverty
PERBLACK - Percent of black citizens in the population
VC100k96 - Rate of violent crimes per 100,00 residents for 1996
SOUTH - SOUTH == 1 indicates a state in the South
DEGREE - An esimate of the proportion of the state population with a
college degree of some kind
State names are included in the data file, though not returned by load.
data = sm.datasets.cpunish.load() # statsmodelの死刑執行数のデータを読み込み
df = pd.DataFrame(data.data)
df.head()
sns.pairplot(df)
r = smf.glm(
"EXECUTIONS ~ INCOME + PERPOVERTY + PERBLACK + VC100k96 + SOUTH + DEGREE",
df,
family=sm.families.Poisson(),
).fit()
df["predict"] = r.predict(df)
r.summary2()
sns.jointplot(x="EXECUTIONS", y="predict", data=df, kind="reg")
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.decomposition import PCA
%matplotlib inline
X = np.array([[1, 2, 3, 4, 5], [15, 20, 34, 42, 58]]).T
pca = PCA(n_components=1)
print(pca.fit_transform(X))
print(pca.get_covariance())
print(pca.explained_variance_ratio_)
s = pca.components_[0]
print(s)
plt.plot(X.T[0], X.T[1], "o")
plt.plot(X.T[0], np.mean(X.T[1]) + s[1] / s[0] * (X.T[0] - np.mean(X.T[0])), "r")
H = np.identity(5) - (1 / 5) * np.ones((5, 5)) # H makes every column centered
XC = H @ X
print("centering: \n H X =", XC)
print("variance= \n", np.cov(XC.T))
print(np.cov(X.T))
COV = (XC.T @ XC) / 4
print(COV)
lambda_, u = np.linalg.eig(COV) # 固有値分解
print(lambda_)
print(u)
u[np.argmax(lambda_)] # 最大固有値に対応する固有ベクトル
米国の州ごとの犯罪数データに対する主成分分析
Number of observations: 51 Number of variables: 8 Variable name definitions:
state
All 50 states plus DC.
violent
Rate of violent crimes / 100,000 population. Includes murder, forcible rape, robbery, and aggravated assault. Numbers for Illinois and Minnesota do not include forcible rapes. Footnote included with the American Statistical Abstract table reads: "The data collection methodology for the offense of forcible rape used by the Illinois and the Minnesota state Uniform Crime Reporting (UCR) Programs (with the exception of Rockford, Illinois, and Minneapolis and St. Paul, Minnesota) does not comply with national UCR guidelines. Consequently, their state figures for forcible rape and violent crime (of which forcible rape is a part) are not published in this table."
murder
Rate of murders / 100,000 population.
hs_grad
Precent of population having graduated from high school or higher.
poverty
% of individuals below the poverty line
white
Percent of population that is one race - white only. From 2009 American Community Survey
single
Calculated from 2009 1-year American Community Survey obtained obtained from Census. Variable is Male householder, no wife present, family household combined with Female household, no husband prsent, family household, divided by the total number of Family households.
urban
% of population in Urbanized Areas as of 2010 Census. Urbanized Areas are area of 50,000 or more people.
data = sm.datasets.statecrime.load()
df = pd.DataFrame(data.data)
df.head()
pca = PCA(n_components=2)
data = df.iloc[:, 1:]
data.head()
data -= data.mean()
data /= data.std()
pca.fit(data) # state 以外を主成分分析
print(pca.explained_variance_ratio_)
print(pca.components_) # 第1成分は暴力犯罪数，殺人数，貧困度を表し，第2成分は都市部の人口から殺人数，貧困度を引いた量を表す
x = pca.transform(data) # データを2次元に変換
# import networkx as nx
# pos ={}
# G=nx.Graph()
# for i, name in enumerate(df.state):
# pos[name] = (x[i,0],x[i,1])
# G.add_node(name)
# nx.draw(G,pos=pos, node_size=5, with_labels=True, font_size=8)
# #plt.savefig("statecrime.pdf")