MyEnigma

とある自律移動システムエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

EMアルゴリズムによる確率分布学習のMATLABプログラム

f:id:meison_amsl:20160113193417p:plain

はじめに

ロボティクスにおいて、

あるデータを複数の確率分布で近似して、

それらのパラメータを推定したい時があります。


例えば、『確率ロボティクス』における

確率ロボティクス (Mynavi Advanced Library)

確率ロボティクス (Mynavi Advanced Library)


レーザレンジファインダのビームモデルは、

一本のビームの観測値を4つの確率分布の

足しあわせによって表現して、観測尤度を計算します。


もしこのビームモデルを使って、

観測尤度を計算したい場合、

それぞれの確率分布のパラメータや、

それぞれの確率分布の重み行列を、

それぞれのセンサに合わせてチューニングする必要があります。



また、名著 PRMLにもある通り、

パターン認識と機械学習 上

パターン認識と機械学習 上

  • 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
  • 出版社/メーカー: 丸善出版
  • 発売日: 2012/04/05
  • メディア: 単行本(ソフトカバー)
  • 購入: 6人 クリック: 33回
  • この商品を含むブログ (20件) を見る
パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

  • 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
  • 出版社/メーカー: 丸善出版
  • 発売日: 2012/02/29
  • メディア: 単行本
  • 購入: 6人 クリック: 14回
  • この商品を含むブログを見る


ある観測データの中から、

それぞれの物体からの観測値を

クラスタリングをしたい時など、

複数の正規分布でデータを近似し、

そのパラメータを推定したい時があります。

混合モデルとEMアルゴリズム(PRML第9章)


以上のような、

あるデータを複数の確率分布で近似し、

その最適なパラメータを推定したい時の

有効な方法の一つに

EMアルゴリズム あります。

EMアルゴリズム - Wikipedia



今回は、このEMアルゴリズムの簡単な説明と

冒頭の図のような、

2つの混合正規分布のパラメータを推定する

EMアルゴリズムのMATLABサンプルプログラムを公開したいと思います。

EMアルゴリズム

EMアルゴリズムは、一言で言うと、

確率分布のパラメータを、

尤度関数を最大化することで計算する方法です。


尤度関数である確率分布の期待値(Expectation)を

最大化(Maximization)することから、

EMアルゴリズムと呼ばれています。


大まかな流れとしては、

求めたいパラメータの初期値を設定して、

その値から尤度(期待値)を計算して、

多くの場合、

尤度関数の偏微分が0になる条件を使って、

繰り返し計算で最大尤度のパラメータを

計算するという方法です。


EMアルゴリズム自体は、

かなり汎用的なものなので、

ロボティクスだけでなく、音声認識や画像処理などにも使われています。


一次元の分布において、

2つのガウス分布をEMアルゴリズムで近似する方法は

下記の計算式で得られます。

f:id:meison_amsl:20160113195121p:plain

参照: Elements of Statistical Learning: data mining, inference, and prediction. 2nd Edition.


今回のMATLABサンプルプログラムも

上記のアルゴリズムを実装した形になります。


一つ注意点として、

混合ガウス分布で近似を実施する場合、

推定するのは各ガウス分布のパラメータだけでなく、

混合率という各ガウス分布をどれだけの割合づつ

組み合わせるのかというのを表したパラメータも推定する必要があります。


EMアルゴリズムの詳細に関しては

下記の資料やPRMLの9章を読むことをおすすめします。

EMアルゴリズム - Wikipedia

MLAC2013 数式を使わずイメージで理解するEMアルゴリズム - Wolfeyes Bioinformatics beta



MATLABサンプルコード

下記がEMアルゴリズムによる確率分布学習のMATLABサンプルプログラムです。

% -------------------------------------------------------------------------
%
% File : EMAlgorithmSample.m
%
% Discription : Sample code to estimate probablistic parameters 
%               with EM algorithm.
%
% Environment : Matlab
%
% Author : Atsushi Sakai
%
% Copyright (c): 2014 Atsushi Sakai
%
% License : GPL Software License Agreement
% -------------------------------------------------------------------------

function EMAlgorithmSample()
close all;
clear all;

%正規分布のパラメータ1
muTrue1=5;
siTrue1=2;
nData1=300;

%正規分布パラメータからサンプルデータを作成
r1 = muTrue1 + siTrue1.*randn(nData1,1);

%正規分布のパラメータ1
muTrue2=10;
siTrue2=1;
nData2=300;

%正規分布パラメータからサンプルデータを作成
r2 = muTrue2 + siTrue2.*randn(nData2,1);

data=[r1;r2];%元データ

%パラメータ初期値
mu=[5 15];%平均値の初期値
sigma=[std(data) std(data)];%標準偏差の初期値 データから計算
z=0.5;%混合率

%EMAlgorithmによる正規分布パラメータの推定
oddsResults=[];
for i=1:100
    burdenRates=EStep(data,mu,sigma,z);
    [mu,sigma,z]=MStep(data, burdenRates);
    odds = CalcLogOdds(data,mu,sigma,z);
    oddsResults=[oddsResults odds];
end

figure(1)
plot(oddsResults)
grid on;

figure(2)

disp('Estimated Mean');
mu
disp('True Mean');
[muTrue1 muTrue2]

disp('Estimated STD');
sigma
disp('True Std');
[siTrue1 siTrue2]

disp('Weight Vector');
z

[f,x]=hist(data,10);hold off;
bar(x,f/trapz(x,f));hold on;

est=[];
for id=0:0.05:max(data)
    est=[est; [id (1-z)*Gauss(id,mu(1),sigma(1))+z*Gauss(id,mu(2),sigma(2))]];
end

plot(est(:,1),est(:,2),'-r');hold on;
grid on;


function [mu,sigma,z]=MStep(data, burdenRates)
% MStepの計算

%平均の計算
mu=[0 0];
sumb1=sum(1-burdenRates);
sumby1=(1-burdenRates)*data;

mu(1)=sumby1/sumb1;

sumb2=sum(burdenRates);
sumby2=(burdenRates)*data;

mu(2)=sumby2/sumb2;

%標準偏差の計算
sigma=[0 0];
sumbyu1=(1-burdenRates)*((data-mu(1)).^2);
sigma(1)=sqrt(sumbyu1/sumb1);

sumbyu2=burdenRates*((data-mu(2)).^2);
sigma(2)=sqrt(sumbyu2/sumb2);

z=sum(burdenRates./length(data(:)));


function burdenRates=EStep(data,mu,sigma,z)
%負担率を計算する関数
burdenRates=[];%負担率

for i=1:length(data(:,1))
    
    %混合確率の計算
    n = z * Gauss(data(i), mu(2), sigma(2));
    d = (1 - z) * Gauss(data(i), mu(1), sigma(1)) + n;
    burdenRates=[burdenRates n/d];
end

function odds = CalcLogOdds(data,mu,sigma,z)
%対数オッズを計算する関数
odds=0;
for i=1:length(data)
    g1=Gauss(data(i),mu(1),sigma(1));
    g2=Gauss(data(i),mu(2),sigma(2));
    odds=odds+log((1-z)*g1+z*g2);
end

function si=UpdateSigma(e,r,mu)
%標準偏差を更新する関数
tmp=0;
for ie=1:length(e)
    tmp=tmp+(e(ie)*(r(ie)-mu)^2);
end
si=sqrt(1/sum(e)*tmp);

function p=Gauss(x,mu,sigma)
%ガウス分布における確率を計算する関数
p=(1/sqrt(2*pi*sigma^2)*exp(-(x-mu)^2/(2*sigma^2)));

下記のGitHubリポジトリでも公開しています。

MATLABRobotics/MachineLearning/EMAlgorithm/EMAlgorithmSample.m at master · AtsushiSakai/MATLABRobotics


このプログラムは、

2つの正規分布の平均と標準偏差、重みパラメータを

EMアルゴリズムを使って学習し、

尤もらしい混合分布のパラメータを学習するものです。


このプログラムを起動すると下記のように

データに対する正規分布のフィッティングの推定結果が表示されます。

f:id:meison_amsl:20160113193417p:plain


この結果を見ると分かる通り、

繰り返し計算の結果、

正確な確率分布のパラメータを学習できていることがわかります。


また下記のように対数尤度の推移をみると、

学習するにつれて対数尤度が増加しているのがわかります。

f:id:meison_amsl:20160113194726p:plain

PythonによるEMアルゴリズムのサンプルプログラム

下記のリンク先のコードが参考になると思います。

qiita.com


MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com