ベイズ推定と機械学習に関する話題

Wolfram Research Asia Limited 　丸山耕司

Join the Conversation @WolframJapan

ベイズ推定

githubで公開中のベイズ推定パッケージを$UserBaseDirectory/Applications/に事前に格納の必要あり
(HyperlinkURL["https://github.com/ssmit1986/BayesianInference", https://github.com/ssmit1986/BayesianInference])

データを元に予想する上での不確定要素

既知

データ

不確定

モデルの選択

モデルを特徴づけるパラメータ

モデルに固有の(予測に付随する)不確定性

簡単な例：データのフィッティング

ListPlot[testdata,PlotStyle->Directive[Red,PointSize[0.01]]]

多項式でフィッティングする場合のモデル例

y= a + b x + ϵ
ϵ \[Distributed] NormalDistribution[0, σ]
θ = {a, b, σ}

y= a + b x + c (x^2) + ϵ
ϵ \[Distributed] NormalDistribution[0, σ]
θ = {a, b, c, σ}

...

通常のフィッティング法(最尤推定)

線形フィットしたときの残差

residuals = testdata[[All,2]]-(a +b testdata[[All,1]])

残差の分布を正規分布 (0, σ)としたときの対数尤度

logLikelihood=Simplify[
LogLikelihood[NormalDistribution[0,σ],residuals]
]

尤度を最大化するa, b, σ.

{max,fit}=Maximize[
{logLikelihood, σ>0},
{a,b,σ}
]
y == a + b x +NormalDistribution[0,σ]/.fit

フィットした線形モデル

Plot[Evaluate[a+b x+InverseCDF[NormalDistribution[0,σ],{0.01,0.5,0.99}]/.fit],{x,-2.5,2.5},Epilog->{Red,PointSize[0.01],Point[testdata]},ImageSize->500,
PlotLabels->{"1%", "median", "99%"},PlotLabel->Style[StringForm["Loglikelihood: `1`", max],Large]
]

リーズナブルな推定に見えるが，a, b, σについての誤差情報は何もなし．

a, b, σにも誤差を考慮することはできるが，予測にまですべてを組み入れることはなかなか難しい．

同じデータを2次曲線でフィッティング

residuals2 = testdata[[All,2]]-(a +b testdata[[All,1]] + c ((testdata[[All,1]])^2));
logLikelihood2=Simplify[
LogLikelihood[NormalDistribution[0,σ],residuals2]
];
{max2,fit2}=Maximize[
{logLikelihood2, σ>0},
{a,b,c,σ}
];
Plot[Evaluate[a+b x+c (x^2)+ InverseCDF[NormalDistribution[0,σ],{0.01,0.5,0.99}]/.fit2],{x,-2.5,2.5},Epilog->{Red,PointSize[0.01],Point[testdata]},ImageSize->500,
PlotLabels->{"1%", "median", "99%"},PlotLabel->Style[StringForm["Loglikelihood: `1`", max2],Large]
]

尤度が高いという理由だけで2次のフィッティングの方が優れているのだろうか?

2次のフィッティングによる外挿を信じられるか?

ベイズ推定とNested Sampling

ベイズ推定では，条件付き確率に関するベイズの定理をもとに議論が構築される．

ベイズの定理

\[ProbabilityPr][Data\[Conditioned]θ,Model]*\[ProbabilityPr][θ\[Conditioned]Model]=\[ProbabilityPr][θ\[Conditioned]Data,Model]*\[ProbabilityPr][Data\[Conditioned]Model]

likelihood *Prior = Posterior * Evidence

尤度：モデルとそのパラメータが仮定されたときに，データが観測される確率

likelihood:  L= \[ProbabilityPr][Data\[Conditioned]θ,Model]

尤度は重要な指標だが，過信は禁物(前ページの例)．
目の前のデータをもっと直接活用して，未知の事象の確率を計算する方が有用！

想定しているモデルがデータを適切に説明しているか

Evidence: Z=\[ProbabilityPr][Data\[Conditioned]Model]=∫L*\[ProbabilityPr](θ\[Conditioned]Model)dθ

モデルのパラメータについての事後分布

Posterior = \[ProbabilityPr][θ\[Conditioned]Data,Model]

観測されたデータとモデルから未来のデータを予測

Prediction = \[ProbabilityPr][new data\[Conditioned]observed data,Model]

すべてベイズの定理を使って計算される．が，一般には計算量が膨大になって非常に厄介．

サンプリング

マルコフ連鎖モンテカルロ法 (MCMC)

ギブスサンプリング

Evidenceの計算 + サンプリング

Nested Sampling (今回のトピック)

Zには多くの情報が詰まっている．

Zを使ったベイズ因子により，異なるモデルの確からしさの比較に使える(古典的な統計学の仮説検定に代わる手法)

Nested samplingの基本アイデア

J. Skilling, “Nested Sampling for General Bayesian Computation”, Bayesian Analysis 1, 833 (2006).

目的: の計算

パラメータθ∈(((Reals))^d)の事前分布はπ(θ)とし，尤度Lはθの関数とする.

基本戦略: 事前確率massとθとの関係により，一次元空間での積分に直す．ただし，θ空間内の点とXの値の対応はわからないので，統計的な近似を行う．

（2次元の場合の尤度分布のイメージ)

(1) Prior mass Xを，尤度がλより高いθの領域での積分として与える．

この逆関数を積分すれば

となってevidenceが得られる．

(2) θ空間で事前分布π(θ)にしたがってランダムにN個の点をサンプリングすると，各点での尤度は容易にわかる．測度は正確にはわからないが，logXに対してほぼ線形に分布することは言える．

(3) よって，θ空間でN個サンプリングし，N個のXを[0,1]内の一様分布からサンプリングする．(Subscript[いちばん低い尤度L, 1])(Subscript[の点といちばん大きなXの点X, 1])が対応，(Subscript[さらにθ空間内でL, 1])以上の領域で新たにひとつの点((Subscript[L, 2]))，X∈[0,(Subscript[X, 1])]でひとつサンプリング((Subscript[X, 2]))，...を繰り返すと，の各項が得られる．

Nested samplingの動作例(2次元θ空間)

下の図はθ空間(モデルパラメータ)

事前分布にしたがってN個の点をサンプル(ここではN=30)し，z=0とする

尤度がもっとも低い点を見つけ，その点での尤度(Subscript[L, 1])とΔ(Subscript[X, 1])の積をzに加算，i.e., z = z + (Subscript[L, 1])*(Subscript[ΔX, 1])

(Subscript[L, 1])よりも高い尤度の領域で新しい点をひとつサンプル，その領域でもっとも低い尤度((Subscript[L, 2]))を用いてzを更新．z = z + (Subscript[L, 2])*(Subscript[ΔX, 2])

同様の作業を繰り返す(zの値が安定するまで)

Nested samplingの実装例

Mathematicaパッケージはここから入手可能
https://github.com/ssmit1986/BayesianInference

回帰曲線モデルの比較

上と同じテストデータ．Predictに与える教師データのように，associationの形{x1->y1,x2->y2,...}に書き直し．

testdata=Table[testdata[[i,1]]->testdata[[i,2]],{i,Length[testdata]}];
ListPlot[List@@@testdata,PlotStyle->Directive[Red,PointSize[0.01]]]

モデルは独立な変数に依存する(任意の)確率分布として与える

obj=defineInferenceProblem[
"Data"-> testdata,
"GeneratingDistribution" -> NormalDistribution[a x + b,σ],
"Parameters"->{{a,-10,10},{b,-10,10},{σ,0.1,10}},
"IndependentVariables"-> {x},
"PriorDistribution"-> {"LocationParameter","LocationParameter", "ScaleParameter"}
]
obj["Data"]

関数predictiveDistributionとregressionPlot1Dにより結果を可視化する

result=nestedSampling[obj]
predictions=predictiveDistribution[takePosteriorFraction[result,0.99],Range[-2,2,0.1]];

regressionPlot1D[result,predictions,"DistributionPercentiles"->{0.95,0.5,0.05}]

result["LogEvidence"]

calculationReport[result]//Quiet

1次，2次，3次式でフィッティングした結果を比較

variables[order_]:={
Sequence@@Table[{Symbol["c"<>ToString[i]],-100,100},{i,0,order}],
{s0,0.01,100}
};
result3=
Association@Table[
order-> With[{
obj=With[{
vars = variables[order]
},
defineInferenceProblem[
"Data"-> testdata,
"GeneratingDistribution" -> NormalDistribution[
vars[[;;-2, 1]] . x^Range[0,order],
vars[[-1,1]]
],
"Parameters"->vars,
"IndependentVariables"-> {x},
"PriorDistribution"-> Append[ConstantArray["LocationParameter",order+1], "ScaleParameter"]
]
]
},
nestedSampling[
obj,
"MonteCarloSteps"->300,
"SamplePoolSize"->25,
"MaxIterations"->3000
]
],
{order,{1,2,3}}]

regressionPlot1D[#,predictiveDistribution[#,Range[-2,2,0.1]]]&/@Map[
takePosteriorFraction[#,0.99]&,
result3
]

ニューラルネットワークによる(近似)ベイズ推定

基本アイデア:

Yarin Gal, Uncertainty in Deep Learning, 2016, PhD thesis.  (HyperlinkURL["http://www.cs.ox.ac.uk/people/yarin.gal/website/thesis/thesis.pdf", http://www.cs.ox.ac.uk/people/yarin.gal/website/thesis/thesis.pdf])

ニューラルネットにDropoutLayersを加え，ネットワーク自体にランダム性をもたせる

L2正規化を用いて訓練 (事前確率分布のパラメータを決めていることに相当)

予測(サンプリング)を行うにあたりDropoutLayersをアクティブにしておき(NetEvaluationMode -> “Train”)，ネットワークを繰り返し実行することで「予測エラー」を生成する

この手法の利点:

ガウス過程の回帰を行う場合のように，式の形を事前に指定する必要なし

Wolfram言語のニューラルネットワークで容易に実装可能

大きいサイズのデータに対しても良好なパフォーマンス

ネットワークの構築

各xの値に対してyの値とエラーバーを与えるネットワークを作りたい

BlockRandom[
SeedRandom[1];
testdata={#,(#+2 RandomReal[])}&/@RandomReal[{-1,1},20];
]

λ2 = 0.01; (* L2正規化係数 *)
pdrop = 0.5; (* dropoutの確率 *)
nUnits = 50; (* 隠れ層の次元 *)
activation =Ramp; (* 非線形要素 *)
regressionNet= NetGraph[
<|
"regression"-> NetChain[{
LinearLayer[nUnits],(* 入力 x を高次元ベクトルへマップ *)
ElementwiseLayer[activation],
DropoutLayer[pdrop],
LinearLayer[nUnits],
ElementwiseLayer[activation],
DropoutLayer[pdrop]
}],
"number1"-> LinearLayer[], (* 出力ベクトルをスカラーに変換 *)
"number2"->LinearLayer[]
|>,
{
NetPort["Input"]->"regression",
"regression"->"number1"-> NetPort["Mean"],
"regression"->"number2"-> NetPort["LogPrecision"]
},
"Input"->"Real",
"Mean"->"Real",
"LogPrecision"->"Real"  (* logτ = Log[1/(σ^2)]*)
]

ネットワークを訓練

訓練では，(観測された)データのnegative log-likelihoodを最小化する  ( = log-likelihoodを最大化 =  観測されたデータyの出現確率を最大化するμとσとなるようにNNを訓練)

FullSimplify[-2LogLikelihood[NormalDistribution[yPredicted,σ],{yObserved}]/.σ-> (1/(Sqrt[Exp[logτ]])),Assumptions->logτ∈Reals] (* τ=(σ^(-2)) is precision, which is large when σ is small (thus concentrated). τ is also (sort of) minimised, thus σ may be made large *)

loss=Function[{yObserved,yPredicted,logPrecision},
(((yPredicted-yObserved))^2)*Exp[logPrecision]-logPrecision
];
lossNet =NetInitialize@ NetGraph[<|
"reg"-> regressionNet,
"negLoglikelihood"->ThreadingLayer[loss]
|>,
{
NetPort["x"]->"reg",
{NetPort["y"],NetPort[{"reg","Mean"}],NetPort[{"reg","LogPrecision"}]}-> "negLoglikelihood"->NetPort["Loss"]
},
"y"->"Real",
"Loss"->"Real"
]

trainedNet = NetTrain[
lossNet,
<|"x"->testdata[[All,1]],"y"->testdata[[All,2]]|>,
LossFunction->"Loss",
Method->{"ADAM","L2Regularization"->λ2},
TimeGoal->20
];

予測

回帰を行うネットワークの部分を抽出

trainedRegressionNet=NetExtract[trainedNet,"reg"]

ネットワークにxの値のリストを一度与えるとこんな感じ:

xValues=Range[-2,2,0.025];
singleSample = trainedRegressionNet[xValues,NetEvaluationMode->"Train"];
Short/@singleSample

errorPlot[mean_List,stdev_List,xValues_List,col:_?ColorQ:Blue]:=ListPlot[
TimeSeries[
Transpose[{
mean-stdev,
mean,
mean+stdev
}],
{xValues},
ValueDimensions->3
],
errorStyle[col]
];
Show[
errorPlot[singleSample["Mean"],Exp[-0.5*singleSample["LogPrecision"]],xValues],
ListPlot[testdata,PlotStyle->Red],
ImageSize->Large
]

DropoutLayerをActiveにしたまま何度もネットワークによりサンプリングすることで，出力(予測)の統計をとることができる

sampleNet[net :(_NetChain|_NetGraph),xvalues_List,sampleNumber_Integer?Positive]:=
With[{regressionNet=NetExtract[net,"reg"]},
With[{
(* Evaluate the stochastic network a given number of times on the xvalues *)
samples = Select[Table[regressionNet[xvalues,NetEvaluationMode->"Train"],{i,sampleNumber}],AssociationQ]
},
<|
"Mean"-> Mean[samples[[All,"Mean"]]],
"StandardDeviation" -> Sqrt[Variance[samples[[All,"Mean"]]]+Mean[Exp[-samples[[All,"LogPrecision"]]]]]
|>
]
];

sampledPredictions = sampleNet[trainedNet,xValues,1000]; (* サンプリングする回数 *)
Short/@sampledPredictions

Show[
errorPlot[sampledPredictions["Mean"],sampledPredictions["StandardDeviation"],xValues],
ListPlot[testdata,PlotStyle->Red],
ImageSize->600,
PlotRange->All,
PlotLabel->"Prediction + 1σ error bands"
]

確率的でないネットワークによる予測と比べてみる(確率的でない = DropoutLayerをアクティブにしない)

defaultPrediction = trainedRegressionNet[xValues];
Short/@defaultPrediction

Legended[
Show[
ListPlot[testdata,PlotStyle->Red],
errorPlot[defaultPrediction["Mean"],Exp[-0.5*defaultPrediction["LogPrecision"]],xValues,Orange],
errorPlot[sampledPredictions["Mean"],sampledPredictions["StandardDeviation"],xValues,Blue],
ImageSize->Large,
PlotRange->All
],
SwatchLegend[{Blue,Orange},{"Sampled stochastic network","Standard prediction"}]
]

確率的でない分，NN中のパラメータについての不確定性が予測に含まれず，外挿部分でのエラーバンドの幅が狭い．

結局，ベイズ推定って何?

人間味ある推定をかなり(比較的)クールにやること

(完全に冷徹にやるわけではない)

Wolfram言語でお手軽に！

参考文献

John Skilling, Nested sampling for general Bayesian computation, Bayesian Anal. Volume 1, Number 4 (2006), 833-859, DOI:10.1214/06-BA127

Yarin Gal, Uncertainty in Deep Learning, 2016, PhD thesis. Available at (HyperlinkURL["http://www.cs.ox.ac.uk/people/yarin.gal/website/thesis/thesis.pdf", http://www.cs.ox.ac.uk/people/yarin.gal/website/thesis/thesis.pdf])

Nested Samplingを用いたベイズ推定パッケージ
(HyperlinkURL["https://github.com/ssmit1986/BayesianInference", https://github.com/ssmit1986/BayesianInference])

Initialisation

$HistoryLength=10;
<<BayesianInference`;
stockdata =Once[FinancialData["GE","Jan. 1, 2010"],"Local"];

BlockRandom[
SeedRandom[1];
testdata={#,(#+2 RandomReal[])}&/@RandomReal[{-1,1},20];
];

sdeTestData=;
residuals = testdata[[All,2]]-(a +b testdata[[All,1]]);
logLikelihood=Simplify[
LogLikelihood[NormalDistribution[0,σ],residuals]
];
errorStyle[col_?ColorQ] := Sequence[
Joined->True,
Filling->{1->{2},3->{2}},
PlotStyle->{Lighter[col],col,Lighter[col]}
];
errorPlot[mean_List,stdev_List,xValues_List,col:_?ColorQ:Blue]:=ListPlot[
TimeSeries[
Transpose[{
mean-stdev,
mean,
mean+stdev
}],
{xValues},
ValueDimensions->3
],
errorStyle[col]
];

nestedSamplingAnimation[result_,logLikelihood_]:=Module[{
samples=Values/@Values@Part[
SortBy[#LogLikelihood&]@result["Samples"],
All,{"Point","LogLikelihood"}],
subset,
nLiveSamples=result["SamplePoolSize"],
dist = result["GeneratingDistribution"]
},
Table[
subset=TakeDrop[Take[samples,UpTo[n]],-nLiveSamples];
ContourPlot[
logLikelihood[a,b],
{a,-10,10},{b,-10,10},
FrameLabel->(Style[#,"Text",16]&)/@{"a","b"},ImageSize->Large,
Epilog->
{
{Black,PointSize[0.01],Point[Most@subset[[2,All,1]]],PointSize[0.025],Point[Last@subset[[2,All,1]]]},
{Red,PointSize[0.01], Point[Most@subset[[1,All,1]]],PointSize[0.025],Point[Last@subset[[1,All,1]]]}

},
Contours->subset[[2,All,2]],
PlotLegends->Grid[{{Graphics[{Red,Disk[]},ImageSize->10],"Alive"},{Graphics[{Black,Disk[]},ImageSize->10],"Dead"}},BaseStyle->"Text",Alignment->Left]
],
{n,nLiveSamples+1,100}
]
];

loglikeFunction=Function[{a,b},Evaluate[logLikelihood/.σ->1]];
animationResult=nestedSampling[
loglikeFunction,
{"LocationParameter","LocationParameter"},
{{a,-10,10},{b,-10,10}},
"MonteCarloSteps"->300,
"SamplePoolSize"->30,
"MaxIterations"->5000
];

animation=nestedSamplingAnimation[animationResult,Function[{a,b},Evaluate[logLikelihood/.σ->1]]];

ListAnimate[animation,AnimationRate->1]