data = List@90, 86, 84, 79, 77, 77, 70, 68, 67, 66, 65, 64, 63, 60, 59, 58, 54, 53, 51, 47, 47,
46, 44, 43, 41, 40, 40, 40, 40, 40, 38, 38, 37, 37, 37, 37, 36, 36, 36, 36, 36, 36,
36, 36, 36, 36, 35, 35, 34, 34, 33, 33, 33, 32, 32, 32, 32, 31, 31, 30, 30, 30, 30,
30, 30, 30, 30, 30, 30, 30, 30, 29, 29, 29, 29, 29, 29, 29, 28, 28, 28, 27, 27, 27,
26, 26, 26, 26, 26, 26, 26, 26, 25, 25, 24, 24, 24, 24, 24, 24, 24, 23, 23, 22, 22,
22, 22, 22, 22, 22, 22, 22, 22, 22, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 19, 19, 19, 19, 19,
18, 18, 18, 18, 17, 17, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 15, 15, 15, 15, 15,
14, 14, 14, 14, 14, 14, 13, 13, 13, 13, 13, 12, 11, 10, 10, 8, 8, 4, 0, 0, 0, 0D;
Hx - ΜL2
1
pNormal@x_, Μ_, Σ_D :=
ExpB2Π Σ
F;
2 Σ2
r1@x_, Λ1_, Λ2_, Μ1_, Σ1_, Μ2_, Σ2_, Μ3_, Σ3_D :=
Λ1 pNormal@x, Μ1, Σ1D
;
Λ1 pNormal@x, Μ1, Σ1D + Λ2 pNormal@x, Μ2, Σ2D + H1 - Λ1 - Λ2L pNormal@x, Μ3, Σ3D
r2@x_, Λ1_, Λ2_, Μ1_, Σ1_, Μ2_, Σ2_, Μ3_, Σ3_D :=
Λ2 pNormal@x, Μ2, Σ2D
;
Λ1 pNormal@x, Μ1, Σ1D + Λ2 pNormal@x, Μ2, Σ2D + H1 - Λ1 - Λ2L pNormal@x, Μ3, Σ3D
r3@x_, Λ1_, Λ2_, Μ1_, Σ1_, Μ2_, Σ2_, Μ3_, Σ3_D :=
1 - r1@x, Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3D - r2@x, Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3D;
updateParam@data_, param_D :=
ModuleB8Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3, Λ1New, Λ2New,
Μ1New, Σ1New, Μ2New, Σ2New, Μ3New, Σ3New, r1List, r2List, r3List<,
Λ1 = param@@1DD; Λ2 = param@@2DD;
Μ1 = param@@3DD; Σ1 = param@@4DD;
Μ2 = param@@5DD; Σ2 = param@@6DD;
Μ3 = param@@7DD; Σ3 = param@@8DD;
r1List = Function@r1@ð, Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3DD ž data;
r2List = Function@r2@ð, Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3DD ž data;
r3List = Function@r3@ð, Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3DD ž data;
Λ1New = Mean@r1ListD; Λ2New = Mean@r2ListD;
MeanAr1List * data2 E
Mean@r1List * dataD
- Μ1New2 ;
; Σ1New =
Μ1New =
Λ1New
Λ1New
MeanAr2List * data2 E
Mean@r2List * dataD
- Μ2New2 ;
; Σ2New =
Μ2New =
Λ2New
Λ2New
MeanAr3List * data2 E
Mean@r3List * dataD
- Μ3New2 ;
; Σ3New =
Μ3New =
1 - Λ1New - Λ2New
1 - Λ1New - Λ2New
List@Λ1New, Λ2New, Μ1New, Σ1New, Μ2New, Σ2New, Μ3New, Σ3NewD
F;
PP@x_, param_D :=
Module@8Λ1, Λ2, Μ1, Σ1, Μ2, Σ2, Μ3, Σ3<,
Λ1 = param@@1DD; Λ2 = param@@2DD;
Μ1 = param@@3DD; Σ1 = param@@4DD;
Μ2 = param@@5DD; Σ2 = param@@6DD;
Μ3 = param@@7DD; Σ3 = param@@8DD;
Λ1 pNormal@x, Μ1, Σ1D + Λ2 pNormal@x, Μ2, Σ2D + H1 - Λ1 - Λ2L pNormal@x, Μ3, Σ3D
D;
logP@x_, param_D := Log@PP@x, paramDD;
L@data_, param_D := Mean@Function@logP@ð, paramDD ž dataD;
2
handout141219-2.nb
do := Module@8<,
param = updateParam@data, paramD;
L@data, paramD
D
H* パラメータの初期値を以下の様に設定する *L
param = 80.6, 0.2, 30, 10, 60, 10, 10, 10<
80.6, 0.2, 30, 10, 60, 10, 10, 10<
H* このときの対数尤度は以下の通り *L
L@data, paramD
- 4.13968
result = 8L@data, paramD<
8- 4.13968<
H* EMのパラメータ更新を 32回実行して対数尤度の変化を見る *L
Do@result = Append@result, doD, 832<D;
result
8- 4.13968, - 3.92836, - 3.92528, - 3.92376, - 3.92246, - 3.92094,
- 3.91852, - 3.91389, - 3.90533, - 3.89402, - 3.88357, - 3.87513,
- 3.86892, - 3.86463, - 3.8614, - 3.85835, - 3.85465, - 3.84929, - 3.84083,
- 3.82812, - 3.81272, - 3.79879, - 3.78935, - 3.7836, - 3.77897, - 3.77374,
- 3.76742, - 3.76166, - 3.75776, - 3.75442, - 3.74925, - 3.73255, - 3.6394<
ListPlot@resultD
-3.7
-3.8
-3.9
-4.0
-4.1
5
10
15
20
25
30
H* 32回更新後のパラメータ値は以下の通り *L
param
80.747341, 0.0919148, 24.559, 9.79087, 66.9869, 12.3513, 20.0245, 0.169311<
H* もう一回更新しようとすると、以下の様に破綻する *L
do
General
::
unfl
:計算中にアンダーフローが起りました． ‡
General
::
unfl
:計算中にアンダーフローが起りました． ‡
General
::
unfl
:計算中にアンダーフローが起りました． ‡
General
::
stop
:この計算中に，General
::
unfl
のこれ以上の出力は表示されません． ‡
- 2.64693
handout141219-2.nb
param = 80.7473413500794908`, 0.09191483338190157`,
24.558993907228547`, 9.790868423792269`, 66.98694624154335`,
12.35134407851044`, 20.024485068488314`, 0.16931052209272257`<
80.747341, 0.0919148, 24.559, 9.79087, 66.9869, 12.3513, 20.0245, 0.169311<
Plot@PP@x, paramD, 8x, 0, 100<, PlotRange ® 80, 0.05<D
0.05
0.04
0.03
0.02
0.01
0
20
40
60
80
100
Plot@PP@x, paramD, 8x, 0, 100<, PlotRange ® 80, 0.50<D
0.5
0.4
0.3
0.2
0.1
0
20
L@data, paramD
- 3.6394
40
60
80
100
3

(* 自然数 x を2進展開で長さnの2進列に変換する *) binarySeq[x_, n_