Pythonで機械学習をしてみた話
この記事は群馬高専 Advent Calendar 2022 20日目の記事です
ちなみに僕の誕生日です、祝ってください←←←
ってことで今回は機械学習の話をします
数学の話やらなんやら色々出てきますがまあなんかてきとうに見てってください()
機械学習
機械学習とは
まずは機械学習って何?って話ですけど、言葉の定義とかはよく知らないんでそれは置いといて←
ここでは機械学習の中の2つの分野について話していきます
1つ目は回帰
これは連続するデータから値を予測することです、つまりグラフのプロットから曲線を予測するみたいなものですね
実験レポートでグラフとかよく書きますよね、その線を引く作業を機械がやる感じです
2つ目は分類
これは言葉のまんまで、データからそのデータの種類を分けるやつです
ただし今回は2値分類(分類先が2つの分類)のみを扱います
って感じの話をします、あと実装はPythonでやっていきます
Pythonはnumpyとかあるので今回みたいなのには便利ですよ
回帰
回帰の理論
回帰です
こいつは簡単に言うと微分を使います
唐突に大事なところを言っても伝わらないのでちゃんと話をします()
まずは数式とかなしで概要的な話を
といっても、回帰っていうのはデータがあってそのデータに沿ったグラフになる数式を見つけようぜという話です
最初にてきとうなパラメータで数式を決めます、これはまあほんとになんも考えなくていいです
そしてその式のデータと実在するデータとで誤差を求めます
ここまで来て、微分を使いたいので微分の話をちょっとします
1年生だと微分はまだやってないと思うんですけど(1年生がこの記事を見てるとは思ってないですけど(というかこの記事読む人なんていrrrr))、簡単に言うとあるグラフ上での各点の傾きを求める関数を求めることです、求めた関数は導関数とも言います
微分の応用として、そのグラフの減少と増加を調べるっていうのがありましたよね、それを使います
(傾きはその点でのyの増加率ともいえるわけです、増加率が正ならxを増やしたときyは増え、負ならxを増やしたときyは減るはずですね)
さて、誤差の関数を最初の式のパラメータで微分してみたらどうでしょう?このときの増加か減少で小さい方へ行くようにパラメータを変えたら、誤差が最小値(正確には極小値)を取るようなパラメータが求められると思いませんか!?求められるんですよぉぉぉ！！！！
…というのが回帰です
この説明に沿って数式を紹介していきます
まずは実在データに沿わせたい関数を決めます
$$
f_\theta(x)=\theta_1x+\theta_0
$$
このθはパラメータです、1次関数になってますがパラメータを増やしてn次に拡張可能です(これを多項式回帰といいます)
また、実在データが高次元(予測したいデータがy、それを決める要素となるデータがx_1,x_2,…,x_m)の場合についても、それぞれのx_i(1<=i<=m)についてのn次関数の和と考えます(これを重回帰といいます)
もちろん、多項式回帰と重回帰を同時に行うこともできます
以下では、f_θ(x)は多項式回帰として扱います
次は回帰の目的関数(誤差の関数)、これは差の2乗について総和を取ります(差は正だったり負だったりするので2乗取って正にします、離れてるかどうかが見たいので)
$$
E(\theta)=\frac12\sum_{k=1}^n(y_k-f_\theta(x_k))^2\\
ただしx_k,y_kはk番目のデータ、nはデータの個数
$$
半分にしてるのは微分したときに楽だからです((x^2)'=2xですからね)
誤差の極小値となるパラメータを求める分にはこの定数は影響しません
このような誤差の関数を使った回帰を最小2乗法とか言います
これを微分します(正確にはθ_iについての偏微分ですが、偏微分が分からない場合は微分する変数以外を定数として扱うものとして思ってもいいです)
$$
\frac\partial{\partial \theta_i}E(\theta)\\
=\frac\partial{\partial \theta_i}\frac12\sum_{k=1}^n(y_k-f_\theta(x_k))^2\\
=\frac12\frac\partial{\partial \theta_i}\sum_{k=1}^n(y_k^2-2y_kf_\theta(x_k)+f_\theta(x_k)^2)\\
=\frac12\sum_{k=1}^n(2f(\theta_i)\frac\partial{\partial \theta_i}f_\theta(x_k)-2y_k\frac\partial{\partial \theta_i}f_\theta(x_k))\\
\frac\partial{\partial \theta_i}f_\theta(x_k)はx_k^iだけが残るので\\
=\sum_{k=1}^n(f_\theta(x_k)-y_k)x_k^i
$$
微分ができたので、これを使ってパラメータの更新式を定めます
$$
\theta_i:=\theta_i-\eta\sum_{k=1}^n(f_\theta(x_k)-y_k)x_k^i
$$
ただし、E(θ)の微分につけ加えたηは学習率といいます(定数です、あとで説明します)
これは、古いθから新しいθを定義するという意味です、:=はプログラミングでいう=だと思っていいです
誤差が単調増加ならθを小さくしたい、また単調増加のとき誤差の微分は正であるので、これを引く、単調減少の場合についても同様に引く、というような考えです
この更新式を実行して、値が収束するまで繰り返します(実際のプログラムではある程度収束したら打ち切ります)
学習率をかける意図としては、θの値が急激に変化することを防ぐためです(精度を上げるためと、あと場合によっては収束しなくなってしまうので)
なので基本的にηは1未満の数で設定します
ここまでの回帰は勾配降下法とか最急降下法といいます、目的は誤差の最小値を求めることですが、極小値に捕らわれてしまうという問題があります(これを局所解といいます、一応極小値の説明をすると、yが周辺より小さくなるx、すなわちグラフが谷になるとこの頂点のx座標のことです)
そこで、誤差の求め方を変えて、毎回ランダムなデータを選んでそれについてのみの誤差を求めるようにします(つまり、E(θ)において総和を取らず、kをランダムな値とします)
これによって、1回の更新で1つのデータについての誤差を小さくすることになります、このデータの選択をランダムに変えていけば全体についての誤差も小さくなります
そして大事なのは、こうすると局所解に捕まりにくくなります、極小値を取るθを超えて値が変化することがあるからです
これを確率的勾配降下法といいます、出てくる式はほぼ同じなので更新式だけ示します
$$
\theta_i:=\theta_i-\eta(f_\theta(x_k)-y_k)x_k^i\\
ただし、kはランダムな値
$$
回帰の説明はこんな感じです、ってことでこれをPythonで確認してみます
回帰の実験方法
Pythonで実際のコードをだす前にちょっとやりたいことがあって
まずz-score正規化っていうのをしたいんですよ
原理はよく知らないので詳しい話はしませんが
こんな感じのzを定義して、xをこのzと置き換えて回帰をします
$$
z=\frac{x-\mu}{\sigma}\\
ただし、\muは平均値、\sigmaは標準偏差
$$
要するにx軸の縮尺を収束しやすいように変更します
収束しやすい理由は知りません()
それと、回帰後の関数の誤差を以下の方法で求めて表示します
$$
\text{MSE}=\frac{1}{n}\sum_{i=1}^{n}\left(y_i-f_\theta(x_i)\right)^2\\
\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y_i-f_\theta(x_i)\right)^2} \\
\text{MAE}=\frac{1}{n}\sum_{i=1}^n\left|y_i-f_\theta(x_i)\right| \\
\text{MAPE}=\frac{1}{n}\sum_{i=1}^n\left|\frac{f_\theta(x_i)-y_i}{y_i}\right|\\
$$
それぞれ平均2乗誤差(MSE)、2乗平均平方根誤差(RMSE)、平均絶対値誤差(MAE)、平均絶対値誤差率(MAPE)といいます
データはテスト用データと学習用データに分けます、回帰を行うのは学習用データで、誤差を求めるのはテスト用データです
大体データ数の比率は学習:テスト=7:3とか8:2ぐらいですかね
今回はデータにノイズ付きで生成した4次関数のデータを使います(これもPythonで作りました())
Pythonで回帰
ということで実際にプログラムを書いていきます
全部は説明しませんが()
まず、モジュールはnumpyとmatplotlibを使います(あとrandom)
numpyは行列が使えるからです、for文書かなくても一度に計算できます
(numpy.arrayについて、スカラーと演算すれば各要素にその演算が行われ、array同士で演算すれば対応するインデックス同士で演算が行われます、行列の積はnumpy.dotを使います)
matplotlibはデータを渡すとグラフを表示してくれるやつです
ちなみに使う回帰は多項式回帰の最急降下法です
生成したデータはcsvファイルに突っ込んでおきます
↓こんな感じで
x,y
2,17
…
ということで、いきなりコード全部貼ります
コメント書きまくってるので許して()
import numpy as np
import matplotlib.pyplot as plt
from random import randint
train=np.loadtxt("data.csv",delimiter=",",skiprows=1) # データの読み込み
train_x,train_y=train[:,0],train[:,1] # x,yを取り出す
test_x,test_y=np.array([]),np.array([]) # テストデータ用
for i in range(int(len(train_x)*3/10)): # データの3割をテスト用データとする
j=randint(0,train_x.shape[0]-1) # ランダムなインデックス
test_x=np.append(test_x,train_x[j])
test_y=np.append(test_y,train_y[j])
train_x=np.delete(train_x,j)
train_y=np.delete(train_y,j)
mu=train_x.mean() # 平均
sigma=train_x.std() # 標準偏差
standardize=lambda x:(x-mu)/sigma # z-score正規化を行う関数
train_z=standardize(train_x) # z-score正規化を通したx
print("最急降下法\nデータはdata.csvより(z-score正規化を通す)\n次数を入力:")
jisu=int(input()) # 次数を受け取る
t=np.array([1]*(jisu+1)) # パラメータ
f=lambda x,theta:np.dot(np.array([x**i for i in range(jisu+1)]).T,theta) # 曲線f(x)
E=lambda theta:0.5*np.sum((train_y-f(train_z,theta))**2) # 誤差E(Θ)
diff=1 # 前回の更新との誤差の差
error=E(t) # 誤差
print("学習率を入力:")
ETA=float(input()) # 学習率η
print("青:学習用データ,水:テスト用データ,緑:得られた関数")
while diff>0.0001: # パラメータΘの更新
t_f=t
t=t-ETA*(np.dot(f(train_z,t)-train_y,np.array([train_z**i for i in range(jisu+1)]).T))
diff=E(t_f)-E(t)
x=np.linspace(3,-3,100)
plt.plot(train_z,train_y,'o',color="b") # 学習データのプロット
plt.plot(standardize(test_x),test_y,'o',color="c") # テスト用データのプロット
plt.plot(x,f(x,t),color="g") # 回帰で得た曲線
MSE=lambda x,y,theta: np.sum((y-f(x,theta))**2)/x.shape[0] # MSEを計算する関数
MAE=lambda x,y,theta:np.sum(np.abs(y-f(x,theta)))/x.shape[0] # MAEを計算する関数
MAPE=lambda x,y,theta:np.sum(np.abs((f(x,theta)-y)/y))*100/x.shape[0]# MAPEを計算する関数
print(f"パラメータΘ={t.round(3).tolist()}")
print(f"平均2乗誤差(MSE)={round(MSE(standardize(test_x),test_y,t),2)}")
print(f"2乗平均平方根誤差(RMSE)={round(MSE(standardize(test_x),test_y[1],t)**0.5,2)}") # RMSEはMSEを0.5乗(平方根)
print(f"平均絶対値誤差(MAE)={round(MAE(standardize(test_x),test_y,t),2)}")
print(f"平均絶対値誤差率(MAPE)={round(MAPE(standardize(test_x),test_y,t),2)}%")
plt.show()これを実行するとこんな感じになります
x座標はz-score通した縮尺のまま表示してます
xを指定してyを取得したいときはfへ渡すxをz-score正規化をしてから渡します
分類
分類の理論(パーセプトロン)
分類は2通りの方法を紹介していきます
まずはパーセプトロンと呼ばれる線形的な手法から行います
ここではベクトルの知識が必要になります、説明がめんどいのでいきなり数式から入っちゃいますね()
グラフ(x_1,x_2)上にあるデータを2つの領域に直線で分類することを考えます
(分類先はyで表します、パーセプトロンではyは1か-1を取ります)
この直線の法線ベクトルをwとします
位置ベクトルxが直線に沿う場合を考えたとき、wと垂直なので内積を取ると
$$
w\cdot x=w_1x_1+w_2x_2=0
$$
という関係が成り立ちますね
ここで、分類を行う関数f_w(x)を内積の性質から定めます
$$
f_w(x)=\\
1\quad (w\cdot x\geq0)\\
-1\quad (w\cdot x<0)
$$
次に、パラメータwの更新式を定めます
これは、i番目のデータでf_w(x_i)とy_iが異なる場合のみwを次のように更新します(yは正解ラベル、すなわち正解となる分類先)
$$
w:=w+y_ix_i
$$
イメージとしては法線ベクトルwが回転していく感じです
分類の実装は次のロジスティクス回帰でやるので今回Pythonでパーセプトロンは実装しませんが、Processingでやったのがあって、こんな感じになります
扱いやすいモデルではあるんですけど、これだけだとあまり複雑な問題を扱えないわけです(ディープラーニングとかに応用されてるらしいので重要な考え方ではあると思います)
直線で分類できない問題のことを線形分離不可能とか言ったりします
分類の理論(ロジスティクス回帰)
パーセプトロンではできなかった線形分離不可能な問題にも対応できるように、分類を回帰によって行おうというのがロジスティクス回帰です
注意:ロジスティクス回帰では正解ラベルは0か1を取ります
しかし分類先は連続ではないですよね
そこで、こういう風に考えます
データがAである確率を調べる関数を作り、Aである確率が0.5以上ならA、以下ならAでないとする
確率は0~1の連続な値なので、回帰が使えるわけです
確率を表す関数として、次のシグモイド関数を使います
$$
f_\theta(x)=\frac{1}{1+\exp(-\theta^Tx)}
$$
θ^Txというのは、パラメータθとxを列ベクトルとみなしてかけているものです(すなわち、θ_0+θ_1*x_1+θ_2*x_2+…+θ_n*x_n)
シグモイド関数はこんな特徴があります
$$
\lim_{\theta^Tx\rightarrow\infty}f_\theta(x)=1\\\lim_{\theta^Tx\rightarrow-\infty}f_\theta(x)=0\\
\theta^Tx=0\leftrightarrow f_\theta(x)=0.5
$$
つまり、値域が(0,1)ということです、また、単調増加で、xが0のとき0.5になります
これは確率を示すにはちょうどいいわけです
そこで、Aである確率を
$$
P(y=1|x)=f_\theta(x)
$$
こう定義して
yを
$$
y=\\
1 \quad (f_\theta(x)\geq0.5)\\
0 \quad (f_\theta(x)<0.5)
$$
とします
回帰は最小値を探すものでしたが、パラメータの増減を逆にすれば最大値を探すこともできるはずです
ということで、確率が1になって欲しいものを掛け合わせた尤度関数というものを目的関数とします、この尤度関数の理想値は最大値(1)です
$$
L(\theta)=\prod_{i=1}^nP(y_i=1|x_i)^{y_i}P(y_i=0|x_i)^{1-y_i}
$$
ただ、これを微分しないといけないわけですが総乗の微分なんてしたくないわけです
総和に変換したいので、logを取ってみたらどうでしょう
$$
\log L(\theta)=\log\prod_{i=1}^nP(y_i=1|x_i)^{y_i}P(y_i=0|x_i)^{1-y_i} \\
=\sum_{i=1}^n\left(y_i\log P(y_i=1|x_i)+(1-y_i)\log P(y_i=0|x_i)\right) \\
P(y_i=0|x_i)=1-P(y_i=1|x_i)=f_\theta(x_i)なので \\
=\sum_{i=1}^n\left(y_i\log f_\theta(x_i)+(1-y_i)\log(1-f_\theta(x_i))\right)
$$
これで微分ができる形になりました
でも、あれ?別の関数の最大値を求めてどうするの?と思うかもしれません
たしかにこれで得られるのはlogを取った対数尤度関数の最大値ですが、欲しいのは最大値を取るパラメータなので、最大値は一致してなくても最大値を取るパラメータが一致してればいいわけです
logは単調増加の関数なので最大値を取るパラメータは変化しません
なのでこれでいいのです
微分の過程は省略して、最終的な結果だけ示します
$$
\frac{\partial}{\partial\theta_j}\log L(\theta)=\sum_{i=1}^n(y_i-f_\theta(x_i))x_{ij}
$$
そして、回帰の考え方でパラメータの更新式を定めれば
$$
\theta_j:=\theta_j+\eta\sum_{i=1}^n(y_i-f_\theta(x_i))x_{ij}
$$
これでロジスティクス回帰で分類ができるようになりました
ロジスティクス回帰の実験方法
ロジスティクス回帰では、回帰で行ったz-score正規化をx_1とx_2の両方について行います
また、表示する直線(決定境界)はθ^Tx=0で表されます、特に、x_1についてn次の多項式回帰を行い、x_2は1次なら
$$
x_2=-\frac{\theta_0+\theta_1x_1+\cdot\cdot\cdot+\theta_nx_n^n}{\theta_{n+1}}
$$
によって決定境界が得られます
x_1,x_2が2次なら円になったりします
(グラフにするのはちょっと面倒ですが、分かりやすいように図示してるだけで分類をするために必要なものではないので問題なし←)
回帰と同じように、3次関数に沿ったランダムなデータを生成しました(ただしノイズは含まない)
最後に分類の評価方法を紹介します
回帰を行ったモデルに対して、データをこの4つに分けます
まずは精度(Accuracy)です、これは単純に全体に対する正解の割合を示すものです
$$
Accuracy=\frac{TP+TN}{TP+TN+FP+FN}
$$
次に、適合率(Precision)と再現率(Recall)です
これはそれぞれ、Aと分類されたものの中で実際にAであった割合と、Aであるべきものの中でAと分類された割合を示す
$$
Precision=\frac{TP}{TP+FP}\\
Recall=\frac{TP}{TP+FN}
$$
また、PrecisionとRecallを総合的に評価する指標としてF値(Fmeasure)がある
$$
Fmeasure=\frac{2}{\frac{1}{Precision}+\frac{1}{Recall}}=\frac{2\cdot Precision\cdot Recall}{Precision+Recall}
$$
どの指標もそれぞれ最良が1、最悪が0となる
(ただし、Aと分類されたものがなかったりAであるべきものがなかったりしたらPrecisionとRecallはそもそも評価のしようがない、そうでない場合において、F値はPrecisionかRecallが0であったときに計算ができないが、これはlim[x->+0]1/xと同等に考えてよい(多分←)、つまり1/0=∞としてF値は2/∞=0になる)
Pythonでロジスティクス回帰
回帰と同じようにいきなりプログラム全体を貼ります←
テストはテスト用データのみで行うのではなく、データ全体で行います
ちなみにこれは前もって作ってなかったので記事書きながら作りました
import numpy as np
import matplotlib.pyplot as plt
from random import randint
data=np.loadtxt("bunruiData.csv",delimiter=",",skiprows=1) # データの読み込み
x1,x2,y=data[:,0],data[:,1],data[:,2] # x1,x2,yを取り出す
x1_0,x2_0,x1_1,x2_1=np.array([]),np.array([]),np.array([]),np.array([])
for i in range(x1.shape[0]): #ラベルごとにx1,x2を取り出す
if y[i]==0:
x1_0=np.append(x1_0,x1[i])
x2_0=np.append(x2_0,x2[i])
else:
x1_1=np.append(x1_1,x1[i])
x2_1=np.append(x2_1,x2[i])
test_x1,test_x2,test_y=np.array([]),np.array([]),np.array([]) # テストデータ用
for i in range(int(x1.shape[0]*3/10)): # データの3割をテスト用データとする
j=randint(0,x1.shape[0]-1) # ランダムなインデックス
test_x1=np.append(test_x1,x1[j])
test_x2=np.append(test_x2,x2[j])
test_y=np.append(test_y,y[j])
x1=np.delete(x1,j)
x2=np.delete(x2,j)
y=np.delete(y,j)
# x1のz-score正規化
mu=x1.mean() # 平均
sigma=x1.std() # 標準偏差
standardize=lambda x:(x-mu)/sigma # z-score正規化を行う関数
x1=standardize(x1)
# x2のz-score正規化
mu2=x2.mean()
sigma2=x2.std()
standardize2=lambda x:(x-mu2)/sigma2
x2=standardize2(x2)
plt.plot(standardize(x1_0),standardize2(x2_0),'o',color="b") # ラベルで色を分けてx1,x2を表示
plt.plot(standardize(x1_1),standardize2(x2_1),'o',color="r")
print("ロジスティクス回帰\nデータはbunruiData.csvより(z-score正規化を通す)\nx1の次数を入力:")
degree=int(input()) # 次数を受け取る
x=np.array([list(x1**i) for i in range(degree+1)]+[list(x2)])
# x、x0=x1^0=1,x1,x1^2,...,x1^degree,x2
t=np.array([1]*(degree+2)) # Θ、最後の要素はx2にかかるパラメータ
sgm=lambda x,theta:1/(np.exp(-np.dot(theta,x))+1) # シグモイド関数
print("学習率を入力:")
ETA=float(input()) # 学習率η
for i in range(10000):
t=t+ETA*np.dot(y-sgm(x,t),x.T)
graphx1=np.linspace(3,-3,100)
graphx2=-np.dot(np.array([graphx1**i for i in range(degree+1)]).T,t[:degree+1])/t[-1]
# x2=-(θ0*x1^0+...+θdegree*x_x1^degree)/θ(degree+1)
plt.plot(graphx1,graphx2,color="g") # 決定境界の表示
# 評価
TP,TN,FP,FN=0,0,0,0
x1=np.append(x1,standardize(test_x1))
x2=np.append(x2,standardize2(test_x2))
y=np.append(y,test_y)
x=np.array([list(x1**i) for i in range(degree+1)]+[list(x2)])
answer=sgm(x,t) # 得られた関数によるデータの分類の結果
for i in range(x1.shape[0]): # TP,TN,FP,FNを数える
if y[i]==1:
if answer[i]>=0.5:
TP+=1
else:
FP+=1
else:
if answer[i]<0.5:
TN+=1
else:
FN+=1
print(f"TP:{TP},TN:{TN},FP:{FP},FN:{FN}")
print(f"Accuracy={(TP+TN)/(TP+TN+FP+FN)}")
Precision=TP/(TP+FP)
print(f"Precision={Precision}")
Recall=TP/(TP+FN)
print(f"Recall={Recall}")
print(f"Fmeasure={2*Precision*Recall/(Precision+Recall)}")
plt.show()実行してみるとこんな感じです
まとめ
ということで今回はPythonで回帰と分類をやってみました
まだ正則化とかもありますけど記事がくっそ長くなってきたのでここまでにしておきます←
僕のnoteは他にも色々記事出してるのでよかったらそれも()
ここまで読んだ人いるんですか???
ということで、来てくれた人、ありがとうございました～
この記事が気に入ったらサポートをしてみませんか？