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

f:id:meison_amsl:20160113193417p:plain

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

posted with カエレバ

C.M. ビショップ丸善出版 2012-04-05

Amazonで探す

楽天市場で探す

Yahooショッピングで探す

はじめに

ロボティクスにおいて、

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

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

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

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

作者: Sebastian Thrun,Wolfram Burgard,Dieter Fox
出版社/メーカー: マイナビ出版
発売日: 2015/04/02
メディア: Kindle版
この商品を含むブログを見る

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

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

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

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

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

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

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

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

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

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

作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
出版社/メーカー: 丸善出版
発売日: 2012/04/05
メディア: 単行本（ソフトカバー）
購入: 6人クリック: 33回
この商品を含むブログ (20件) を見る

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

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

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

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

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

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

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

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

以上のような、

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

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

有効な方法の一つに

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

EMアルゴリズム - Wikipedia

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

冒頭の図のような、

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

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

EMアルゴリズム

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

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

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

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

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

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

大まかな流れとしては、

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

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

多くの場合、

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

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

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

EMアルゴリズム自体は、

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

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

一次元の分布において、

２つのガウス分布を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

数式を使わずイメージで理解するEMアルゴリズム from yag ays

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)));