LoginSignup
3
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@tomson784

Juliaでベイズ分類器を実装（アヤメ分類）

Posted at

Juliaでベイズ分類器を実装（アヤメ分類）

プログラミング言語Juliaを使って、アヤメをベイズで分類していきます。

Julia Version 1.2.0

参考

ガウシアンナイーブベイズを自作する
図解・ベイズ統計「超」入門 (サイエンス・アイ新書)

ベイズの定理

トランプの例でベイズの定理を軽く説明します
ジョーカーを除いたトランプの総数 $13 \times 4=52$

クローバーを引き当てる確率　$\frac{1}{4}$

その条件の下でキングを引く確率 $\frac{1}{13}$

\frac{1}{4} \times \frac{1}{13} = \frac{1}{52}

クローバーのキングを引き当てる確率 $\frac{1}{52}$

よって、以下のような関係になっていることがわかります


p(y,x) = p(x|y)p(y) \\

p(y,x) = p(y|x)p(x) \\

p(y|x) = \frac{p(x|y)p(y)}{p(x)} \\

分類に必要なデータ

JuliaのRDatasetsからirisのデータを取得します。
RDatasetsは外部パッケージなのでインポートする必要があります。
必要なパッケージは適宜インポートしてください。

using RDatasets
iris = dataset("datasets","iris")

データセットは以下の表のようになっています

150×5 DataFrames.DataFrame
│ Row │ SepalLength │ SepalWidth │ PetalLength │ PetalWidth │ Species      │
│     │ Float64     │ Float64    │ Float64     │ Float64    │ Categorical… │
├─────┼─────────────┼────────────┼─────────────┼────────────┼──────────────┤
│ 1   │ 5.1         │ 3.5        │ 1.4         │ 0.2        │ setosa       │
│ 2   │ 4.9         │ 3.0        │ 1.4         │ 0.2        │ setosa       │
⋮
│ 148 │ 6.5         │ 3.0        │ 5.2         │ 2.0        │ virginica    │
│ 149 │ 6.2         │ 3.4        │ 5.4         │ 2.3        │ virginica    │
│ 150 │ 5.9         │ 3.0        │ 5.1         │ 1.8        │ virginica    │

表から、データがSepalLength、SepalWidth、PetalLength、PetalWidthという説明変数と Speciesというラベルデータの150セットのデータで構成されていることがわかります。

ここではラベルデータを$y$、説明変数を$\boldsymbol{x}$で表現します

出現するラベルの種類を取得

labels = unique(iris.Species)

# 3-element Array{String,1}:
#  "setosa"    
#  "versicolor"
#  "virginica"

よってラベルと説明変数のデータは以下のように表現します

y = {setosa, versicolor, virginica} \\
\boldsymbol{x} = (x_{SepalLength}, \ x_{SepalWidth}, \ x_{PetalLength}, \ x_{PetalWidth}) \\
X = (\boldsymbol{x_0},\boldsymbol{x_1},\boldsymbol{x_2},....\boldsymbol{x_i})

これらのデータを用いてアヤメを分類していきます

ベイズ分類に必要なもの

p(y|\boldsymbol{x}) = \frac{p(\boldsymbol{x}|y)p(y)}{p(\boldsymbol{x})}

このベイズの式によって、あるデータ$\boldsymbol{x}$が与えられた時の$y$の確率を
各ラベルごとに求めて比較し、もっとも確率の高かったラベルに分類します

p(y=setosa|\boldsymbol{x})  \quad\quad  データ\boldsymbol{x}が与えられた時、y=setosaである確率 \\ 
p(y=versicolor|\boldsymbol{x})  \qquad  データ\boldsymbol{x}が与えられた時、y=versicolorである確率 \\
p(y=virginica|\boldsymbol{x})  \qquad  データ\boldsymbol{x}が与えられた時、y=virginicaである確率

これらを比較し、最も数値の高くなったラベルに分類します。
この分類器を構築するために必要な要素を揃えていきます。

尤度

とあるラベル$y$における$\boldsymbol{x}$の確率は$p(\boldsymbol{x}|y)$で表現されます

とあるラベル$y$の時に変数$\boldsymbol{x}$がとある数値になる確率を求めるのですが、$\boldsymbol{x}$は離散値であり、どのような確率分布に従うかわかりません。

とりあえずデータの中身を確認してみます。確認のために各ラベルごとの各次元のデータをグラフ化します。

using Plots

a = []
for i in labels
	for j in 1:4
		push!(a, histogram(iris[iris.Species.==i, j] , bins=15))
	end
end

plot(
	a[1],a[2],a[3],a[4],
	a[5],a[6],a[7],a[8],
	a[9],a[10],a[11],a[12],layout=(3,4),legend=false
	)

ヒストグラムでプロットしたグラフです。
iris_param_histgram.png
なんとなく真ん中にピークが存在していて、ガウス分布に従っている気がします。
怪しいデータもありますがここでは母集団がガウス分布に従うものとして扱います。

データの特徴から最尤推定により、確率分布の形状を決定するパラメータを計算します

母集団がガウス分布と仮定して最尤推定で平均と分散を求めます。
ガウス分布の最尤推定に関しては(https://mathtrain.jp/mle)を参照してください

ガウス分布の式は以下のようになります

p(x|y) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\Biggl[ -\frac{(x-\mu)^2}{2\sigma^2} \Biggr]

gaussian(x, μ, σ²) = exp.((x - μ).^2 ./ 2σ²) ./ sqrt.(2*π*σ²)

説明変数の次元が多次元なので、多次元ガウス分布の公式でやる必要がありますが
共分散を求めたりなど、複雑な式が必要となってしまいます。
これを避けるために、各変数を独立なものとして扱うガウシアンナイーブベイズで考えます。
一応これでも、そこそこいい結果が得られるのでガウシアンナイーブベイズは使われるらしいです。

$\boldsymbol{x}$が独立の変数だということで以下のように分解できます。（独立の変数でなければ多次元ガウス分布として考える必要があります）

p(\boldsymbol{x}|y) = p(x_{SepalLength}|y)p(x_{SepalWidth}|y)p(x_{PetalLength}|y)p(x_{PetalWidth}|y)

簡潔に表現すると以下のようになります

p(\boldsymbol{x}|y) = \prod_{i=1}^{n} {p(x|y)}

ラベルごとの各説明変数のガウス分布を求める必要があるので
juliaで以下のように記述します

using Statistics

function average_label(label)
	mean_label = [
		mean(iris[iris.Species.==label, :SepalLength])
		mean(iris[iris.Species.==label, :SepalWidth])
		mean(iris[iris.Species.==label, :PetalLength])
		mean(iris[iris.Species.==label, :PetalWidth])
		]
	return mean_label
end

function variance_label(label)
	var_label = [
		var(iris[iris.Species.==label, :SepalLength])
		var(iris[iris.Species.==label, :SepalWidth])
		var(iris[iris.Species.==label, :PetalLength])
		var(iris[iris.Species.==label, :PetalWidth])
		]
	return var_label
end

means = [average_label(labels[1]), average_label(labels[2]), average_label(labels[3])]
# 3-element Array{Array{Float64,1},1}:
#  [5.006, 3.428, 1.462, 0.24600000000000002]           
#  [5.935999999999999, 2.77, 4.26, 1.3259999999999998]  
#  [6.587999999999999, 2.9739999999999998, 5.552, 2.026]

variance = [variance_label(labels[1]), variance_label(labels[2]), variance_label(labels[3])]
# 3-element Array{Array{Float64,1},1}:
#  [0.12424897959183674, 0.14368979591836734, 0.030159183673469387, 0.011106122448979593]
#  [0.2664326530612245, 0.09846938775510206, 0.22081632653061226, 0.03910612244897959]   
#  [0.4043428571428572, 0.10400408163265304, 0.30458775510204084, 0.07543265306122447]

これでラベルごとの各説明変数に対する平均と分散を求めたので、
各ガウス分布も導出することができます

事前分布

そもそも$y$自体はなんの条件もなしで、どれくらいの確率で発生するのかを確率$p(y)$で表現します

単純に特定のラベルの数を全体のラベル数で割り算すれば求まります

$p(\boldsymbol{x})$は$p(y|\boldsymbol{x})$を確率として扱うための正規化定数です。
今回のベイズ分類ではあるデータxが与えられた時の $y=setosa,versicolor,virginica$を
比較できれば良いので、ラベルごとの

p(y|\boldsymbol{x}) = \frac{p(\boldsymbol{x}|y)p(y)}{p(\boldsymbol{x})} \\ 
    \propto p(x|y)p(y)

を求めて大きさを比較すればよいです
ここでは、ラベルごとの事前分布は以下のようにします

prior_probability(label) = sum(iris.Species.==label)/length(iris.Species)

p_y = [
	prior_probability(labels[1]) # "setosa"
	prior_probability(labels[2]) # "versicolor"
	prior_probability(labels[3]) # "virginica"
]

# 3-element Array{Float64,1}:
#  0.3333333333333333
#  0.3333333333333333
#  0.3333333333333333

Juliaで実装

juliaの計算では行列の要素ごとの計算や比較、代入では
通常の計算演算子に $.$ をつけることで要素ごとの計算ができます

ex) .==   .*   ./ .=


A*B 行列積
A.*B アダマール積

以下のコードはベイズのアヤメ分類をひとまとめにしたものです。

using RDatasets
using Statistics
using Plots

iris = dataset("datasets", "iris")

labels = unique(iris.Species)

prior_probability(label) = sum(iris.Species.==label)/length(iris.Species)
p_y = [
	prior_probability(labels[1]) # "setosa"
	prior_probability(labels[2]) # "versicolor"
	prior_probability(labels[3]) # "virginica"
]

function average_label(label)
	mean_label = [
		mean(iris[iris.Species.==label, :SepalLength])
		mean(iris[iris.Species.==label, :SepalWidth])
		mean(iris[iris.Species.==label, :PetalLength])
		mean(iris[iris.Species.==label, :PetalWidth])
		]
	return mean_label
end

function variance_label(label)
	var_label = [
		var(iris[iris.Species.==label, :SepalLength])
		var(iris[iris.Species.==label, :SepalWidth])
		var(iris[iris.Species.==label, :PetalLength])
		var(iris[iris.Species.==label, :PetalWidth])
		]
	return var_label
end

means = [average_label(labels[1]), average_label(labels[2]), average_label(labels[3])]
variance = [variance_label(labels[1]), variance_label(labels[2]), variance_label(labels[3])]

gaussian(x, μ, σ²) = exp.((x - μ).^2 ./ 2σ²) ./ sqrt.(2*π*σ²)
label_probability = zeros(length(labels))

true_label = zeros(Int64, size(iris)[1])

for i in 1:size(iris)[1]
	for j in [1,2,3]
		x = [iris[i, :SepalLength], iris[i, :SepalWidth], iris[i, :PetalLength], iris[i, :PetalWidth]]
		label_probability[j] = p_y[j] * prod(gaussian(x, means[j], variance[j]))
	true_label[i] = findall(label_probability .== minimum(label_probability))[1]
	end
end

count_true = 0
for i in 1:1:size(iris)[1]
	global count_true
	if labels[true_label[i]] == iris.Species[i]
		count_true += 1
	end
end

println(count_true/length(iris.Species))

結果、精度は $0.94$ と表示されました。
Juliaは初めてだったのであまりきれいに記述することができませんでした。
そのうちどのようなデータにも対応できるような関数にしてみたいです。

3
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?