kees-study/cds-monte-carlo-methods
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

04: Try to create transition matrix P

Browse Source
This commit is contained in:
Kees van Kempen
2022-10-03 14:17:01 +02:00
parent 5152fbfc31
commit fb12323b0e
1 changed files with 39 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

48
Exercise sheet 4/exercise_sheet_04.ipynb
View File

@ -459,7 +459,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 10,
"id": "5de68351", "id": "5de68351",
"metadata": { "metadata": {
"deletable": false, "deletable": false,
@ -479,22 +479,39 @@
"def chain_P_histogram(graph,start,n):\n", "def chain_P_histogram(graph,start,n):\n",
" '''Produce a histogram of the states visited (excluding initial state) \n", " '''Produce a histogram of the states visited (excluding initial state) \n",
" by the P Markov chain in the first n steps when started at start.'''\n", " by the P Markov chain in the first n steps when started at start.'''\n",
" # YOUR CODE HERE\n", " n = graph.number_of_nodes()\n",
" raise NotImplementedError()\n", " number_of_visits = np.zeros(n)\n",
" x = start\n",
" \n",
" for _ in range(k):\n",
" x = sample_next_state(graph, x)\n",
" number_of_visits[x] += 1\n",
" \n",
" return number_of_visits\n",
"\n", "\n",
"def transition_matrix_P(graph):\n", "def transition_matrix_P(graph):\n",
" '''Construct transition matrix Q from graph as numpy array.'''\n", " '''Construct transition matrix Q from graph as numpy array.'''\n",
" # YOUR CODE HERE\n", " n = graph.number_of_nodes()\n",
" raise NotImplementedError()\n", " P = np.zeros((n, n))\n",
" Q = transition_matrix_Q(graph)\n",
" \n",
" for x in range(n):\n",
" for y in range(n):\n",
" if x == y:\n",
" P[y, x] = 1 - np.sum([Q[x, z]*acceptance_probability(graph, x, z) if z != x else 0 for z in range(n)])\n",
" else:\n",
" P[y, x] = Q[x, y]*acceptance_probability(graph, x, y)\n",
" \n",
" return P\n",
"\n", "\n",
"# plotting\n", "# plotting\n",
"# YOUR CODE HERE\n", "# YOUR CODE HERE\n",
"raise NotImplementedError()" "#raise NotImplementedError()"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 11,
"id": "4b5cc504", "id": "4b5cc504",
"metadata": { "metadata": {
"deletable": false, "deletable": false,
@ -521,7 +538,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 12,
"id": "f4e9d8aa", "id": "f4e9d8aa",
"metadata": { "metadata": {
"deletable": false, "deletable": false,
@ -538,7 +555,20 @@
"task": false "task": false
} }
}, },
"outputs": [], "outputs": [
{
"ename": "AssertionError",
"evalue": "11133.0 != 2222 within 180 delta (8911.0 difference)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
"Input \u001b[0;32mIn [12]\u001b[0m, in \u001b[0;36m<cell line: 2>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(chain_P_histogram(example_graph,\u001b[38;5;241m0\u001b[39m,\u001b[38;5;241m100\u001b[39m)) \u001b[38;5;241m==\u001b[39m example_graph\u001b[38;5;241m.\u001b[39mnumber_of_nodes()\n\u001b[0;32m----> 2\u001b[0m \u001b[43massert_almost_equal\u001b[49m\u001b[43m(\u001b[49m\u001b[43mchain_P_histogram\u001b[49m\u001b[43m(\u001b[49m\u001b[43mexample_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m20000\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m8\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m2222\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mdelta\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m180\u001b[39;49m\u001b[43m)\u001b[49m\n",
"File \u001b[0;32m/opt/jupyter-conda/lib/python3.9/unittest/case.py:891\u001b[0m, in \u001b[0;36mTestCase.assertAlmostEqual\u001b[0;34m(self, first, second, places, msg, delta)\u001b[0m\n\u001b[1;32m 885\u001b[0m standardMsg \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m != \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m within \u001b[39m\u001b[38;5;132;01m%r\u001b[39;00m\u001b[38;5;124m places (\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m difference)\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;241m%\u001b[39m (\n\u001b[1;32m 886\u001b[0m safe_repr(first),\n\u001b[1;32m 887\u001b[0m safe_repr(second),\n\u001b[1;32m 888\u001b[0m places,\n\u001b[1;32m 889\u001b[0m safe_repr(diff))\n\u001b[1;32m 890\u001b[0m msg \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_formatMessage(msg, standardMsg)\n\u001b[0;32m--> 891\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfailureException(msg)\n",
"\u001b[0;31mAssertionError\u001b[0m: 11133.0 != 2222 within 180 delta (8911.0 difference)"
]
}
],
"source": [ "source": [
"assert len(chain_P_histogram(example_graph,0,100)) == example_graph.number_of_nodes()\n", "assert len(chain_P_histogram(example_graph,0,100)) == example_graph.number_of_nodes()\n",
"assert_almost_equal(chain_P_histogram(example_graph,0,20000)[8],2222,delta=180)" "assert_almost_equal(chain_P_histogram(example_graph,0,20000)[8],2222,delta=180)"