Skip to content

Commit caf666c

Browse files
KosaroKosaro
committed
Finished viterbi
1 parent b9fadd8 commit caf666c

File tree

1 file changed

+174
-24
lines changed

1 file changed

+174
-24
lines changed

‎HMM.ipynb‎

Lines changed: 174 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
},
2424
{
2525
"cell_type": "code",
26-
"execution_count": 118,
26+
"execution_count": 87,
2727
"outputs": [],
2828
"source": [
2929
"import numpy as np\n",
@@ -52,7 +52,7 @@
5252
},
5353
{
5454
"cell_type": "code",
55-
"execution_count": 119,
55+
"execution_count": 88,
5656
"outputs": [],
5757
"source": [
5858
"def get_pb1():\n",
@@ -73,6 +73,39 @@
7373
}
7474
}
7575
},
76+
{
77+
"cell_type": "markdown",
78+
"source": [
79+
"#### Emission"
80+
],
81+
"metadata": {
82+
"collapsed": false
83+
}
84+
},
85+
{
86+
"cell_type": "code",
87+
"execution_count": 89,
88+
"outputs": [],
89+
"source": [
90+
"def fair_die_emission():\n",
91+
" return 1/6\n",
92+
"\n",
93+
"def loaded_die_emission(value):\n",
94+
" if value == 6:\n",
95+
" return .5\n",
96+
" else:\n",
97+
" return .1\n",
98+
"\n",
99+
"def emission(value):\n",
100+
" return np.asarray((fair_die_emission(), loaded_die_emission(value)))"
101+
],
102+
"metadata": {
103+
"collapsed": false,
104+
"pycharm": {
105+
"name": "#%%\n"
106+
}
107+
}
108+
},
76109
{
77110
"cell_type": "markdown",
78111
"source": [
@@ -86,7 +119,7 @@
86119
},
87120
{
88121
"cell_type": "code",
89-
"execution_count": 120,
122+
"execution_count": 90,
90123
"outputs": [
91124
{
92125
"name": "stdout",
@@ -108,18 +141,6 @@
108141
}
109142
],
110143
"source": [
111-
"def fair_die_emission():\n",
112-
" return 1/6\n",
113-
"\n",
114-
"def loaded_die_emission(value):\n",
115-
" if value == 6:\n",
116-
" return .5\n",
117-
" else:\n",
118-
" return .1\n",
119-
"\n",
120-
"def emission(value):\n",
121-
" return np.asarray((fair_die_emission(), loaded_die_emission(value)))\n",
122-
"\n",
123144
"def viterbi(sequence):\n",
124145
" a = np.asarray([[.95,.05],[.05,.95]]) # transition probability\n",
125146
" b = emission(sequence[0]) # Emission probability\n",
@@ -163,31 +184,160 @@
163184
}
164185
},
165186
{
166-
"cell_type": "markdown",
187+
"cell_type": "code",
188+
"execution_count": 91,
189+
"outputs": [
190+
{
191+
"name": "stdout",
192+
"output_type": "stream",
193+
"text": [
194+
"0.39891275119173036\n"
195+
]
196+
}
197+
],
167198
"source": [
168-
"2. Download the dataset hmm_pb2.csv from Canvas. It represents a sequence of\n",
169-
"10000 dice rolls x from the Dishonest casino model but with other values for the a and\n",
170-
"b parameters than those from class. Having so many observations, you are going to\n",
171-
"learn the model parameters.\n"
199+
"def forward(sequence):\n",
200+
" tran = np.asarray([[.95,.05],[.05,.95]]) # transition probability\n",
201+
" b = emission(sequence[0]) # Emission probability\n",
202+
" p = np.asarray((.5,.5)) # Start probability\n",
203+
" a = np.ndarray((sequence.size, 2))\n",
204+
" a[0] = b*p\n",
205+
" for i in range(1,sequence.size):\n",
206+
" a[i] = emission(sequence[i]) * np.sum(a[i-1]*tran, axis=1)\n",
207+
" return a\n",
208+
"\n",
209+
"a = forward(get_pb1())\n",
210+
"print(a[124,0]/a[124,1])"
172211
],
173212
"metadata": {
174-
"collapsed": false
213+
"collapsed": false,
214+
"pycharm": {
215+
"name": "#%%\n"
216+
}
217+
}
218+
},
219+
{
220+
"cell_type": "code",
221+
"execution_count": 137,
222+
"outputs": [
223+
{
224+
"name": "stdout",
225+
"output_type": "stream",
226+
"text": [
227+
"3.856448201261194\n"
228+
]
229+
}
230+
],
231+
"source": [
232+
"def backward(sequence):\n",
233+
" tran = np.asarray([[.95,.05],[.05,.95]]) # transition probability\n",
234+
" B = np.ndarray((sequence.size, 2))\n",
235+
" B[-1] = np.ones(2)\n",
236+
" for i in range(sequence.size-2, -1, -1):\n",
237+
" B[i] = np.sum(tran*B[i+1]*emission(sequence[i+1]), axis=1)\n",
238+
" B[i] *=6 # multiply by constant to avoid overflow\n",
239+
" return B\n",
240+
"\n",
241+
"B = backward(get_pb1())\n",
242+
"print(B[124,0]/B[124,1])"
243+
],
244+
"metadata": {
245+
"collapsed": false,
246+
"pycharm": {
247+
"name": "#%%\n"
248+
}
175249
}
176250
},
177251
{
178252
"cell_type": "markdown",
179253
"source": [
254+
"2. Download the dataset hmm_pb2.csv from Canvas. It represents a sequence of\n",
255+
"10000 dice rolls x from the Dishonest casino model but with other values for the a and\n",
256+
"b parameters than those from class. Having so many observations, you are going to\n",
257+
"learn the model parameters.\n",
258+
"\n",
180259
"Implement and run the Baum-Welch algorithm using the forward and backward\n",
181260
"algorithms that you already implemented for Pb 1. You can initialize the $\\pi,a,b$ with\n",
182261
"your guess, or with some random probabilities (make sure that $\\pi$ sums to 1 and that\n",
183262
"$a_{ij}, b^i_k$\n",
184263
"sum to 1 for each $i$). The algorithm converges quite slowly, so you might need\n",
185264
"to run it for up 1000 iterations or more for the parameters to converge.\n",
186-
"Report the values of $\\pi,a,b$ that you have obtained. (4 points)\n",
187-
"\n"
265+
"Report the values of $\\pi,a,b$ that you have obtained. (4 points)\n"
188266
],
189267
"metadata": {
190-
"collapsed": false
268+
"collapsed": false,
269+
"pycharm": {
270+
"name": "#%% md\n"
271+
}
272+
}
273+
},
274+
{
275+
"cell_type": "code",
276+
"execution_count": 151,
277+
"outputs": [
278+
{
279+
"name": "stderr",
280+
"output_type": "stream",
281+
"text": [
282+
"c:\\users\\okosa\\pycharmprojects\\regressionasignment\\venv\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide\n",
283+
" \n"
284+
]
285+
},
286+
{
287+
"name": "stdout",
288+
"output_type": "stream",
289+
"text": [
290+
"pi [0.71477891 0.28522109]\n",
291+
"a [[2.72155818e-179 4.85324561e-181]\n",
292+
" [3.87594275e-180 2.70099412e-179]]\n",
293+
"b [[nan nan]\n",
294+
" [nan nan]\n",
295+
" [nan nan]\n",
296+
" [nan nan]\n",
297+
" [nan nan]\n",
298+
" [nan nan]]\n",
299+
"b sum [nan nan]\n"
300+
]
301+
}
302+
],
303+
"source": [
304+
"def baum_welch(sequence):\n",
305+
" for _ in range(20):\n",
306+
" # Expection step\n",
307+
" a = np.asarray([[.95,.05],[.05,.95]]) # transition probability\n",
308+
" A = forward(sequence)\n",
309+
" B = backward(sequence)\n",
310+
" emissions = np.asarray([emission(i) for i in range(1,7)])\n",
311+
" gamma = A*B / np.sum(A*B, axis=1)[np.newaxis].T\n",
312+
" eta = np.ndarray((sequence.size,2,2))\n",
313+
" for t in range(sequence.size):\n",
314+
" for i in range(2):\n",
315+
" for j in range(2):\n",
316+
" temp = A[t,i]*a[i,j]*B[t,j]*(emissions[int(sequence[t]-1)][j])\n",
317+
" eta[t,i,j] = temp\n",
318+
"\n",
319+
" # maximization step\n",
320+
" pi = gamma[0]\n",
321+
" for i in range(2):\n",
322+
" for j in range(2):\n",
323+
" a[i,j] = np.sum(eta[:,i,j])\n",
324+
" for i in range(2):\n",
325+
" for k in range(1,7):\n",
326+
" emissions[k-1,i] = np.sum(gamma[sequence==k, i])/np.sum(gamma)\n",
327+
" print(\"pi\", pi)\n",
328+
" print(\"a\", a)\n",
329+
" print(\"b\",emissions)\n",
330+
" print(\"b sum\",np.sum(emissions, axis=0))\n",
331+
"\n",
332+
"\n",
333+
"\n",
334+
"baum_welch(get_pb2())\n"
335+
],
336+
"metadata": {
337+
"collapsed": false,
338+
"pycharm": {
339+
"name": "#%%\n"
340+
}
191341
}
192342
}
193343
],

0 commit comments

Comments
 (0)