Skip to content

Instantly share code, notes, and snippets.

@trycycle

trycycle/pymc_irt.ipynb Secret

Last active April 27, 2023 05:18
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save trycycle/84bdf4ab40807a01386e1876f6c13b87 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save trycycle/84bdf4ab40807a01386e1876f6c13b87 to your computer and use it in GitHub Desktop.
Download ZIP
PyMCで項目反応理論を実装する練習
Display the source blob
Display the rendered blob
Raw
pymc_irt.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "478037e6-7e6e-449a-bb28-5a0cc5582d3f",
"metadata": {},
"source": [
"# PyMCによる項目反応理論の実装する練習\n",
"\n",
"[項目反応理論（IRT: Item Response Theory）](https://ja.wikipedia.org/wiki/%E9%A0%85%E7%9B%AE%E5%BF%9C%E7%AD%94%E7%90%86%E8%AB%96)は，TOEFLなどのような試験の各設問の難易度，および回答者の能力を推測するためのモデルである．基本的な考え方は，回答者の能力と設問の難易度をパラメータとするロジスティック関数をある設問の正解確率と見なし，試験の回答状況は定義された正解確率をパラメータとするベルヌーイ分布に従う，というものである．\n",
"\n",
"パラメータの推定方法は最尤推定法が一般的であるが，バージョン4にアップデートされた[PyMC](https://www.pymc.io/welcome.html)の練習をかねて，MCMCで項目反応理論のパラメータ推定をやってみた．"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "8a808602-6799-4bc8-b08f-12f5b83815d8",
"metadata": {},
"outputs": [],
"source": [
"import arviz as az\n",
"import pandas as pd\n",
"import numpy as np\n",
"import pymc as pm"
]
},
{
"cell_type": "markdown",
"id": "b8bb1dbf-670c-4013-8237-67690609db8e",
"metadata": {},
"source": [
"## データセットの読み込み\n",
"練習用データとして，20ある設問に対して2000名の受験者の回答結果を記したデータを用意した．\n",
"各行が各受験者の回答状況に対応しており，設問1（Q1）〜設問20（Q20）に対する回答結果が，正解1，不正解0として記されている．\n",
"このデータは疑似データであり，[試験の数理](https://qiita.com/takuyakubo/items/43d56725952e67032b49)を参考に生成した．\n",
"\n",
"[データ](https://raw.githubusercontent.com/trycycle/pymc_irt/main/dataset/item_response.tsv)の中身は以下の通りである．"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "0141ac3f-587f-41a4-ae57-8141e66a8442",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Q1</th>\n",
" <th>Q2</th>\n",
" <th>Q3</th>\n",
" <th>Q4</th>\n",
" <th>Q5</th>\n",
" <th>Q6</th>\n",
" <th>Q7</th>\n",
" <th>Q8</th>\n",
" <th>Q9</th>\n",
" <th>Q10</th>\n",
" <th>Q11</th>\n",
" <th>Q12</th>\n",
" <th>Q13</th>\n",
" <th>Q14</th>\n",
" <th>Q15</th>\n",
" <th>Q16</th>\n",
" <th>Q17</th>\n",
" <th>Q18</th>\n",
" <th>Q19</th>\n",
" <th>Q20</th>\n",
" </tr>\n",
" <tr>\n",
" <th>participant</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14 Q15 \\\n",
"participant \n",
"0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 \n",
"1 1 1 0 1 1 1 0 0 1 0 1 1 1 1 0 \n",
"2 1 0 0 1 1 0 1 1 1 1 0 1 1 0 0 \n",
"3 1 1 0 1 1 1 1 0 0 0 1 1 1 1 1 \n",
"4 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1 \n",
"\n",
" Q16 Q17 Q18 Q19 Q20 \n",
"participant \n",
"0 0 0 1 1 0 \n",
"1 1 0 0 1 0 \n",
"2 0 0 1 0 1 \n",
"3 1 0 1 1 1 \n",
"4 1 1 1 0 1 "
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.read_table('https://raw.githubusercontent.com/trycycle/pymc_irt/main/dataset/item_response.tsv', sep='\\t', header=0, index_col=0)\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"id": "68d02155-8bb1-4a06-a082-fb30d12a7400",
"metadata": {},
"source": [
"この疑似データは，項目反応理論の3パラメータロジスティックモデルで生成している．すなわち，回答者$i$の設問$j$に対する正解確率$p_{i,j}$は回答者の能力パラメータ$\\theta_i$と設問パラメータ$a_j$（識別力），$b_j$（難易度），$c_j$（当て推量）で以下のように表現できると仮定している．\n",
"\n",
"$$\n",
"p_{i,j} = c_j + (1-c_j) \\frac{1}{1+exp(- a_j(\\theta_i - b_j))}\n",
"$$\n",
"\n",
"生成に用いた設問のパラメータは`dataset/item_params.tsv`，回答者の能力パラメータは`dataset/participant_params.tsv`で確認できる．"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "75c21001-31ce-469f-9c04-f1d4ba645f26",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>b</th>\n",
" <th>c</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>Q1</th>\n",
" <td>0.530947</td>\n",
" <td>-1.126938</td>\n",
" <td>0.262671</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q2</th>\n",
" <td>0.778774</td>\n",
" <td>0.440435</td>\n",
" <td>0.016508</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q3</th>\n",
" <td>0.851570</td>\n",
" <td>1.997521</td>\n",
" <td>0.078272</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q4</th>\n",
" <td>0.496084</td>\n",
" <td>-1.853463</td>\n",
" <td>0.310627</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q5</th>\n",
" <td>0.592267</td>\n",
" <td>-0.531258</td>\n",
" <td>0.134943</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q6</th>\n",
" <td>0.827886</td>\n",
" <td>1.180413</td>\n",
" <td>0.169449</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q7</th>\n",
" <td>0.686584</td>\n",
" <td>0.559992</td>\n",
" <td>0.363648</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q8</th>\n",
" <td>0.735757</td>\n",
" <td>1.890666</td>\n",
" <td>0.340324</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q9</th>\n",
" <td>0.366514</td>\n",
" <td>0.357466</td>\n",
" <td>0.198492</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q10</th>\n",
" <td>0.504268</td>\n",
" <td>-1.591176</td>\n",
" <td>0.140980</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q11</th>\n",
" <td>0.800319</td>\n",
" <td>1.106090</td>\n",
" <td>0.331420</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q12</th>\n",
" <td>0.704560</td>\n",
" <td>-1.174947</td>\n",
" <td>0.192516</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q13</th>\n",
" <td>0.927804</td>\n",
" <td>1.391537</td>\n",
" <td>0.204997</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q14</th>\n",
" <td>0.744368</td>\n",
" <td>-0.900265</td>\n",
" <td>0.013639</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q15</th>\n",
" <td>0.464219</td>\n",
" <td>1.375144</td>\n",
" <td>0.398396</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q16</th>\n",
" <td>0.433757</td>\n",
" <td>0.358018</td>\n",
" <td>0.399029</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q17</th>\n",
" <td>0.598347</td>\n",
" <td>1.633047</td>\n",
" <td>0.206599</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q18</th>\n",
" <td>0.675967</td>\n",
" <td>0.070730</td>\n",
" <td>0.262606</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q19</th>\n",
" <td>0.431652</td>\n",
" <td>-1.348594</td>\n",
" <td>0.368219</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q20</th>\n",
" <td>0.537283</td>\n",
" <td>-0.668843</td>\n",
" <td>0.251276</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a b c\n",
"Q1 0.530947 -1.126938 0.262671\n",
"Q2 0.778774 0.440435 0.016508\n",
"Q3 0.851570 1.997521 0.078272\n",
"Q4 0.496084 -1.853463 0.310627\n",
"Q5 0.592267 -0.531258 0.134943\n",
"Q6 0.827886 1.180413 0.169449\n",
"Q7 0.686584 0.559992 0.363648\n",
"Q8 0.735757 1.890666 0.340324\n",
"Q9 0.366514 0.357466 0.198492\n",
"Q10 0.504268 -1.591176 0.140980\n",
"Q11 0.800319 1.106090 0.331420\n",
"Q12 0.704560 -1.174947 0.192516\n",
"Q13 0.927804 1.391537 0.204997\n",
"Q14 0.744368 -0.900265 0.013639\n",
"Q15 0.464219 1.375144 0.398396\n",
"Q16 0.433757 0.358018 0.399029\n",
"Q17 0.598347 1.633047 0.206599\n",
"Q18 0.675967 0.070730 0.262606\n",
"Q19 0.431652 -1.348594 0.368219\n",
"Q20 0.537283 -0.668843 0.251276"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 設問のパラメータを確認\n",
"item_df = pd.read_table('https://raw.githubusercontent.com/trycycle/pymc_irt/main/dataset/item_params.tsv', sep='\\t', header=0)\n",
"item_df.index = [f'Q{qid}' for qid in range(1, 21)]\n",
"item_df"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "1529b040-eeef-4463-99dc-a9a883a457c0",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>theta</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>-0.846466</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.215250</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0.367516</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0.283451</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>2.386955</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>0.300458</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>-0.204278</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>-0.249474</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>0.080900</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>1.191649</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" theta\n",
"0 -0.846466\n",
"1 -0.215250\n",
"2 0.367516\n",
"3 0.283451\n",
"4 2.386955\n",
"5 0.300458\n",
"6 -0.204278\n",
"7 -0.249474\n",
"8 0.080900\n",
"9 1.191649"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 回答者の能力パラメータ\n",
"participant_df = pd.read_table('https://raw.githubusercontent.com/trycycle/pymc_irt/main/dataset/participant_params.tsv', sep='\\t', header=0)\n",
"participant_df.head(10)"
]
},
{
"cell_type": "markdown",
"id": "fa877271-dbc4-4342-9f73-42b0f7888388",
"metadata": {},
"source": [
"上記回答者の能力パラメータ$\\theta$と設問パラメータ$a$，$b$，$c$が未知と仮定して，それらを回答結果データ（`df`）とPyMCを使って推定するのが今回のタスクとなる．"
]
},
{
"cell_type": "markdown",
"id": "6c0eebe6-0660-42be-9af8-dcbaf3acb2e1",
"metadata": {},
"source": [
"## データの前処理"
]
},
{
"cell_type": "markdown",
"id": "621d947a-8e85-48f4-a119-6bc44ad902eb",
"metadata": {},
"source": [
"生の回答結果データ`df`は横持ちのデータになっているので，PyMCで処理しづらい．そこで`df`に格納されたデータを整然データ（tidy data）化しておく．"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "b6e4bb63-23e9-4be7-949f-83e03afc7bf3",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>participant</th>\n",
" <th>question</th>\n",
" <th>response</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>Q1</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>Q1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>Q1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>Q1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>Q1</td>\n",
" <td>1</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" participant question response\n",
"0 0 Q1 0\n",
"1 1 Q1 1\n",
"2 2 Q1 1\n",
"3 3 Q1 1\n",
"4 4 Q1 1"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# tidy data化\n",
"response_df = pd.melt(df.reset_index(),\n",
" id_vars='participant', var_name='question', value_name='response')\n",
"response_df.head()"
]
},
{
"cell_type": "markdown",
"id": "e9eec3ca-09d1-4ecd-874c-caa904611bb6",
"metadata": {},
"source": [
"次に，`response_df`を構成する回答者情報，設問名をコード化して取り出しておく．この処理を行っておくことで，PyMCでモデルを書くときにfor文などを使わずに各種パラメータをベクトルとしてうまく扱えるようになる．"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "b4f2ca68-d5a4-4a75-a31d-74b19bef2630",
"metadata": {},
"outputs": [],
"source": [
"# participants[participant_idx] is equal to response_df['participant']\n",
"participant_idx, participants = pd.factorize(response_df['participant'])\n",
"\n",
"# questions[question_idx] is equal to response_df['question']\n",
"question_idx, questions = pd.factorize(response_df['question'])"
]
},
{
"cell_type": "markdown",
"id": "3fc1a6b8-30a9-4cfc-a966-5d5afbe326af",
"metadata": {},
"source": [
"`participant_idx`および`quesiton_idx`は，それぞれ`response_df`の各レコード（行）に対応する回答者および質問のIDが格納されている．\n",
"`participants`および`questions`はnumpyの配列であるため，ファンシーインデックスとしてが使える．"
]
},
{
"cell_type": "markdown",
"id": "b360a237-8179-46dd-8d7b-5528151ba572",
"metadata": {},
"source": [
"## モデリング"
]
},
{
"cell_type": "markdown",
"id": "54e4572f-7f86-4e11-8cf5-82b968f7a067",
"metadata": {},
"source": [
"データの準備が整ったので，モデルの定義，事後分布からのサンプリングを行う．今回用いるデータは3パラメータロジスティックモデルに従うようデータを生成したが，そのことは知らないふりをして **2パラメータロジスティックモデル** を仮定してモデリングを行う．\n",
"\n",
"項目反応理論では基本的なモデルとして，1パラメータ，2パラメータ，3パラメータのロジスティックモデルが提案されている．パラメータが増えるごとに複雑な状況を想定できる．2パラメータのロジスティックモデルでは，\n",
"* 回答者\n",
" - 能力パラメータ $\\theta$（実数）\n",
"* 設問\n",
" - 識別力パラメータ $a$（正の値）\n",
" - 難易度パラメータ $b$（実数）\n",
"\n",
"の存在を仮定している．その上で，回答者$i$の能力を$\\theta_i$，設問$j$の識別力，難易度，当て推量を$a_j$，$b_j$としたとき，回答者$i$の設問$j$に対する回答の正誤$r_{i,j}$は以下に従うと考える：\n",
"\n",
"$$\n",
"p_{i,j} = \\frac{1}{1 + exp(-a_j(\\theta_i - b_j))} \\\\\n",
"r_{i,j} \\sim Bernoulli(p_{ij})\n",
"$$\n",
"\n",
"これから，$a$，$b$，$\\theta$の事後分布をPyMCによって推定する．\n",
"パラメータの事前分布は適当に\n",
"\n",
"$$\n",
"a \\sim HalfNormal(1) \\\\\n",
"b \\sim N(0,1) \\\\\n",
"\\theta \\sim N(0, 1)\n",
"$$\n",
"\n",
"とすることにする．上記の仮定をPyMCでモデル化すると，下記のように書ける．"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "aafd2de6-0eba-44b4-8df0-21f072c845a2",
"metadata": {},
"outputs": [],
"source": [
"coords = {'participant': participants, 'question': questions}\n",
"model = pm.Model(coords=coords)\n",
"\n",
"with model: \n",
" # 個人の能力（回答者ごとに割り当て）\n",
" theta = pm.Normal('theta', mu=0, sigma=1, dims='participant')\n",
" # 問題の識別力（設問ごとに割り当て）\n",
" a = pm.HalfNormal('a', sigma=1, dims='question')\n",
" # 問題の難易度（設問ごとに割り当て）\n",
" b = pm.Normal('b', mu=0, sigma=1, dims='question')\n",
" \n",
" # 2パラメータロジスティックモデル（2PLM）\n",
" p = 1 / (1 + pm.math.exp(-1 * a[question_idx] * (theta[participant_idx] - b[question_idx])))\n",
"\n",
" y_obs = pm.Bernoulli('y_obs', p=p, observed=response_df.response)"
]
},
{
"cell_type": "markdown",
"id": "9f5813f3-3398-4f7b-bbfd-404d86bc7a19",
"metadata": {},
"source": [
"参考として，定義したモデルをDAGとして可視化しておく．"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "836b0e7c-6b0f-4f11-be9d-613ec3371eee",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
"<!-- Generated by graphviz version 3.0.0 (20220226.1711)\n",
" -->\n",
"<!-- Pages: 1 -->\n",
"<svg width=\"381pt\" height=\"260pt\"\n",
" viewBox=\"0.00 0.00 381.00 259.91\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
"<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 255.91)\">\n",
"<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-255.91 377,-255.91 377,4 -4,4\"/>\n",
"<g id=\"clust2\" class=\"cluster\">\n",
"<title>clusterquestion (20)</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M142,-129.95C142,-129.95 353,-129.95 353,-129.95 359,-129.95 365,-135.95 365,-141.95 365,-141.95 365,-231.91 365,-231.91 365,-237.91 359,-243.91 353,-243.91 353,-243.91 142,-243.91 142,-243.91 136,-243.91 130,-237.91 130,-231.91 130,-231.91 130,-141.95 130,-141.95 130,-135.95 136,-129.95 142,-129.95\"/>\n",
"<text text-anchor=\"middle\" x=\"320.5\" y=\"-137.75\" font-family=\"Times,serif\" font-size=\"14.00\">question (20)</text>\n",
"</g>\n",
"<g id=\"clust3\" class=\"cluster\">\n",
"<title>cluster40000</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M153,-8C153,-8 259,-8 259,-8 265,-8 271,-14 271,-20 271,-20 271,-109.95 271,-109.95 271,-115.95 265,-121.95 259,-121.95 259,-121.95 153,-121.95 153,-121.95 147,-121.95 141,-115.95 141,-109.95 141,-109.95 141,-20 141,-20 141,-14 147,-8 153,-8\"/>\n",
"<text text-anchor=\"middle\" x=\"246\" y=\"-15.8\" font-family=\"Times,serif\" font-size=\"14.00\">40000</text>\n",
"</g>\n",
"<g id=\"clust1\" class=\"cluster\">\n",
"<title>clusterparticipant (2000)</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M20,-129.95C20,-129.95 110,-129.95 110,-129.95 116,-129.95 122,-135.95 122,-141.95 122,-141.95 122,-231.91 122,-231.91 122,-237.91 116,-243.91 110,-243.91 110,-243.91 20,-243.91 20,-243.91 14,-243.91 8,-237.91 8,-231.91 8,-231.91 8,-141.95 8,-141.95 8,-135.95 14,-129.95 20,-129.95\"/>\n",
"<text text-anchor=\"middle\" x=\"65\" y=\"-137.75\" font-family=\"Times,serif\" font-size=\"14.00\">participant (2000)</text>\n",
"</g>\n",
"<!-- theta -->\n",
"<g id=\"node1\" class=\"node\">\n",
"<title>theta</title>\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"72\" cy=\"-198.43\" rx=\"41.94\" ry=\"37.45\"/>\n",
"<text text-anchor=\"middle\" x=\"72\" y=\"-209.73\" font-family=\"Times,serif\" font-size=\"14.00\">theta</text>\n",
"<text text-anchor=\"middle\" x=\"72\" y=\"-194.73\" font-family=\"Times,serif\" font-size=\"14.00\">~</text>\n",
"<text text-anchor=\"middle\" x=\"72\" y=\"-179.73\" font-family=\"Times,serif\" font-size=\"14.00\">Normal</text>\n",
"</g>\n",
"<!-- y_obs -->\n",
"<g id=\"node4\" class=\"node\">\n",
"<title>y_obs</title>\n",
"<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"197\" cy=\"-76.48\" rx=\"48.17\" ry=\"37.45\"/>\n",
"<text text-anchor=\"middle\" x=\"197\" y=\"-87.78\" font-family=\"Times,serif\" font-size=\"14.00\">y_obs</text>\n",
"<text text-anchor=\"middle\" x=\"197\" y=\"-72.78\" font-family=\"Times,serif\" font-size=\"14.00\">~</text>\n",
"<text text-anchor=\"middle\" x=\"197\" y=\"-57.78\" font-family=\"Times,serif\" font-size=\"14.00\">Bernoulli</text>\n",
"</g>\n",
"<!-- theta&#45;&gt;y_obs -->\n",
"<g id=\"edge3\" class=\"edge\">\n",
"<title>theta&#45;&gt;y_obs</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M94.09,-166.44C103.31,-154.45 114.54,-140.98 126,-129.95 134.12,-122.14 143.43,-114.49 152.57,-107.55\"/>\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"154.94,-110.16 160.89,-101.39 150.77,-104.53 154.94,-110.16\"/>\n",
"</g>\n",
"<!-- b -->\n",
"<g id=\"node2\" class=\"node\">\n",
"<title>b</title>\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"315\" cy=\"-198.43\" rx=\"41.94\" ry=\"37.45\"/>\n",
"<text text-anchor=\"middle\" x=\"315\" y=\"-209.73\" font-family=\"Times,serif\" font-size=\"14.00\">b</text>\n",
"<text text-anchor=\"middle\" x=\"315\" y=\"-194.73\" font-family=\"Times,serif\" font-size=\"14.00\">~</text>\n",
"<text text-anchor=\"middle\" x=\"315\" y=\"-179.73\" font-family=\"Times,serif\" font-size=\"14.00\">Normal</text>\n",
"</g>\n",
"<!-- b&#45;&gt;y_obs -->\n",
"<g id=\"edge1\" class=\"edge\">\n",
"<title>b&#45;&gt;y_obs</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M294.28,-165.76C285.86,-153.95 275.62,-140.78 265,-129.95 257.56,-122.37 249,-114.94 240.51,-108.17\"/>\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"242.48,-105.27 232.44,-101.9 238.19,-110.8 242.48,-105.27\"/>\n",
"</g>\n",
"<!-- a -->\n",
"<g id=\"node3\" class=\"node\">\n",
"<title>a</title>\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"197\" cy=\"-198.43\" rx=\"58.88\" ry=\"37.45\"/>\n",
"<text text-anchor=\"middle\" x=\"197\" y=\"-209.73\" font-family=\"Times,serif\" font-size=\"14.00\">a</text>\n",
"<text text-anchor=\"middle\" x=\"197\" y=\"-194.73\" font-family=\"Times,serif\" font-size=\"14.00\">~</text>\n",
"<text text-anchor=\"middle\" x=\"197\" y=\"-179.73\" font-family=\"Times,serif\" font-size=\"14.00\">HalfNormal</text>\n",
"</g>\n",
"<!-- a&#45;&gt;y_obs -->\n",
"<g id=\"edge4\" class=\"edge\">\n",
"<title>a&#45;&gt;y_obs</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M197,-160.79C197,-149.38 197,-136.65 197,-124.63\"/>\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"200.5,-124.31 197,-114.31 193.5,-124.31 200.5,-124.31\"/>\n",
"</g>\n",
"<!-- y_obs&#45;&gt;y_obs -->\n",
"<g id=\"edge2\" class=\"edge\">\n",
"<title>y_obs&#45;&gt;y_obs</title>\n",
"<path fill=\"none\" stroke=\"black\" d=\"M235.69,-98.89C250.65,-99.88 263.08,-92.41 263.08,-76.48 263.08,-64.4 255.94,-57.18 245.97,-54.83\"/>\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"245.92,-51.32 235.69,-54.06 245.4,-58.3 245.92,-51.32\"/>\n",
"</g>\n",
"</g>\n",
"</svg>\n"
],
"text/plain": [
"<graphviz.graphs.Digraph at 0x29c407df0>"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pm.model_to_graphviz(model)"
]
},
{
"cell_type": "markdown",
"id": "d03f4721-34ad-4406-8524-45d5bd2923e6",
"metadata": {},
"source": [
"## サンプリング\n",
"ここまでくれば，あとはMCMCでサンプリングするだけ．以下のコードでMCMCを実行．"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "c8597f4a-8b01-4f2a-8d37-523ad383fc7c",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [theta, a, b]\n"
]
},
{
"data": {
"text/html": [
"\n",
"<style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
"</style>\n"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [16000/16000 04:01<00:00 Sampling 4 chains, 0 divergences]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 248 seconds.\n"
]
}
],
"source": [
"with model:\n",
" # デフォルトだとサンプリング数が1000なので，多めに\n",
" idata = pm.sample(draws=3000)"
]
},
{
"cell_type": "markdown",
"id": "bef65c66-bef8-44d0-bb2a-d00ba30030ad",
"metadata": {},
"source": [
"事後分布のフォレストプロットを表示．回答者は2000名もいるので，設問のパラメータのみ表示する．"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "6d181291-ce86-4836-8fbe-2801f54de350",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([<AxesSubplot:title={'center':'95.0% HDI'}>,\n",
" <AxesSubplot:title={'center':'r_hat'}>], dtype=object)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAIeCAYAAAC86qZ+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABY7klEQVR4nO3de5wcVZ3+8c9DFAIYjVwCIUAAATFhFRRcQFQuwqqAhuW2EF3A34qKoEaBVWEhQlQUIUIQRXSJF9jlIgRBEJCbclNBkJUQAYGYAAkEiAnkRpLv749THXo6PZmenu706e7n/XrNa9JVp6q/E4acrqpznqOIwMzMzDrXGq0uwMzMzJrLnb2ZmVmHc2dvZmbW4dzZm5mZdTh39mZmZh3Onb2ZmVmHc2dvWZK0p6Q7JS2U9KKkn0naqKLNFpKil6+hNbzHGpK+IukpSYsk/VnSQVXafUnSTEmzJX1T0hoV+/9Z0nxJW9T4s40vanxdlX1bF/uOKtt2VMXP9kpR89WSDpWkKucJSRNqqcesm0iaLGlmA883RtIXG3W+ZlnpHxuzVpP0XuAm4EbgIGB9YAJwi6R3RcTiikO+CfyyYtv8Gt7qDOAE4GTgfuDfgCsk7R8R1xe17AWcCXy2OOeFwF+BycX+QcD3gW9ExFP9+kH77xBgJrAWsDmwH/A/wDGSDoiIhU1+fzNb2RjgA8A5La5jldzZW45OA6YDYyJiKYCkR4A/Av8PuKCi/RMRcW9/3kDSMFJHf2ZEfKfYfJukrUmd+/XFtg8BN0fED4vj3ldsm1zs/wywNlA6RzM9GBGPl73+maQrgCuAbwPHr4YazLIlaa0qFwOGb+NbnnYhdbBLSxsi4j7gBeDABr3HvwBrAj+v2P5z4J8kbVm8XhMov2JeAAwGKB4rnAEcGxGvNqiufomIXwDXAJ+UtE4rajBrhbLHYdtLulHSy8Dl/Th+R0m/k7RA0mOSPl2xf0NJF0p6tGgzQ9KlkkaUtZkMHAmMKHvM9lSDfsSGcmdvOVoGLKmyfTGwfZXt35S0VNI/JP1S0j/V8B6ji/M9XrH94eL7qOL774EPSHpncdV/CFC6i/Ad4FcRcVsN71fNIEmvK/8CBtVxnutJt/Z3qrMOs3Z2DXAH8BFgYo3HvBG4lPTh/qOku4bfl7RnWZv1gEXAV4APAicC2wB3SRpctDmD9P/f88CuxVejLkgayrfxLUd/JV3dryBpJDAcKL+CXkx6hn4T6X+27YCvAndLendEPLKK91gPmBsrLw7xYtl+gMtIz+TuL17fBpwn6f3A/sV71mvRAI4t9/fi+/AGnc+snZwXEef285ghpDtytwFI+i3pbt/hpP/HiYi/Ap8vHVCMz7mL9P/bh4CrI+Jvkp4HlvT3UeLq5it7y9G5wLslTZA0TNJ2wM+A5cUXABHxbER8OiKuiojfRcRFwPuAIA26G7CIWBYRhwIjgJERsRfprsP3gFMiYrakz0t6QtJzkn4gae0aT78LsHPFVz1XBaXR+F7VyrrR1XUcs6D8jlzxnP9R0sDXFSR9ppil8zKwlNc+WL+13mJbxVf2lp2IuKTo4Esj5YN0hX091W/jlx87Q9KdpI5zVV4ChkpSxdV96Yr+xfLGEfFM2csvkK7Kvy9pH9KtvPcBT5NmEHwV+K8+3h/g/vJxCQCS5tZwXKXNiu/P1nGsWbur5/f+pSrbFlOMxwGQdDxwHmmU/YnFMWuQHuMNrnJ81nxlb1mKiP8CNgDeDgyPiMNJz8vurPUUfex/mPSc+y0V20vP6qdWO0jSpsApwGciYjnpWd7NEfFgRDwPXFxsW532I334uL+vhmYdqFl3tP4NuCUivhQRN0XEH4HnmvReTefO3rIVEa9ExP8Vt8o/SHo+/oNVHSNpc2B34A99nP7XpOf/Yyu2fwz4S0Q82ctx3wUuKf7HL1m37M9v4LXb6k1XhAB9BPhBRCxYXe9r1gXWoecYIYCjq7RbTJp+mzXfxrfsSNqRNADmT8Wm3Um30b4dEXeXtTub9IH1HtIAvbeSRs4uB75ecc6lwE8i4v8BRMRzks4BviJpfvFehwF7kTrPanX9C/Beej6v+w3weUnHAs+Q5rpPrvdn78MOkjYgTQfcnDRA8BDgZtLPbWaN82vgPyV9lXTxsBdwcJV2U4H1JH0GuA9YFBH/t/rKrI07e8vREuDDwEmkW+2PAJ+OiIsr2j1MCrU5inRF/QJwK/C1YiRtuUGsPK3tZOBl0ojbjUmzAA6NiOsqC5K0FnA+cGJEzC1tj4gbin8Mvkq6EphCSvtrhiuK74tItxP/RLrVeGWVWQVmNjCnA0OBcaRn9HeQRuw/UdHuR6TBtt8o2k8HtlhNNdZM/jfCzMyss/mZvZmZWYfzbXwzM2t7xeqPq0ygrJzq2k18ZW9mZp3gSNLo+VV9dS0/szczs7YnaX1gy1W1KRbU6kru7M3MzDpcxz6z/+AHPxhz5sxpdRlmZnW5//77b4yI1Z3GWDP/G5un3n5vOvnKvmN/MDPrCqstibFO/jc2T1V/bzxAz8zMrMO5szczM+tw7uzNzMw6nDt7MzOzDufO3szMrMO5szczM+twA+7sJW0hKYqvaY0oqpf3GV/2Pic0633MzMw6TSOv7D8I7F6+QdJgSf8l6RFJiyS9KOk6Sf9c0W64pEslTZO0TNLkKuf/DjAcmNnAms06xnPzF3H/9Jd4bv6iVpdiXe6SSy5ly222Y41Bg9hym+245JJLW11S12tkgt4LEbEiTknSmsBNwFuAk4A7gfWA44HfSTooIq4tmq8FzAHOBI6pdvKIeBl4WdKyBtZs1hGmPPA04699mJHrrcP0Fxcw/oDRjNlxRKvLsi50ySWX8plxJ7LO3sex2ZhRLJw5lc+MOxGAsWOPaHF13aumBD1JHwROBrYnpSb9EfhCRDwiaQvgSWDn8kUGJJ1E6rx3iog/VZxvCrAbsEVELKjYdx0wJyKO6qWWp4DzI+I7fZTtdCfraIddeA8AS5YuZ9qs+fzyuPewzUZDeGz2fD76vbt460ZDWPN1a3DZp3ZtcaVWp7ZM0Ntym+1YuNORDB759hXbFk1/iLXv+wlPPta0J732mgEl6K0LfBd4N7AH8A/g2uLqvTdjgd9UdvSFs4ANgX1qfH8z68XipcvYZOhgttloCADbbDSETYauzeKlvglmq9/0Jx5jrU1H9di21qajmP7EYy2qyKDG2/gR8Yvy15KOBuaROv/enqFvC9zey76pxfe31vL+Zray0hX7c/MXsdd37uCx2fNXXNnPeXkxN417H8OGDG5xldZtRm61DQtnTu1xZb945lRGbrVNC6uymjp7SW8BzgD+mXRFvkbxtTkDGzC3ZADHmhkwbMhgJozZnkMuvKfHM3t39NYKE8afmp7R730ca206isUzp7LglvM5e+JZrS6tq9U6QO86Uqf+KeBpYCnp6nxVt/EfBUb1sm9UWRszG6AxO45gt63XZ8aLC9lsvbXd0VvLlAbhnTL+dKZf/hgjt9qGsyee5cF5LdZnZy9pfWA74NiIuK3Y9s4ajr0EOFPSO6s8tz8JeAa4uf8lm1k1w4YMdidvWRg79gh37pmp5cr+JdK0uE9KmgGMIA2wW9rHcd8FDgB+WYzMvwt4M/A50pz8f4mIV0uNJe1Q/PGNwPLi9ZKImIqZmZnVrc/OPiKWSzoMOA/4C/A48CXgF30ct0TSPsB/AqcBWwKvJ31w2CEiHqk45IGK1wcA04Et+v4xzMzMrDc1Tb2LiFsjYvuIGFx8vzEi3hARk/s4blFEfC0i3hoRa5Lm1q8FHFWlrap8bdH/H8nMzMzKNTIu949Fbv0uvTWIiHuADwALJW1V64nL8/eBkQ2o1czMrGs0orOfCfwBuBLYCbi/tEPSv0i6RdI8SQsl/RnYFTg9Ip4o2mwh6ceSnijaPCHpm5LWLnuPGcDbgB8DjwE/akDdZl3JGfrWbM7Gz8+As/EjYqmkhcDsiCjv6I8FJgFnA58HXgH2Bb4F7AIcXjTdDhgEfIbUkb8N+CGwPkVOfkQsA6ZJmkkatDd3oHWbdSNn6FuzORs/TzVl4/d5Eul24C8RcVzxelPgb8D3I+ILFW3HAFcDh0bEFb2c71jgjIhYv2L7eODgiNi+hrKcjW9dr5SfD9Uz9D9y/l1st3HK0C9xln42nI1v9RhQNn5/HUIK3Pl25Y6ImEK6gl/VR7w3kqb8mVmDVM/QH+wMfWsoZ+PnqZFL3JbbFpgXEc/0sv8ResnFlzQSOAH4RpNqM+sa5Vfp1TL0Z89bzK0nvN9hPNYwzsbPU7M6e+j7NvpKufiSNgJ+TUrWm9iMosy6VbUM/QljtndHbw3lbPw8NauzfxR4k6QREfF0lf2jqAjRkbQxcCspuOfj0YjBBGbWgzP0rdmcjZ+nVgzQOxC4CtgvIq4vtg0HbgMeBg6LiKpRvB6gZ2ZdpC0H6FnLVf29acqVfUTMlDQOmCRpCfATYAGwD2nQ3kVlHf0mpHXvnwG+AGwgraj1+WLanZmZmdWpac/sI+ICSU8CJ5Lm0L+h2HVCRJxd1nRfYJvi6+8Vp9kSeKpZNZqZmXWDRt3Gnwu8qXi5a0TcW6XNusCvgA2APSPi+X6cfwvgyeLl4oio5UGjbzGZWTvzbXyrR9Pn2S8n3ap/7R3L4nJJq92tT7pSf19ZmzUk/VLS3yUtkvSspJ9Lqoz1WlC8h5mZmfVDozr7qcClwDsoRtkXKXjXk7LydyONwD8f2B04uOL4W4FDSXPvDwK2IqXslcwszv09fFvfbJWcfW+t5mz8/DTqmf0S4B8R8TisGI0/EZgUESeVtbtQ0mzgaklXRcQVEbEc+G5Zm+mSzgSukTS4WCZ3KfC4pBeBqiP1zczZ99Z6zsbPU7Om3o0DzgFGVEvRk/Qo8HBEHFhl33rA94GREbFLxb7xeOqdWQ+l/Ptasu+de99W2vKZvbPxW261ZuP3Oy5X0rckvQK8AGwO7N+k2sw6krPvLQfOxs9TTnG5Z5HWqx8JnAb8XNKHnKRntmqlq3Vn31sOnI2fp2ziciNiDmnE/qOSHgFmkAbz/a5JNZp1FGffWw6cjZ+nZnX2VwBnkgJ1vlC+o4jL3Rr4/CqOLz1eWKsZxZl1KmffW6s5Gz9PTRmgV2w7FpgEnM3KcbmXR8QxRbtdgXcCdwJzgbcAZwCbAG+LiEVl5xyPB+iZWXdoywF61nKrLxsf+hWXu5A07/50YF3gWdIyt4eVd/RmZmZWn2YO0CMibgBugB5xuUdL+mkpLjciHgT2bGYdZmZm3awZ2fifiIiLe2m3LulZ/bSIuKof598NuKt46Wx8M+sGvo1v9WjqPPupwJXArsDlK96xLBtf0kLgbuAVYErVCqXBkv4sKSTtVLbrvuLcP8ZxuWZmZv3SqM5+CTA7Iu6NiFeg12z8C4CvAZf0cp7vkHLwe4iIJcVKejNxXK5ZVc7Et1w4Gz8/TXlm359s/LJjPkp6dn8w8OFm1GXWqZyJb7lwNn6essjGLz4c/AH4EPAP0tr1O0fEfRXHjcdT78xWOOzCe1iydDl/nT2faz5bPRPfefhtqy2f2Tsbv+XyzMaXNIh0W//siPhzk+ox61iLly5jxNC1nYlvWXA2fp5yyMb/avHnc5pYi1lHuuxTu/Lc/EXsO/G3zsS3LDgbP085ZOPvDbwXeFXqcffhXkmXRcTYJtVo1hGGDRnM+ANGOxPfsuBs/DzlkI1/NCk5r2QT4EZgLK/NrTezVXAmvuXC2fh5akpnHxEzi0F6kyQtYeVs/Isi4vqi7ZPlx0p6ufjj3yJipWl4ZlbdsCGD3clbFsaOPcKde2ZyyMY3MzOzJmrWaHwgZeNHxF4RMYTU2d9BysbfcBXHPBURqpx2Z2ZmZvVxNr6ZWZ7acp69tVy+2fiSniry8Mu/zixr4mx8MzOzOuWUjX86MLzsa0Jph7PxzernzHxb3ZyNn59ssvGB+RExqxn1mHUrZ+bb6uZs/Dzlko3/FDAYeD0wgzRP/6yIWFJx3HicjW/Wq8MuvGfFn5csXc60WfP55XGvZeYfMOlOhq7zekau/1q0hbPzs9WWz+ydjd9yeWbjF84DDietenc+MI50y9/M6rR46TI2GTq4R2b+8KFrs3S5Pwdb8zgbP085ZOMTEeW5+A9JmgdcJuk/I+KFplRn1oHKr9Kfm7+Ivb5zR4/M/GfnLuRju4zklP1HreIsZvVzNn6ecsjGr+b3xfetAXf2ZnUYNmQwE8Zs3yMz/8yD3u5n9tZUzsbPUw7Z+NXsUHx/tgm1mXUNZ+bb6uZs/Dw1ZYBese1YYBJwNitn418eEccU7XYFdgFuA/4B7EwayX9fRHy04n3G4wF6ZtYd2nKAnrVc1d+bHLLxFwOHAacBawHTgYtIHwrMzMxsgJo5QI+IuAG4AVZE5f6KlI3/04h4vmjzJ9KVvZmZmTWBs/HNzPLk2/hWj3yz8cva3iNpgaS5km4t2+1sfDMzszplkY0vaQzwv8DPgB15rWMHnI1vVg9n4lurOBs/Py3Pxpc0iJSgd1JEXFTW9pFm1GbWDZyJb63ibPw8tTwbX9K7SSE6R5Ge528CPAT8Z0Q8UHHceDz1zmwlfWXif+T8u9hu4yGs+bp0M895+G2hLZ/ZOxu/5bLNxt+q+H468A1gP9Lt+tslDW9SfWYdq1om/iZDB7N46bIWV2bdwNn4ecohG7/0gePrEXElgKRjgA8A/w58qznlmXWOvjLxZ89bzK0nvN8JetZ0zsbPUw7Z+KVI3KmlnRGxVNJjwOZNqs+sY1XLxJ8wZnt39LZaOBs/Tzlk499PStF7K3Bn0WYN4C3AjU2qz6yjORPfWsXZ+HlqSmcfETOLQXqTJC1h5Wz8iyLi+qLtPEk/AL4maSZpHv1xwJtJU/HMrA7Dhgx2J28tMXbsEe7cM5NDNj5Fm9KHgnWAPwF7RoRXvTMzMxugZo3GB1I2fkTsFRFDSJ39HaRs/A0r2r0aESdFxMYR8caI2KPIzDczM7MBapds/M1Jq+EBLI+IQTUc5nn2ZtbO2nKevbVcntn4kvaQFL18HVI0e6Y49wTA2Z9mZmb9kEM2/t3A8IqvbwIvUyyPGxFLi2z8v+FPk2Y9OAPfcuNs/Py0PBs/IpYAsyqOPxj4n4h4uRn1mXUKZ+BbbpyNn6eWZ+NX2bcHcBuwc0TcV7HvKOD8iHhD5XFV+A6AdaRSDv6Spcv56+z5XPNZZ+B3qLZ8Zu9s/JbLNhu/0jHAg5UdvZn1tHjpMkYMXdsZ+JYVZ+PnKYds/BUkrQ/8K/DFplRk1gFKV+vPzV/EvhN/6wx8y4qz8fOUQzZ+uX8HltFzAJ+ZVTFsyGDGHzDaGfiWFWfj5ymHbPxy/wFcERH/aFJdZh3FGfiWG2fj56kpA/SKbccCk4CzWTkb//KIOKbiHLsDvwN2j4i7enmfo/AAPTPrDm05QM9arurvTS7Z+ACfBB7praM3MzOz+mSRjV+0PTIiRq10EjMzMxuQdsnG3wJ4snjpbHwz6wa+jW/1yDMbv2i3raQpkuZImi/pXkkfLGsyozj3j4HHG1SzWdtzVK7lyHG5+WnUM/sV2filDRUD9D5P6uT3Bb4F7AIcXnb8dcATwN5Fu08D10gaFRF/i4hlQOkDwKsNqtmsrTkq13LkuNw8NSsud1PSojXfj4gvVLQdA1wNHBoRV0jaAHge2CsibivavA5YDBwWEVeWHTseODgitq+hLN9iso5SismFFJU7bdZ8fnlc9ahcx+R2hLa8je+43JZbrXG5hwBrkqbZ9RARU4DHgNJHvBdI8bkfl/QGSYNIkbnzAY/MN6ti8dJlbDJ0sKNyLTuOy81Ts6be1ZyNHxEhaR/S1f48YDnwIvChiHi2SfWZtZ3yq/Xn5i9ir+/c4ahcy47jcvPU8mx8SSKtc/8C8F5gISlJ7xeSdu4lbtesqw0bMpgJY7Z3VK5lx3G5ecohG38v4ABgvYiYW2w7trjaPxqY0KQazdqao3ItR47LzVMO2fjrFN+XV5xjOU0O/TFrd8OGDHYnb9kZO/YId+6ZaUpnHxEzJY0DJklawsrZ+BdFxPVF83tIz+gvlnQ66Tb+J4GtSFPyzMzMbABano0fEXOK+fNfB24FXk8awDcmIv7UrPrMzMy6RTMH6BERNwA3wIqo3F+RsvF/GhHPl7W7D/iXZtZiZmbWrdolG383XptzvzgianlI6VAdM2tnbRmqYy2XdTb+OyXdLGmupBck/VBS+Zr19/FaNv5TDarZrK04B9/ahbPx89Oozn5FNn5EvAIrsvGvB+4HdiNNt7sA+BpwSelASZsAvyFl4/8z8EFgNDC51CYilhS5+zOBpQ2q2axtTHngafad+FtOv/Zh9p34W6Y84PgJy1MpG3/hTkey2RevYuFOR/KZcSe6w2+xHLLxjwG+AWxULHiDpH8CHgK2iYjHy44dj7PxrUuUsvCdg9+12vI2vrPxWy7bbPy1gFdLHX1hYfF99ybVZ9Y2nINv7cTZ+HlqeTY+abrdOZK+DJwDrEsK5AEY3qT6zLJXumJ3Dr61E2fj56mZCXU1ZeNHxMPAkaSkvYXALOBJYDYrp+qZdZ3yHPyPnn8nh1x4j3PwLVsTxp/KglvOZ9H0h4hlS1k0/SEW3HI+E8af2urSuloO2fhExKXApZI2Io3WD+CLpEF7Zl3POfjWLpyNn6dWDNA7ELgK2K8sMrfyfJ8AJgEjyhbH8QA9M+smbTlAz1qu6u9NDtn4SDqOlJE/v2hzFvDl8o7ezMzM6tPybPzCu0nz798ATAM+FRE/a1ZtZmZm3aSm2/jFbfr3Fy93LQJuKvevuI3fyzlK2fgbAHuWZ+PXVKh0FFCK4f3eqt6r4FtMZtbOfBvf6jHgefYXk6bC3b/ijEUcLmk+/GclLZN0tqQe55V0MnATsDMpHe+9K1UnbS7pWkmvSJoj6TxJa5b2AecDr/LaHHwzMzOrQX86+wURMSsiXoWV4nD3BfYETgP+H2VxuIW1SIPyvgtQuQiOpEGkq/4hpA8ChwMHA6Xb/c8AO5BG8f+5HzWbdR1n6FurORs/P3U9sy9G208EJkXESWW7bpf0F+BqSVdFxBUAEXFqcdzBvZxyX9IV/8iImFG0PQn4kaSTI2Ie8HixfXE9NZt1gykPPM34ax9m5HrrMP3FBYw/YDRjdhzR6rKsi5Sy8dfZ+zg2GzOKhTOn8plxJwJ4+l0L9eeZffnUunGktLsR1VLyJD0KPBwRB1ZsPxi4IiJUsf104KCIGF22bUPgOWCviLitt1pWwc+TrCvUmqEPOEe/vbTlM3tn47dcQ7Px+xOHW4uNSYl55eYAy4p9ZtYHZ+hbDpyNn6eBTL2rKQ7XzJrLGfqWE2fj56nezr5fcbg1mAW8p2LbBsCgYp+Z9aE8Q7/0zN4Z+ra6TRh/anpGv/dxrLXpKBbPnMqCW87n7Ilntbq0rlZvZ38FaWW6E0kL2KxQxOFuDXy+H+e7BzhF0qYRMbPYtg+wmLKpfma2as7Qt1ZzNn6e6hqgV2w7lpRffzYrx+FeHhHHlLXdHFgP+AApCnfHYtfjEfFyMfXuQeB54EvA+sU5r4qI4/uqpRceoGdm7awtB+hZyzU2G7+fcbink5axLSnd4t8TuD0ilknaD7gAuIsUnHNJcW4zMzMbgLqv7Ku0GVAcbq18ZW9mXcJX9laPAU2924EUhxuSjq7WICJeAfYDLqVKHO5ASZogKUgZ/fs1+vxmZmadqtbOfipwJbArcHlpYykbX9I8SQuBu4FXgCnlB0s6WdJdRe591U+Dks6VdJ+kRZKeqtLknOL9ryfNBjAzM7Ma1NrZLwFmR8S9xRV8ZTb+bqTpdheQlqpdZTb+Kmr5CfDTajsj4sVitb3nSQvimBnOwrf8OBs/P43Oxr9Q0mz6n41PadS9pBNIWflm1gdn4VtunI2fpyyy8SvanAAcFxFb9LJ/MrBBROzfR9kePGIdp5SDD9Wz8D/6vbt460YpC985+G2vLQfoORu/5bLOxjezfqqehb+2s/CtpZyNnydn45u1kfKr9WpZ+HNeXsxN497n5DxrGWfj5ymXbHwz66dqWfjjDxjtjt5aytn4ecolG9/M6uAsfMuNs/HzVFdnHxEzi0F6kyQtYeVs/Isi4vpS+7Js/C2K1zsUux6PiJeLbVuTInc3AdYsazM1IvxIwKwXw4YMdidvWRk79gh37pnJIhu/+POPSOl4lW22BJ6qt04zM7Nu147Z+JPx1Dsz63xtOfXOWq7q702tnf1c4E3Fy09ExMW9tFuX9Kx+WkRcVV+dvdZwPHBe8fK5iNioj0P8i2hm7cydvdVjQPPsm5qNL+kdkv5H0gxJCyX9VdJJksrr+wmvZeM/VGPdZmZmXS+XbPx3kTLvPw6MBk4D/gv4cqlBRMwry8ZfXGPdZtYLZ+pbszgbPz9ZZONHxH9XbHpC0juBg4Bv1FOjmfXOmfrWLM7Gz1N22fhlbX8AbB0RH6jYPhkP0DOrWXmePsD0F15h7oJXufb43Vck733k/LvYbuOUqV/ibP2Wa8tn9s7Gb7n2ycYvruqPAr5f7znMrLqly4PhQ9euyNQf7Ex9awhn4+cpu2x8SW8lTeH7bkT8op5zmNlrKq/QJ1w3lZ/fO71Hpv7seYu59YT3O5zHBszZ+HnKKhtf0nbAbcD/RsSX+2pvZv13yv6j2H7Em3pk6k8Ys707emsIZ+PnKZtsfEmjgFuByyNiXJ11mVkNnKlvzeJs/DzVnaBXTL2bBJzNytn4l0fEMWVtS9n4HwDOAnYsdj0eES9LGk3q6G+j4sNDRMyqqGUyHqBnZp2vLQfoWctV/b3JJRv/EGAYcFjx1WfhZmZmVhtn45uZ5Sn3Cx3/G5snZ+ObmbURd/ZWj6yz8TeUdKOkZyQtLjLyvyfpTWXNnI1vZmZWh1yy8ZcDVwMHkAJ7jgL2Bi4qNXA2vll1zri33DgbPz+5ZOO/APygbNN0SRcAX6mnPrNu4Yx7y42z8fOUZTa+pE2AS4FZEfFvFfsm4wF61qXKc+6XLF3OtFnz+eVx76mace9s+7bXls/snY3fcvln4xdr2i8AngbmA0fXWZ9Zx1u8dBmbDB3sjHvLirPx85RbNv440jP/bYFvkp7xf6qO85h1pPKr9efmL2Kv79zhjHvLirPx85RVNn6RljcLmCbpReB3kiZExIw66zTrWMOGDGbCmO2dcW9ZcTZ+nrLJxq+i9IhhrQGex6xjOePecuNs/DzV1dlHxMxikN4kSUtYORv/ooi4vtS+LBt/i+L1DsWuUjb+/sD6pGl8LwOjSRn690bE4/XUaNYthg0Z7E7esjJ27BHu3DOTSzb+IuDTwNtIV/IzSPPuz6y3PjMzM0ucjW9mlqe2nHpnLedsfDOzNuLO3uqRbzZ+RdsNJD0tKSRtULbL2fhmZmZ1yCUbv9zFwIOVG52Nb9Yczta3RnM2fn6yyMYvO+/ngXWArwMfrqc2M6uds/Wt0ZyNn6dssvEl7Uga4LczsA1wG7BhRMypaDcZD9Az67fyXH2onq1/wKQ7GbrO6xm5/ro92jpnvyXa8pm9s/FbLt9s/GJg3/8Cx/eSyGdmDVYtW3/40LVZutyfk61+zsbPUy7Z+OcBd0bELwZQj5mtQuXVebVs/WfnLuRju4zklP1H9XIWs1VzNn6ecsnG3xvYTFIpeKd0G2KWpG9FxMl11mlmvaiWrX/mQW/3M3sbEGfj5ymXbPx9gTXLXu8M/DewB+B7P2ZN4mx9azRn4+cpi2z8iHi0/Pxl8+unVQ7QM7PGcra+NZqz8fOTSza+mZmZNYmz8c3M8tSWU++s5ZyNb2bWRtzZWz3yzsYvsvArvz5d1sTZ+Nb1HG1r7cBxufmp9Zn9imz80oYiG38ScDbpav4V0qj6bwG7AIeXHV/Kxr8d+Ooq3ueTwHVlr/9R+kNEzAPulfQ86VGBWVdxtK21A8fl5qneuNxNgb8B34+IL1S0HQNcDRxaysYv27equNwADomIK/uoZTJ+Zm9d5LAL72HJ0uX8dfZ8rvnsa9G2Hzn/LrbbeAhrvm4Nx9l2pra8je+43JZraFzuIaR58d+u3BERU0hz4+v5CHeupDmS/ijp05Lqrc+soyxeuowRQ9fuEW27ydDBLF66rMWVmfXkuNw81Tv1rqHZ+IVTSYvfvExK1DubdLt+Qp01mnWEyz61K8/NX8S+E3/bI9p29rzF3HrC+z1H3rLiuNw85ZKNT0ScUfbyQUmDgJNxZ2/GsCGDGX/A6B7RthPGbO+O3rLjuNw85ZKNX83vgTdK2igiZg/wXGZtz9G21g4cl5unXLLxq9kBWATMHeB5zDqGo22tHTguNz9ZZONLOgDYGLgHWEiK0T0d+GFELK6nRjMzM0tyycZ/FTgWOIc0Q+AJ0oC979Vbn5mZmSXOxjczy1NbzrO3lnM2vplZG3Fnb/XIOxu/aPcxSQ9KWlSE6/y0bLez8c0qOCvfcuRs/Pxkk40v6XPAV0hjAO4F1iaF9wDOxjer5Kx8y5Gz8fOURTa+pKHA08CYiLi5j1om42f21qUOu/AeAJYsXc60WfP55XHOyu9gbXkb39n4LZd1Nv6+wCBgI0lTJT0t6WpJW9VZn1lHW7x0GZsMHeysfMuOs/HzlEs2/lakDx6nkEJ6XqTIypf0tohYUGedZh2ldMX+3PxF7PWdO5yVb9lxNn6ecsnGXwN4PfC5iLgJQNJYYBZwAHBZXRWadahhQwYzYcz2zsq37DgbP0+5ZOM/W3yfWtoQEf+Q9AyweZ01mnU0Z+VbjpyNn6dcsvHvKr6/FZhZnOcNwHBgep01mnU8Z+VbjpyNn58ssvEj4lFJ1wDnSvoU8BLwNeA54Lp6ajQzM7Mkl2x8gI+TsvGvJU0duBPY24PzzMzMBqbPefbFHPv3Fy93LQ/WqWjXtGx8SVsATxYv5wJ3eZ69mXW4tpxnby03oHn2s0gr0/1G0s7Qa1Tur4D/Ad674l37iMqVdJSk6OVr56LZDGBMUcPQGms2MzMzau/sf00aYb8D8FARlXs9cD+wW7HvAuBkYPuKRXBKUbnf7eXcl5EG4pV//Zy0zO19ABGxDLgJ+AEpsOeTNdZt1pGciW85czZ+fmp9Zv9KRDwOK6JyJwKTIuKksjYXSpoNXC3pqlJUbkScWhx3cLUTR8RCYGHptaR1SHPrvx1lzxgiYqGkF4ElEfHsymcy6w7OxLecORs/T7U+sy/PxR9HGkg3olqCnqRHgYcj4sCK7VVz8ascfxRwEbBZRMyq2DceODgitl/1jwX4eZJ1kFoz8QHn4neOtnxm72z8lmtYNn6jo3IrHQNcV9nRm5kz8S1/zsbPU71T7xoZlbuCpNGkNev3q+d4s07lTHxrF87Gz1M9nX2jo3LLHUMaef/rOo8362jOxLfcORs/T/V09o2Oyi0dO5gUrHNeRCyvoy6zruBMfMuZs/Hz1O/OvtFRuWWnPhh4E/Df/f8xzLqLM/EtZ87Gz0+92fiNjsqFNHf+xoj4ez01mZmZWXX1jMYHICJuiIi9ImIIqbO/Azha0oYV7Y6KCFX5ur2i3fsj4sP11mNmZmbV1drZH1sWYbtL5c6IeIU0gv5SyqJyG0XS5kUk72nAaEleCc/MzKxGtXT2Y4E/AlcCI0kRuUDPfHxgDnAYsJmkHuftKx+/aLOzpN9Imlt83SLp3cXuZ4B3kqblXdvfH9LMzKyb9dnZF9PrFgKzI+LvEfEqwCry8b8GXFJxmlXm40t6A2m63TPALqRO/VngRklDImJpRDxSrLj3Yn9/SLNu5Px8axVn4+enrgF6jc7HB7Yjjdg/LSKeLNr+F+muwlspFsQxs9o4P99axdn4eeozGx+an48vaQjwN+CHwBnF5tNInf12xWI5pbaTgQ28nr3Za0rZ+eD8/A7ibHyrR8Oy8aHB+fgRMR/YAziUNGd/Aen5/z7lHb2Z9c35+dZKzsbPU73Z+NDAfHxJa5PCdO4lXc0PAk4ArpG0UzHa38x6UX6l7vx8ayVn4+ep3s6+0fn4RwBvAd4TEcsAJB0BvAQcCPy8zjrNuo7z862VnI2fp3o7+0bn469DulNQnom/vNhWd/CPWbdyfr61irPx81RvXG6j8/FvBs4CLpB0HqmD/zKwDLi1nhrNup3z861VnI2fn7qf2TcyHz8ipkk6gDQC/x7SFf2DwIciYma9NZqZmVmdU+96abMu8CtgA2DPiHi+UUVWvM9kPPXOzDpfW069s5ar+ntTa2c/l7T8LMAnIuLiXtqtS3pWPy0irqqvzl5rOB44r3j5XERs1Mch/kU0s3bmzt7qMaB59lNJ2fi7ApevOGNZNn6xUM3dwCvAlB7vXFs2/t6S7pY0X9IsSd+SVP6Y4SfF+18PPFRj3WZmZl2v1s5+CSkb/97SnPcGZ+O/ozjXzcCOpECdj5BG/AMQEfOKbPzngcU11m3WVZyHbzlwNn5+csnGPwyYGhGnFa8fl3QScLmkrxUJe2a2Cs7Dtxw4Gz9PuWTjnw3sFhG7lm37AOlKf8+IuL1s+2Q8QM8MeC0Tv688fGfht6W2fGbvbPyWyzcbH7gR+GdJH5P0OkkjgFOLfcPrrNGsazgP33LhbPw8ZZGNHxE3SToB+B4wmfRM/gzgvfRM1TOzMqUrdufhWy6cjZ+neq/sV2Tj97J/VNGmZhFxDjAU2Jw0V/+aYtcTddZo1jXK8/A/ev6dHHLhPc7Dt5aYMP5UFtxyPoumP0QsW8qi6Q+x4JbzmTD+1L4PtqbJJRsfgEgDCJ4pznM4MAP4U501mnUV5+FbDpyNn6dcsvGRdCLwa9Jt+38lZeMfWloFz8z65jx8y4Gz8fOTRTZ+8ecPASeT5uT/GfhoRNxQb31mZmaWOBvfzCxPbTn1zlpuQFPvdgA+KykkHV2tQZGstx9wKWkUfUNJ+kYRtXskMLrR5zczM+tUTc/Gl7SFpB9LekLSwuL7NyWtXf4GkjaXdG2Rnz9H0nmS1ixrclbx/g8At9Xzw5qZmXWj1ZGNvx0wiPRcfzRwPPDvwLmlBpIGkR4BDCHdFTgcOBhY8ew/Il4qsvHnkQYDmhnOw7f8OBs/P03Pxo+IX5NG2Zc8IenrpNCcY4pt+5I+CIyMiBnFe5wE/EjSyRExr546zTqd8/AtN87Gz9NqzcYv2/9l4D8iYuvi9enAQRExuqzNhsBzwF4RcVvZ9h61rIIHj1jHOuzCe1iydDl/nT2faz7rPPwO1ZYD9JyN33J5ZONLGgmcQLrlX7IxMLui6RxgWbHPzCosXrqMEUPXdh6+ZcXZ+Hlardn4kjYi3dK/mfQYwMzqcNmnduW5+YvYd+JvnYdvWXE2fp7q7exXZONHxNNV9o/iteAcACRtDNwK/AX4ePR8fjALeE/FOTYgDeybVWeNZh1t2JDBjD9gNIdceM+KZ/bOw7dWmzD+1PSMfu/jWGvTUSyeOZUFt5zP2RPPanVpXW21ZONLGk6aLvcwcHhELK043z3AKZI2jYiZxbZ9SKvf3V9njWYdz3n4lhtn4+ep7gS9YurdJNL0uMps/Msj4pii3SakSNxngI8Dr5ad+vmIWFZMvXsQeB74ErB+cc6rIuL4vmrphQfomVk7a8sBetZyVX9vVkc2/r7ANsXX3ytOsyXwVNHh70catHcXsJA0V//EeuszMzOzpB2z8fuspeBPnWbWznxlb/UY8NS7YyS9LGnnajtXQzb+WEkvN+PcZmZmnazWK/sRQCnLfkZELG5qVdVrGAJsVLycGxFz+jjEnzrNrJ35yt7qUfX3pqbOvk117A9mZl3Bnb3Vo6EJemZmZtYm3NmbmZl1OHf2ZmZmHc6dvZmZWYerqbOXdLuk81ex/ylJUXw1ZZU6SXuUvcd1zXgPMzOzTtTIK/vTgeGkNehXkHSEpHuKOfqvSPq9pI9VHizpXEn3SVok6akq57+7OP/lDazZzMys4zWys58fEbMiYnlpg6RvARcD1wDvAnYErgJ+LOnMKrX8BPhptZNHxJKImEWK0rUBeG7+Iu6f/hLPzV/U6lLMrANdcsmlbLnNdqwxaBBbbrMdl1xyaatL6nr9ycZ/naRzgX8vXv8I+M/yzr2cpHcDJwHjIuK7Zbu+JWkxMFHSVRHxB4DSgjeSTiDl6VsTTHngacZf+/CKJVHHHzCaMTuOaHVZZtYhLrnkUj4z7kTW2fs4NhszioUzp6Ylb8Er37VQf7Lx3wVMBr4HvB24CPhaRJxT3HY/PyK+U3bMucAngPUjYknF+dYCXgAuiohxFftOAI6LiC16qWUysEFE7N9H2dkHPhx24T0NP+fMlxbw9Nzer9gHv24Nrj1+d7bZaAiPzZ7PAZPuZNHSqp/XGDF0MJu+eZ26a7nsU7vWfayZtWeozpbbbMfCnY5k8Mi3r9i2aPpDrH3fT3jysWmrrbguNuBQnWeBz0XEtIi4HDgL+OIq2m8LPFHZ0QMUcbt/A97aj/e3Bhg+dG222WgIANtsNIThQ9fu4wgzs9pNf+Ix1tp0VI9ta206iulPPNaiigz6dxv/3uh5G+Ae4AxJbxzA+6/0QaCbNOPKd+LNj3LuLb3/T/Xs3IU8Nnv+iiv7Z+f2PgTi4Hdtxrh9tm14jWbWuUZutQ0LZ07tcWW/eOZURm61TQursrrXs6/Bo8B7Ja1VuXBOcRv/LcCNTXz/rjRun21X2UFPeeBpDrnwnhXP7M886O1+Zm9mDTNh/KnpGf3ex7HWpqNYPHMqC245n7MnntXq0rpafzr7f5aksqv7XYBnImKeVPURwSXA54DPAN+t2HcssA69jLy35hmz4wh223p9Zry4kM3WW5thQwa3uiQz6yClQXinjD+d6Zc/xsittuHsiWd5cF6L9aez3wT4rqQLgH8CTgQm9NY4Iv4g6duk0fdrAVeTBnQcCJwBnBIRfym1l7Q18IbifdaUtEOxa2q15/5Wv2FDBruTN7OmGTv2CHfumelPZ38JMAj4PanT/jEwcVUHRMR/Svo/4DjgNKA0GuzQiLiiovmPgPeXvX6g+L4l8FQ/6jQzM7MyDVnPvtrUu17aDQNuI6XsfTgi+h2Q00lT78zMVqEtp95ZyzV1PfuNgbOK3PoP9dYoIp4D9gRuB2oeii5pi1IuPnDkQIs1MzPrJo3q7B8AriR14LeWNkr6F0m3SJonaaGkPwOHA2dExK1Fmy0k/VjSE0WbJyR9U1L5BPAZpNv5PwYeAz7ZoLrNzMw6XqM6+8XA7Ii4tzTNTtKxwPXA/cBuwCjgAuBrpOf/JduRxgJ8BhgNHE+K5D231CAilkXEU8BMYElEPNugus3aktc3sJw5Gz8/TZlnL2lT0uC9SRFxUtmuCyXNBq4ucvGviIhfA78ua/OEpK+TRuwf04z6zNqZ1zewnDkbP0+NGqB3O/CXiDiueD0OOAcYERHPVGn/KPBwRBzYy/m+DPxHRGxdsX08cHBEbF9DWR48Yh2hfA2FJUuXM23WfH553HtWpCB+5Py72G7jIaz5unSjzmsSdIy2HKDnbPyWa+oAvUrbAvOqdfSFR+glF1/SSOAE0i1/MyuzeOkyNhk6uMf6BpsMHczipctaXJlZ4mz8PDUzLrevK+uVgnIkbUS6pX8zfczhN+sW5Vfqz81fxF7fuaPH+gaz5y3m1hPe76Aky4Kz8fPUrM7+UeBNkkZExNNV9o/itdAcACRtTBrJ/xfg49GI5wtmHWbYkMFMGLN9j/UNJozZ3h29ZcPZ+HlqVmd/BXAmKVL3C+U7JB0IbA18vmzbcFLYzsPA4RGxtEl1mbU9r29gOXM2fp6aMkCv2HYsMAk4G/gJsADYB/g2cHlEHFO024QUsvMM8HHg1bJTPx8Ry8rOOR4P0DOz7tCWA/Ss5ar+3jTtmX1EXCDpSdLV/WdIi9wAnBARZ5c13RfYpvj6e8VpnItvZmY2QM0ajQ9ARNwQEXtFxBBSZ38HcLSkDcvaTI4I9fL1VDPrMzMz6waNuo0/F3hT8fITEXFxL+3WJT2rnxYRV/Xj/LsBdxUvF0dELQ8pfYvJzNqZb+NbPZo6z34qr2XjX77iHSuy8YG7gVeAKWVt1pD0S0l/l7RI0rOSfi6pPBLsvuLcP8a39c3MzPqlUZ39El7Lxn8F+pWND2nK3aGkoJ2DgK2Aq0s7I2JJRNxLysb3SH0zGxCvLdBczsbPTw7Z+MuB75a1mS7pTOAaSYMjwv83mlnDeG2B5nI2fp6yy8aXtB7wfWBkROxSsW88nnpn1jXK1wXor5kvLeDpuStfKwx+3Rpce/zuKxIID5h0J4uWLq96jhFDB7Ppm9epuwYY0FoFbfnM3tn4LZd3Nr6kb0l6BXgB2BzYv0m1mVkXGz507R5rCwwfunaLK+oszsbPU07Z+GeRBuCNBE4Dfi7pQ47NNeteA1nBb+LNj3LuLSt3MM/OXdhjbYFn5y7s9RwHv2szxu2zbd01dCNn4+cpm2z8iJgDzAEelfQIMAPYHfhdk2o0sw42bp9tq3bUUx54usfaAmce9HY/s28gZ+PnKYts/CpKjxfWakZxZta9vLZAczkbP085ZOPvCrwTuBOYC7wFOAPYBHhb+Wh8D9Azsy7SlgP0rOWyzcZfCBwMnA6sCzxLWtP+ME+7MzMzG7hmDtAjIm4AboAVUbm/ImXj/zQini/aPAjs2cw6zMzMupmz8c3M8uTb+FaPPLPxe1QoDZb0Z0khaaeyXc7GNzMzq1Mu2fgl3yHl3/fQbdn4zu02s3bmbPz8tDwbv+yYj5Ke3R8MfLgZdbUD53abWTtzNn6essjGLz4c/AH4EPAP4Elg54i4r+K48TR56t1Asrj70ltWd7n+5HaXtDi/28yaoy2f2Tsbv+XyzMaXNIh0W//siPhzk+ppG87tNrN25mz8POWQjf/V4s/nNLGWmjXzCre3rO5y/cntLnF+t5nlwtn4ecohG39v4L3Aq1KPuw/3SrosIsY2qcbVrres7nLO7TazduZs/DzlkI1/NCk5r2QT4EZgLK/Nre8azu02s3bmbPw8NaWzj4iZxSC9SZKWsHI2/kURcX3R9snyYyW9XPzxbxGx0jS8bjBsyGB38mbWtsaOPcKde2ZyyMY3MzOzJmrWaHwgZeNHxF4RMYTU2d9BysbfcBXHPBURqpx2Z2ZmZvVxNr6ZWZ7acp69tVy+2fiSniry8Mu/zixr4mx8MzOzOuWUjX86MLzsa0JpR7dl41vzef0Bs+ZxNn5+ssnGB+ZHxKxm1GNWzusPmDWPs/HzlEs2/lPAYOD1wAzSPP2zImJJxXHjaXI2vrW/Va1vMP2FV5i74NWq6w/0tsaA1w2wFmnLZ/bOxm+5PLPxC+cBh5NWvTsfGEe65W/WUEuXh9cfMGsiZ+PnKYdsfCKiPBf/IUnzgMsk/WdEvNCU6qxjrepKfMJ1U/n5vdOrrj/gNQbMBs7Z+HnKIRu/mt8X37cG3Nlbw5yy/yi2H/Emrz9g1iTOxs9TDtn41exQfH+2CbVZl/P6A2bN42z8PDVlgF6x7VhgEnA2K2fjXx4RxxTtdgV2AW4D/gHsTBrJf19EfLTifcbjAXpm1h3acoCetVzV35scsvEXA4cBpwFrAdOBi0gfCszMzGyAmjlAj4i4AbgBVkTl/oqUjf/TiHi+aPMn0pW9mZmZNYGz8c3M8uTb+FaPfLPxy9reI2mBpLmSbi3b7Wx8swZzZLA1i+Ny89Oo2/grsvFLGyoG6H2e1MnvC3yLdNv+8LK2Y4CLgZOBo0gfQt5Z2l8k6d0r6YM4G99swBwZbM3iuNw8NSsud1Pgb8D3I+ILFW3HAFcDh0bEFZIGAU8CZ0TERX28z3g8Gt+sJr3FBtcTGVzi6ODVqi1v4zsut+VWa1zuIcCaVBlRHxFTgMeA0ke8dwGbAUsk/UnSLEk3SdqxSbWZdTVHBlszOS43T80ajd+fbPytiu+nA18iXeV/Frhd0nYR4WAdszr0dhXuyGBrJsfl5imHbPzS3YWvR8SVAJKOAT4A/DvpGb+ZNYgjg62ZHJebpxyy8UtX7lNLOyNiqaTHgM2bVJ9ZV3NksDWL43Lz1Kxn9leQrtxPrNxRlo0/udh0PylF761lbdYA3kJK0zOzJhg2ZDDvGvlmd/TWcGPHHsGTj01j+bJlPPnYNHf0GWjKlX1EzJQ0DpgkaQkrZ+NfFBHXF23nSfoB8DVJM0nz6I8D3gz8rBn1mZmZdZMcsvEp2pQ+FKwD/AnY04PzzMzMBq5Zt/GBlI0fEXtFxBBSZ38HKRt/w4p2r0bESRGxcUS8MSL2KDLzzczMbIDaJRt/c157fr88IgbVcJhDdcysnbVlqI61XJ7Z+JL2kBS9fB1SNHumOPcEwEHeHcYZ7Wadxdn4+ckhG/9uYHjF+T4HHE+xPG5ELCVl42+HP012FGe0m3UWZ+PnqeXZ+L2c71Hg9og4pmL7UcD5EfGGasdV8IeCBuotZ71WM19awNNzV75yH/y6NapmtJf0ldVeydnt1kHa8ja+s/FbLtts/B4k7QFsA/ywSbVZRpzRbtZZnI2fpxyy8SsdAzwYEfc1pTKry0CvmCfe/Cjn3rLy/+zPzl1YNaO9xFntZu3F2fh5yiEbfwVJ6wP/CnyxKRVZy4zbZ9uqnfaUB552RrtZB3E2fp5yyMYv9+/AMuCSJtVlmXFGu1lncTZ+nloxQO9A4Cpgv1Jkbtm+h4E/RsRRvbzPUXiAnpl1h7YcoGctV/X3puXZ+Cuqk3YnXfEfU3k+MzMzq18u2fgAnwQeiYi7mlWTmZlZN2pGXO6u5eE6ZW3WBX4FbEBa5Ob5fpx/C+DJ4qXjcs2sG/g2vtWj6fPsl5Nu1b/2jmVxucAcYH3SErbvq2i3raQpkuZImi/pXkkfLGsyozj3sgbWa2Zm1hUamY1/KfAOilH2RVzu9cD9wG6k5/HnA7sDB1ccfx0wGNgb2BG4E7hG0lsAImJZce4LgMcbVLNZW/DaAdZunI2fn0Zm4/8jIh6HFaPxJwKTIuKksnYXSpoNXC3pqoi4QtIGpMS8T0XEn4vjvwyMI3X8fwOIiMclvQi82qCazbLntQOs3TgbP0/Nmno3DjgHGFEtRa/Ivn84Ig6UJOBh4F7SAjgLgU8B3wDeFhHPlh03Hjg4IravoSw/T7LsrWrNgekvvMLcBa9WXTtgVWsGeH2AjtGWz+ydjd9yqzUbv+a43EifNvYBtgfmAYuB8cCHyjt6s26zdHl47QBrO87Gz1PL43KLK/sLgBeA95Ku7P8D+IWknXtJ4DPrCKu6Cp9w3VR+fu/0qmsHeM0Ay5Wz8fOUQ1zuXsABwHoRMbfYdqykfYCjgQlNqtEsa6fsP4rtR7zJawdYW3E2fp6a1dlfAZxJCtT5QvmOIi53a+DzxabSg8fl9LSc5j1mMGsLXjvA2o2z8fOUQ1zuPcCLwMWSTifdxv8ksBVpSp5ZVxs2ZLA7eWsrY8ce4c49My2Py42IOUWAzteBW4HXkwbwjYmIPzWrPjMzs27RzAF6RMQNwA3QIy73aEk/LY/LjYj7gH9pZi1mZmbdqhnZ+J+IiIt7abcu6Vn9tIi4qh/n3w0oLZCzOCJquafpefZm1s7acp69tVxT59lPBa4EdgUuX/GOZdn4khYCdwOvAFN6VCa9U9LNkuZKekHSDyWVr1l/X3HuH5Oy9c3MzKxGjerslwCzI+LeiHgFes3GvwD4GnBJ6UBJmwC/AZ4A/hn4IDAamFxqExFLipX0ZgJLG1SzdSlnzZs1l7Px89OUZ/b9ycYH9idNszu2WPAGSZ8GHpK0dSlv36wRnDVv1lzOxs9TDtn4xwNfjYjhZfu3Bh4Djo6IyWXbx+Ns/I60qoz4/pr50gKenlv9qn3w69aomjVfsqrM+WqcQ29N1JbP7J2N33J5ZuOTptttIOnLktaU9GZSIA/A8KpHm9XJWfNmzeVs/Dy1PBs/Ih6WdCTpTsDXSc/kzwNms3KqnnWoRl4hT7z5Uc69pfo/LM/OXVg1a77EmfNmA+Ns/DzlkI1PRFwKXCppI9Jo/QC+SBq0Z9Yv4/bZttcOe8oDTztr3qyJnI2fpxyy8VeIiNlFm08Ai4Cbm1SfdSlnzZs1l7Px85RDNj6SjiNl5M8v2pwFfLlsFTyzhnHWvFlzORs/Py3Pxi+8mzT//g3ANOBTEfGzZtVmZmbWTXLJxv/3ZtZhZmbWzWqaeifpdklRfO3SS7NjJL0saedqO4tkvf2AS4H39qdISZtLWgScBoyWdH5/jjczM+tm/ZlnfzFp3vv9pQ2l7HvgXaRpcjOA90rqcV5JJ0u6C3gO+Hq1RXCKDv1aSa9ImiPpPElrFrufKd5j1/L3NzMzs771p7NfEBGzIuJVWCn7flfgbcB3gVMpy74vrAVcVexfiaRBpFv8Q0hX/YcDBwNnA0TE0oh4uMjHf7kfNVuTOWfezCo5Gz8/dT2z72f2PRFxanHcwb2ccl/S4jcjI2JG0fYk4EeSTo6IefXUac3lnHkzq+Rs/DzVlI0/kOz7iu0HA1dEhCq2nw4cFBGjy7ZtSLrtv1dE3NZbLavgbPxCI3Lnq+XN95UzD/3PmgfnzZsVnI1v9WhoNn5/su9rsTEpHrfcHGBZsc8y5Jx5M6vkbPw8DWTqXU3Z99Z6jbhSrpY331fOPDhr3qzbOBs/T/V29v3Kvq/BLOA9Fds2AAYV+6zFquXNO2fezCo5Gz9P9Xb2dWXfr8I9wCmSNo2ImcW2fYDFeKpdtpwzb2aVnI2fp7oG6BXbjgUmkabHVWbfXx4Rx5S13RxYD/gAKfd+x2LX4xHxcjH17kHgeeBLwPrFOa+KiOP7qqUXHqBnZu2sLQfoWctV/b2p+5l9P7PvTweOLHtdusW/J3B7RCyTtB9wAXAXsJA0V//EeuszMzOzpO4r+yptStn3GwB7lmffN5Kv7M2sS/jK3uoxoKl3OwCfLbLxj67WYCDZ97WQNEFSAO8v3sfMzMxqUGtnPxW4khSLe3lpYykbX9I8SQuBu4FXgCnlB5ey8Yvc+6qfBiWdK+k+SYskPVWlyTnF+19Pmg1gZmZmNai1s18CzI6Ie4sr+Mps/N1I0+0uIK1L369s/LJafgL8tNrOiHixyMZ/Hni1xrrNGsprAZj1zdn4+cklG5/SqHtJJ5Cy8s2y4rUAzPrmbPw8ZZGNX9HmBOC4iNiil/2TgQ0iYv8+yvbgEetTLesGlNYF6G0tgFXl/zvn3wagLQfoORu/5bLOxjfLntcCMOubs/Hz5Gx862q1XHmX1gXobS0A5/+bvcbZ+HnKJRvfLFuldQG8FoBZ35yNn6dcsvHNsue1AMz65mz8PNXV2UfEzGKQ3iRJS1g5G/+iiLi+1L4sG3+L4vUOxa7HI+LlYtvWpMjdTYA1y9pMjQg/ErAsDBsy2J28WR/Gjj3CnXtmssjGL/78I1I6XmWbLYGn6q3TzMys27VjNv5kPPXOzDpfW069s5ar+ntTa2c/F3hT8fITEXFxL+3WJT2rnxYRV9VXZ681HA+cV7x8LiI26uMQ/yKaWTtzZ2/1GNA8+6Zm40t6h6T/kTRD0kJJf5V0kqTy+n7Ca9n4D9VYt5mZWdfLJRv/XaTM+48Do4HTgP8CvlxqEBHzyrLxF9dYt2XOWfNmncfZ+PnJIhs/Iv67YtMTkt4JHAR8o54aLX/OmjfrPM7Gz1N22fhlbX8AbB0RH6jYPhkP0MtOLRnz8FrOPNBr1jywyrz5EufOW4dry2f2zsZvufbJxi+u6o8Cvl/vOSx/zpo36zzOxs9Tdtn4kt5KmsL33Yj4RT3nsNWv1qvsUs480GvWPDhv3qxdORs/T1ll40vaDrgN+N+I+HJf7a39lHLmAWfNm3UgZ+PnKZtsfEmjgFuByyNiXJ11WRtx1rxZ53E2fp7qTtArpt5NAs5m5Wz8yyPimLK2pWz8DwBnATsWux6PiJcljSZ19LdR8eEhImZV1DIZD9Azs87XlgP0rOWq/t7kko1/CDAMOKz46rNwMzMzq42z8c3M8pT7hY7/jc2Ts/HNzNqIO3urR9bZ+BtKulHSM5IWFxn535P0prJmzsY3ayJHF1ujOC43P7lk4y8HrgYOIAX2HAXsDVxUauBsfLPmmfLA0+w78becfu3D7Dvxt0x5oNqMWrO+leJyF+50JJt98SoW7nQknxl3ojv8Fqs3LndT4G/A9yPiCxVtx5A67kNL2fhl+/oTl/s54CsRMbxi+2T8zN6sT7VEGJfii1cVXQy1xReXOMa4YdryNr7jcluuoXG5hwBrkqbZ9RARU4DHgLonVUraBPhX4I56z2FmtXN0sTWK43LzVO/Uu6Zk40v6H+CjwNrAdcDRddZn1vVqucIuxRevKroYHF9stXNcbp5yy8YfR3rmvy3wTdIz/k/VcR4zq0EpvtjRxdYojsvNU1bZ+EVa3ixgmqQXgd9JmhARM+qs08xq4OhiaxTH5eYpm2z8KkrjCdYa4HnMrAbDhgx2J28NMXbsEe7cM1NXZx8RMyWNAyZJWsLK2fgXRcT1pfZl2fhbFK93KHaVsvH3B9YnTeN7GRhNytC/NyIer6dGMzMzS3LJxl8EfBp4G+lKfgZp+t6Z9dZnZmZmST1xubsW4TaVbZqajS9pD9KqeOC4XDPrfG05z95absDz7JeTbtW/dsayuFxgDulW/FPA+yrarTIut6LtBpKelhSSNqj4ARbgXzAzM7N+6U82/qXAOyhuwfcSl3s+sDtwcMXxfcXllrsYeLDK9nuL978aZ+PbADgD3qy5nI2fn1qf2S8B/lEaLFfE5U4EJkXESWXtLpQ0G7ha0lWluNyIOLU4rvJDQA+SPg+sA3wd+HD5vohYCDwuaT4eoW91mvLA04y/9uEV88nHHzDa88nNGqiUjb/O3sex2ZhRLJw5Nc27B4/Qb6F6s/HHAecAI6ql6El6FHg4Ig6s2N5rNr6kHUnP/HcGtiE9n98wIuZUtJuMs/FtFXrLhJ/+wivMXfBqjwz4j5x/F9ttPIQ1X/faTS5nu1sm2vKZvbPxW66h2fgNjcstBvf9L3B8LyE9ZgO2dHmslAG/ydDBLF66rMWVmXUOZ+PnKZe43POAOyPiFwOoxwzo/cp8wnVT+fm903tkwM+et5hbT3i/w2TMGsTZ+HnKJS53b2AzSaW5+KXbELMkfSsiTq6zTrMVTtl/FNuPeFOPDPgJY7Z3R2/WQM7Gz1Mucbn7kpbMLdkZ+G9gD9JyuWYN4Qx4s+ZyNn6esojLjYhHy89fNr9+WuUAPbOBcga8WXM5Gz8/ucTlmpmZWZPUNfWulzZNjcste5/JeOqdmXW+tpx6Zy1X9femnmz8T0TExb20W5f0rH5aRFxVX5291nA8adQ+OBvfzDqfO3urx4Dm2U8FrgR2BS5fccaybHxJC4G7gVeAKT3euYZs/CILv/Lr02VNflK8//U4LtfMzKxmtXb2S4DZEXFvRLwCvWbjXwB8Dbik4vhas/E/CQwv+/pJaUdEzCtW23seWFxj3dbGnGFv1p6cjZ+fugboNSsbH5gbEbPqqck6izPszdqTs/HzlFM2fgDPkO4CPAn8GPhhRCyvaDcZD9BbrXrLmh+omS8t4Om51a/aB79ujR4Z9gdMupNFS5ev1G7E0MFs+uZ1+vW+zr63NtGWz+ydjd9y+WbjF04FDgM+QMrJPxv4ap31WZurzLAfPnTtFldkZrVwNn6ecsnGJyLOKHv5oKRBwMnAhP4WZo3VrCvhiTc/yrm3VP8H4Nm5C3tk2D87d2HVdge/azPG7bNtU+ozs/5zNn6ecsnGr+b3wBslbRQRswd4LsvQuH227bWjnvLA0z0y7M886O1+Zm/WBpyNn6dcsvGr2QFYBMwd4HmsDTnD3qw9ORs/T1lk40s6ANgYuAdYSIrRPZ00QM/T7LqUM+zN2pOz8fOTSzb+q8CxpBH+awBPkAbsfa/e+szMzCxxNr6ZWZ7acuqdtZyz8c3M2og7e6tH3tn4RbuPSXpQ0iJJcyT9tGy3s/HNzMzqkE02vqTPAWcB3wG2Jz3Pv6a039n47cW59mbdy9n4+ckiG1/SUOCbwJiIuLls1//VU5+1lnPtzbqXs/HzlEU2vqRDgZ8C/0GKyH0T8AfgSxHxREXbyXiA3oA0Iuu+Vbn2Jc63ty7Qls/snY3fclln429V1HIK8EXgQOD1wG2S6usNrGWca2/WvZyNn6dcsvHXIHXun4uImwAkjQVmAQcAl9VVoVXViKti59qbWTXOxs9TLtn4zxbfp5Y2RMQ/JD0DbF5njdZEzrU3s2qcjZ+nXLLx7yq+vxWYWZznDcBwYHqdNVqLONferHs5Gz9PdSfoFVPvJpHWna/Mxr88Io4pa1vKxv8AaXrdjsWuxyPi5aLNFNKHhE8BL5Gm8O0MjIqIBWXnmowH6JlZ52vLAXrWclV/b3LJxgf4OGmE/7VFsXcCe5d39GZmZtZ/zsY3M8uTr+ytHvVNvSs6+vcDn5W0QNLO1doVyXr7AZcC762/zqo1bCEpiqjdI/s8wMzMzFaodZ79ZaRc+ndS5NJX5uJL+jMpFOfM8kVw+srFl3RUqSOv8lX6YDED2BL4MfAY8Mm6f2IzM7MuU2tnP6fIxZ8WEYsbnIt/GWnUffnXz0lr2t8HEBHLIuIp0kj9JRHxbPVTmVlfvG6BNZuz8fPT7wF6jc7Fj4iFwIrUlSIx7wDg21HLgAIzq5nXLbBmczZ+nvocoNfsXPwqxx8FXARsFhGzKvaNBw6OiO1X/WMBHjxiXaCWdQ7K1zFY1boF9a5V4HUKmqYtB+g5G7/lGpaN3+hc/ErHANdVdvRmNnBet8Cazdn4eap3nn0jc/FXkDSaNBBwv3qON+s2tVxVl69jsKp1C7xWgTWCs/HzVE9n3+hc/HLHkEbe/7rO482sQvk6Bl63wJrN2fh5qqezb3QufunYwaQUvfMiYuXFz81swLxugTWbs/Hz1O/OPiJmFoP0Jklawsq5+BdFxPWl9mW5+FsUr3codq3IxS8cDLwJ+O/+/xhmVqthQwa7k7emGjv2CHfumanrmX0TcvEhBeXcGBF/r6cmMzMzq66e0fgARMQNEbFXRAwhdfZ3AEdL2rCi3VERoSpft1e0e39EfLjeeszMzKy6Wjv7Y8sibHep3NnMXHxIjwIkLQROA0ZLuq7R72FmZtapaunsxwJ/BK4ERpIicoGe+fjAHOAwYDNJPc7bVz5+0WZnSb+RNLf4ukXSu4vdz5By+XclLYFrZmZmNeqzsy+m1y0EZkfE3yPiVYBG5uNLegNput0zwC6kTv1Z4EZJQyJiaUQ8EhH3Ai/294e07uDMd7M8OBs/P3UN0Gt0Pj6wHWnE/mkR8WTR9r9IdxXeSrEgjllvnPlulgdn4+epz2x8aH4+vqQhwN+AHwJnFJtPI3X22xWL5ZTaTgY2iIj9+yjb2fhtqpa8d6gt872vvHfnulvGnI1v9WhYNj40OB8/IuYDewCHkubsLyA9/9+nvKM3640z383y4Gz8PNWbjQ8NzMeXtDYpTOde0tX8IOAE4BpJOxWj/a1L1Hq1XUvmu/PezVYvZ+Pnqd7OvtH5+EcAbwHeExHLACQdAbwEHAj8vM46rYM5890sP87Gz1O9nX2j8/HXId0pKM/EX15sqzv4x7qHM9/N8uBs/DzVG5fb6Hz8m4GzgAsknUfq4L8MLANuradG6z7OfDfLg7Px81P3M/tG5uNHxDRJB5BG4N9DuqJ/EPhQRMyst0YzMzOrc+pdL23WBX4FbADsGRHPN6rIiveZjKfemVnna8upd9ZyVX9vau3s55KWnwX4RERc3Eu7dUnP6qdFxFX11dlrDccD5xUvn4uIjfo4xL+IZtbO3NlbPQY0z34qKRt/V+DyFWcsy8YvFqq5G3gFmNLjnWvLxt9b0t2S5kuaJelbksofM/ykeP/rgYdqrLtrOTrWzFrFcbn5qbWzX0LKxr+3NOe9wdn47yjOdTOwIylQ5yOkEf8ARMS8Ihv/eWBxjXV3pSkPPM2+E3/L6dc+zL4Tf8uUB6rNjjQza7xSXO7CnY5ksy9excKdjuQz4050h99i9cblbkqKt/1+RHyhou0Y4Grg0FI2ftm+3uJyv0EajLdj2bYDSHcRhhUJe6Xtk2mjZ/a1Rr/2R3lMbDW9RcdW01ecbC0cOWvWFG15G99xuS3X0LjcQ4A1SdPseoiIKcBjpKCcWq0FVPZeC4HBwLvqK7F7OTrWzFrFcbl5qnfqXUOz8YEbgXGSPgb8L7ARcGqxb3idNWahGVe95TGx1fQWHVuN42TNrJEcl5unLLLxI+ImSScA3wMmk57JnwG8l56pekbPmNhqHB1rZq3iuNw85ZKNT0ScI2ki6Ur+JVLa3jeBJ+qssWs5OtbMWsVxuXlqxgC9A0kj7/crj8wt9lUdoNfLe54OHAVsWVocp9g+mTYaoGdmVqe2HKBnLVf19yaXbHwknQj8mnTb/l9J2fiHlnf0ZmZm1n9ZZOMXf/4QcDJpZP6fgY9GxA311mdmZmaJs/HNzPLk2/hWjwHNs98B+KykkHR0tQZFst5+wKWkUfQNJekbRdTukcDoRp/fzMysUzU9G1/SFpJ+LOkJSQuL79+U1CPpRdLmkq4t8vPnSDpP0pplTc4q3v8B4LZ6fljrP2fsm1l/ORs/P7U+s1+RjV/aUGTjTwLOJq109wqwL/AtYBfg8KLpdsAg0nP9x4C3AT8E1geOKc41iPQI4AXSXYH1SYP+BBwPEBEvAfdKmkcaDGhNNuWBpxl/7cMr5uuPP2C05+ub2SqVsvHX2fs4NhszioUzp6Z59+Dpdy20WrPxy9ocC5wREesXrz9E6uxHRsSMYtvHgB+RsvHn9VbLKnTk86RGZu3nlLHvfH2zlbTlM3tn47dcVtn4byQF55TsCjxS6ugLN5JG5jsbv0WcsW9m/eVs/Dyt9mx8SSOBE4BvlG3eGJhd0XQOsKzYZ4VGXgE7Y9/MGs3Z+Hlardn4kjYiBefcDEwcwHtbAzhj38wazdn4eVpt2fiSNgZuBf4CfDx6DhaYBbyn4hwbkAb2zaqzRhsgZ+ybWX85Gz9PqyUbX9Jw0nS5h4HDImJpxTGlAXqbR8TMYtsRwH/jAXpm1p3acoCetVxrsvElbUKKxH0G+AKwgbSilueL7PubSB8EfirpS6Spd2cV51nR0ZuZmVn/rY5s/H2BbYqvv1ecZkvgqYhYJmk/4ALgLmAhcElxbjMzMxuAdszG77OWot2vi1r6YwPSLADrP//dDYz//gamE//+5kTEB1tdRG9q/De2E/+7tFItf59Vf2/609nvRhphv2dE/LGXduuS0vSmRcRVfZ64HySNBS4E1iaNFejrmX0973FfROzU6PN2A//dDYz//gbGf3958n+XxhrI32ett/HHkjpZgBm9NSoWw/lGb/sH6JfA74s/z23Se5iZmXWcmjr7XqbXrVYRMR+Y3+o6zMzM2k29cbmd6oetLqCN+e9uYPz3NzD++8uT/7s0Vt1/nzU9szczM7P25St7MzOzDufO3szMrMO5s18FJTdICkkHt7qediBpPUmTJE2TtFDSDEnfl7R+q2vLlaRjJT0paZGk+yW9t9U15U7SVyT9UdI8Sc9LulbS9q2uy0DS+yT9UtLTxb+dR7W6pnbWqN91d/ar9iVgeauLaDObACOAk4B/Aj4GvA/4n1YWlStJhwHnkqas7gjcDdwgafOWFpa/PUiJm7sBewFLgd9IWq+VRRmQ0lT/Qspc6X1dbKvVHjTgd90D9HohaWfSgj7vAmYDh0TEla2tqj1J+jBwHTDUax30JOn3wEMR8cmybY8BV0bEV1pXWXuR9AbgH8CYiLi21fVYIull4LiImNzqWjpFvb/rvrKvQtIQ4FLgmIh4rtX1dIA3AotJiyVZQdKapA+TN1Xsuon0Kd5qN4T079lLrS7ErMnq+l13Z1/dD4BfR8QNrS6k3UkaCpxBWsFwaR/Nu80GwCDSnaNys4GNV385be1c4EHgnhbXYdZsdf2ud01nL2lCMVhkVV97SPo48A684l4Ptf79VRzzBuBa4GnSM3yzhpN0DrA7cFCxZLZZRxrI73rdS9y2oe8CP++jzd+Bo4BRwMuSyvddJumeiNi9KdXl77vU9vcHrOjory9e7h8Ri5pUVzubAywDNqrYvhEwa/WX034kTQT+jbRA1xOtrsesWQb6u941nX1EzKGGpRYlnQx8p2Lz/wEnANc0obS2UOvfH6wY83ADIOCDEfFyM2trVxGxRNL9wD7AFWW79gF+0Zqq2oekc4HDSP/4TWt1PWbN0ojf9a7p7GtVLPrTY+Gf4gp/hq8c+lZ09DeRBuWNAdYtlj4GeDEilrSqtkydA/xM0h+Au4BPk6Yv/qClVWVO0veAj5N+x16SVBrj8LI/XLZWcVdv6+LlGsDmknYg/f//914PtKoa9bvuqXc1kBR46l1Niuf2t/Wye8+IuH21FdMmJB1LGtMwnDQ/eVxE/La1VeWt+H+ymq9FxPjVWYv1tIp/A34SEUet1mI6QKN+193Zm5mZdbiuGY1vZmbWrdzZm5mZdTh39mZmZh3Onb2ZmVmHc2dvZmbW4dzZm5mZdTh39mZmZh3Onb2ZmVmHc2dvZmbW4f4/NatR2MAjJPQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 576x648 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"az.plot_forest(idata, var_names=['a', 'b'],\n",
" combined=True, hdi_prob=0.95, r_hat=True)"
]
},
{
"cell_type": "markdown",
"id": "05ec12bb-5ca7-4a26-b39a-49df7092a37a",
"metadata": {},
"source": [
"推定結果の統計情報を表示してみる．"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "447c3bc4-188a-4076-9100-b7aaa143dec1",
"metadata": {},
"outputs": [],
"source": [
"summary = az.summary(idata, round_to=3)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "fd74a40a-2347-4ab7-9a52-263a40d4b1fc",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>sd</th>\n",
" <th>hdi_3%</th>\n",
" <th>hdi_97%</th>\n",
" <th>mcse_mean</th>\n",
" <th>mcse_sd</th>\n",
" <th>ess_bulk</th>\n",
" <th>ess_tail</th>\n",
" <th>r_hat</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>a[Q1]</th>\n",
" <td>0.443</td>\n",
" <td>0.065</td>\n",
" <td>0.320</td>\n",
" <td>0.560</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4748.626</td>\n",
" <td>6279.968</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q2]</th>\n",
" <td>0.705</td>\n",
" <td>0.085</td>\n",
" <td>0.544</td>\n",
" <td>0.866</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4115.761</td>\n",
" <td>6039.690</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q3]</th>\n",
" <td>0.499</td>\n",
" <td>0.071</td>\n",
" <td>0.364</td>\n",
" <td>0.628</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4478.687</td>\n",
" <td>5916.186</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q4]</th>\n",
" <td>0.460</td>\n",
" <td>0.063</td>\n",
" <td>0.350</td>\n",
" <td>0.583</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4952.672</td>\n",
" <td>6566.816</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q5]</th>\n",
" <td>0.613</td>\n",
" <td>0.081</td>\n",
" <td>0.466</td>\n",
" <td>0.768</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>3974.533</td>\n",
" <td>6185.484</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q6]</th>\n",
" <td>0.542</td>\n",
" <td>0.076</td>\n",
" <td>0.398</td>\n",
" <td>0.684</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4140.866</td>\n",
" <td>5837.839</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q7]</th>\n",
" <td>0.413</td>\n",
" <td>0.066</td>\n",
" <td>0.292</td>\n",
" <td>0.538</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4685.018</td>\n",
" <td>5705.328</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q8]</th>\n",
" <td>0.374</td>\n",
" <td>0.072</td>\n",
" <td>0.237</td>\n",
" <td>0.509</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4649.939</td>\n",
" <td>5844.028</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q9]</th>\n",
" <td>0.249</td>\n",
" <td>0.062</td>\n",
" <td>0.138</td>\n",
" <td>0.365</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4458.675</td>\n",
" <td>4787.609</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q10]</th>\n",
" <td>0.505</td>\n",
" <td>0.072</td>\n",
" <td>0.368</td>\n",
" <td>0.639</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4382.538</td>\n",
" <td>6283.526</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q11]</th>\n",
" <td>0.390</td>\n",
" <td>0.070</td>\n",
" <td>0.260</td>\n",
" <td>0.524</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4375.445</td>\n",
" <td>6215.533</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q12]</th>\n",
" <td>0.749</td>\n",
" <td>0.092</td>\n",
" <td>0.586</td>\n",
" <td>0.929</td>\n",
" <td>0.002</td>\n",
" <td>0.001</td>\n",
" <td>3230.900</td>\n",
" <td>5453.242</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q13]</th>\n",
" <td>0.391</td>\n",
" <td>0.071</td>\n",
" <td>0.258</td>\n",
" <td>0.524</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4076.469</td>\n",
" <td>4133.558</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q14]</th>\n",
" <td>0.814</td>\n",
" <td>0.095</td>\n",
" <td>0.640</td>\n",
" <td>0.994</td>\n",
" <td>0.002</td>\n",
" <td>0.001</td>\n",
" <td>3202.875</td>\n",
" <td>5410.097</td>\n",
" <td>1.002</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q15]</th>\n",
" <td>0.289</td>\n",
" <td>0.060</td>\n",
" <td>0.176</td>\n",
" <td>0.399</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4309.694</td>\n",
" <td>5209.750</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q16]</th>\n",
" <td>0.375</td>\n",
" <td>0.062</td>\n",
" <td>0.264</td>\n",
" <td>0.496</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4622.702</td>\n",
" <td>5722.293</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q17]</th>\n",
" <td>0.329</td>\n",
" <td>0.069</td>\n",
" <td>0.203</td>\n",
" <td>0.459</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4279.217</td>\n",
" <td>4921.235</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q18]</th>\n",
" <td>0.544</td>\n",
" <td>0.077</td>\n",
" <td>0.401</td>\n",
" <td>0.688</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>3672.073</td>\n",
" <td>4754.482</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q19]</th>\n",
" <td>0.473</td>\n",
" <td>0.066</td>\n",
" <td>0.350</td>\n",
" <td>0.596</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>3944.499</td>\n",
" <td>5503.309</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>a[Q20]</th>\n",
" <td>0.409</td>\n",
" <td>0.064</td>\n",
" <td>0.290</td>\n",
" <td>0.528</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>4301.130</td>\n",
" <td>5135.241</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \\\n",
"a[Q1] 0.443 0.065 0.320 0.560 0.001 0.001 4748.626 6279.968 \n",
"a[Q2] 0.705 0.085 0.544 0.866 0.001 0.001 4115.761 6039.690 \n",
"a[Q3] 0.499 0.071 0.364 0.628 0.001 0.001 4478.687 5916.186 \n",
"a[Q4] 0.460 0.063 0.350 0.583 0.001 0.001 4952.672 6566.816 \n",
"a[Q5] 0.613 0.081 0.466 0.768 0.001 0.001 3974.533 6185.484 \n",
"a[Q6] 0.542 0.076 0.398 0.684 0.001 0.001 4140.866 5837.839 \n",
"a[Q7] 0.413 0.066 0.292 0.538 0.001 0.001 4685.018 5705.328 \n",
"a[Q8] 0.374 0.072 0.237 0.509 0.001 0.001 4649.939 5844.028 \n",
"a[Q9] 0.249 0.062 0.138 0.365 0.001 0.001 4458.675 4787.609 \n",
"a[Q10] 0.505 0.072 0.368 0.639 0.001 0.001 4382.538 6283.526 \n",
"a[Q11] 0.390 0.070 0.260 0.524 0.001 0.001 4375.445 6215.533 \n",
"a[Q12] 0.749 0.092 0.586 0.929 0.002 0.001 3230.900 5453.242 \n",
"a[Q13] 0.391 0.071 0.258 0.524 0.001 0.001 4076.469 4133.558 \n",
"a[Q14] 0.814 0.095 0.640 0.994 0.002 0.001 3202.875 5410.097 \n",
"a[Q15] 0.289 0.060 0.176 0.399 0.001 0.001 4309.694 5209.750 \n",
"a[Q16] 0.375 0.062 0.264 0.496 0.001 0.001 4622.702 5722.293 \n",
"a[Q17] 0.329 0.069 0.203 0.459 0.001 0.001 4279.217 4921.235 \n",
"a[Q18] 0.544 0.077 0.401 0.688 0.001 0.001 3672.073 4754.482 \n",
"a[Q19] 0.473 0.066 0.350 0.596 0.001 0.001 3944.499 5503.309 \n",
"a[Q20] 0.409 0.064 0.290 0.528 0.001 0.001 4301.130 5135.241 \n",
"\n",
" r_hat \n",
"a[Q1] 1.000 \n",
"a[Q2] 1.000 \n",
"a[Q3] 1.001 \n",
"a[Q4] 1.000 \n",
"a[Q5] 1.000 \n",
"a[Q6] 1.000 \n",
"a[Q7] 1.000 \n",
"a[Q8] 1.000 \n",
"a[Q9] 1.001 \n",
"a[Q10] 1.000 \n",
"a[Q11] 1.000 \n",
"a[Q12] 1.001 \n",
"a[Q13] 1.001 \n",
"a[Q14] 1.002 \n",
"a[Q15] 1.001 \n",
"a[Q16] 1.000 \n",
"a[Q17] 1.001 \n",
"a[Q18] 1.000 \n",
"a[Q19] 1.001 \n",
"a[Q20] 1.001 "
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# パラメータa\n",
"summary['a[Q1]':'a[Q20]']"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "c5e19802-41de-4d4c-bbf7-b9cc9985613b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>sd</th>\n",
" <th>hdi_3%</th>\n",
" <th>hdi_97%</th>\n",
" <th>mcse_mean</th>\n",
" <th>mcse_sd</th>\n",
" <th>ess_bulk</th>\n",
" <th>ess_tail</th>\n",
" <th>r_hat</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b[Q1]</th>\n",
" <td>-2.543</td>\n",
" <td>0.363</td>\n",
" <td>-3.224</td>\n",
" <td>-1.895</td>\n",
" <td>0.005</td>\n",
" <td>0.004</td>\n",
" <td>4923.049</td>\n",
" <td>6378.865</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q2]</th>\n",
" <td>0.324</td>\n",
" <td>0.079</td>\n",
" <td>0.174</td>\n",
" <td>0.468</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>7776.026</td>\n",
" <td>6934.538</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q3]</th>\n",
" <td>2.357</td>\n",
" <td>0.321</td>\n",
" <td>1.793</td>\n",
" <td>2.968</td>\n",
" <td>0.005</td>\n",
" <td>0.003</td>\n",
" <td>4478.772</td>\n",
" <td>5965.396</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q4]</th>\n",
" <td>-3.139</td>\n",
" <td>0.409</td>\n",
" <td>-3.914</td>\n",
" <td>-2.418</td>\n",
" <td>0.006</td>\n",
" <td>0.004</td>\n",
" <td>4925.314</td>\n",
" <td>6661.555</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q5]</th>\n",
" <td>-0.884</td>\n",
" <td>0.131</td>\n",
" <td>-1.125</td>\n",
" <td>-0.638</td>\n",
" <td>0.002</td>\n",
" <td>0.001</td>\n",
" <td>5210.915</td>\n",
" <td>6823.221</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q6]</th>\n",
" <td>0.765</td>\n",
" <td>0.134</td>\n",
" <td>0.525</td>\n",
" <td>1.016</td>\n",
" <td>0.002</td>\n",
" <td>0.001</td>\n",
" <td>5311.182</td>\n",
" <td>6160.830</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q7]</th>\n",
" <td>-1.650</td>\n",
" <td>0.279</td>\n",
" <td>-2.166</td>\n",
" <td>-1.164</td>\n",
" <td>0.004</td>\n",
" <td>0.003</td>\n",
" <td>4916.670</td>\n",
" <td>5786.015</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q8]</th>\n",
" <td>0.212</td>\n",
" <td>0.132</td>\n",
" <td>-0.034</td>\n",
" <td>0.460</td>\n",
" <td>0.001</td>\n",
" <td>0.001</td>\n",
" <td>10463.261</td>\n",
" <td>7028.130</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q9]</th>\n",
" <td>-0.964</td>\n",
" <td>0.311</td>\n",
" <td>-1.573</td>\n",
" <td>-0.470</td>\n",
" <td>0.005</td>\n",
" <td>0.004</td>\n",
" <td>5363.273</td>\n",
" <td>4611.401</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q10]</th>\n",
" <td>-2.049</td>\n",
" <td>0.288</td>\n",
" <td>-2.599</td>\n",
" <td>-1.545</td>\n",
" <td>0.004</td>\n",
" <td>0.003</td>\n",
" <td>4585.992</td>\n",
" <td>6487.782</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q11]</th>\n",
" <td>-0.550</td>\n",
" <td>0.158</td>\n",
" <td>-0.857</td>\n",
" <td>-0.280</td>\n",
" <td>0.002</td>\n",
" <td>0.002</td>\n",
" <td>6290.421</td>\n",
" <td>5415.573</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q12]</th>\n",
" <td>-1.642</td>\n",
" <td>0.185</td>\n",
" <td>-1.989</td>\n",
" <td>-1.305</td>\n",
" <td>0.003</td>\n",
" <td>0.002</td>\n",
" <td>3419.851</td>\n",
" <td>5278.103</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q13]</th>\n",
" <td>0.882</td>\n",
" <td>0.203</td>\n",
" <td>0.529</td>\n",
" <td>1.256</td>\n",
" <td>0.003</td>\n",
" <td>0.002</td>\n",
" <td>5377.008</td>\n",
" <td>4709.581</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q14]</th>\n",
" <td>-0.921</td>\n",
" <td>0.108</td>\n",
" <td>-1.133</td>\n",
" <td>-0.729</td>\n",
" <td>0.002</td>\n",
" <td>0.001</td>\n",
" <td>3846.852</td>\n",
" <td>5899.100</td>\n",
" <td>1.002</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q15]</th>\n",
" <td>-1.506</td>\n",
" <td>0.345</td>\n",
" <td>-2.180</td>\n",
" <td>-0.925</td>\n",
" <td>0.005</td>\n",
" <td>0.004</td>\n",
" <td>4605.548</td>\n",
" <td>4840.131</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q16]</th>\n",
" <td>-1.908</td>\n",
" <td>0.328</td>\n",
" <td>-2.514</td>\n",
" <td>-1.331</td>\n",
" <td>0.005</td>\n",
" <td>0.004</td>\n",
" <td>4833.553</td>\n",
" <td>5919.292</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q17]</th>\n",
" <td>0.719</td>\n",
" <td>0.213</td>\n",
" <td>0.353</td>\n",
" <td>1.123</td>\n",
" <td>0.003</td>\n",
" <td>0.002</td>\n",
" <td>5679.487</td>\n",
" <td>5144.036</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q18]</th>\n",
" <td>-1.182</td>\n",
" <td>0.177</td>\n",
" <td>-1.516</td>\n",
" <td>-0.868</td>\n",
" <td>0.003</td>\n",
" <td>0.002</td>\n",
" <td>4175.659</td>\n",
" <td>4587.220</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q19]</th>\n",
" <td>-2.717</td>\n",
" <td>0.366</td>\n",
" <td>-3.407</td>\n",
" <td>-2.059</td>\n",
" <td>0.006</td>\n",
" <td>0.004</td>\n",
" <td>4063.181</td>\n",
" <td>5383.458</td>\n",
" <td>1.001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>b[Q20]</th>\n",
" <td>-2.268</td>\n",
" <td>0.356</td>\n",
" <td>-2.965</td>\n",
" <td>-1.665</td>\n",
" <td>0.005</td>\n",
" <td>0.004</td>\n",
" <td>4519.021</td>\n",
" <td>4649.440</td>\n",
" <td>1.000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \\\n",
"b[Q1] -2.543 0.363 -3.224 -1.895 0.005 0.004 4923.049 \n",
"b[Q2] 0.324 0.079 0.174 0.468 0.001 0.001 7776.026 \n",
"b[Q3] 2.357 0.321 1.793 2.968 0.005 0.003 4478.772 \n",
"b[Q4] -3.139 0.409 -3.914 -2.418 0.006 0.004 4925.314 \n",
"b[Q5] -0.884 0.131 -1.125 -0.638 0.002 0.001 5210.915 \n",
"b[Q6] 0.765 0.134 0.525 1.016 0.002 0.001 5311.182 \n",
"b[Q7] -1.650 0.279 -2.166 -1.164 0.004 0.003 4916.670 \n",
"b[Q8] 0.212 0.132 -0.034 0.460 0.001 0.001 10463.261 \n",
"b[Q9] -0.964 0.311 -1.573 -0.470 0.005 0.004 5363.273 \n",
"b[Q10] -2.049 0.288 -2.599 -1.545 0.004 0.003 4585.992 \n",
"b[Q11] -0.550 0.158 -0.857 -0.280 0.002 0.002 6290.421 \n",
"b[Q12] -1.642 0.185 -1.989 -1.305 0.003 0.002 3419.851 \n",
"b[Q13] 0.882 0.203 0.529 1.256 0.003 0.002 5377.008 \n",
"b[Q14] -0.921 0.108 -1.133 -0.729 0.002 0.001 3846.852 \n",
"b[Q15] -1.506 0.345 -2.180 -0.925 0.005 0.004 4605.548 \n",
"b[Q16] -1.908 0.328 -2.514 -1.331 0.005 0.004 4833.553 \n",
"b[Q17] 0.719 0.213 0.353 1.123 0.003 0.002 5679.487 \n",
"b[Q18] -1.182 0.177 -1.516 -0.868 0.003 0.002 4175.659 \n",
"b[Q19] -2.717 0.366 -3.407 -2.059 0.006 0.004 4063.181 \n",
"b[Q20] -2.268 0.356 -2.965 -1.665 0.005 0.004 4519.021 \n",
"\n",
" ess_tail r_hat \n",
"b[Q1] 6378.865 1.000 \n",
"b[Q2] 6934.538 1.000 \n",
"b[Q3] 5965.396 1.001 \n",
"b[Q4] 6661.555 1.000 \n",
"b[Q5] 6823.221 1.001 \n",
"b[Q6] 6160.830 1.000 \n",
"b[Q7] 5786.015 1.000 \n",
"b[Q8] 7028.130 1.000 \n",
"b[Q9] 4611.401 1.001 \n",
"b[Q10] 6487.782 1.000 \n",
"b[Q11] 5415.573 1.000 \n",
"b[Q12] 5278.103 1.001 \n",
"b[Q13] 4709.581 1.001 \n",
"b[Q14] 5899.100 1.002 \n",
"b[Q15] 4840.131 1.000 \n",
"b[Q16] 5919.292 1.000 \n",
"b[Q17] 5144.036 1.000 \n",
"b[Q18] 4587.220 1.000 \n",
"b[Q19] 5383.458 1.001 \n",
"b[Q20] 4649.440 1.000 "
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# パラメータb\n",
"summary['b[Q1]':'b[Q20]']"
]
},
{
"cell_type": "markdown",
"id": "fa0b7fb2-76ba-4449-8694-3c54a792c6c2",
"metadata": {},
"source": [
"期待値を見ると，最も難易度が高いと推定される設問はQ3．時点はQ13，Q17あたり．逆に最も簡単と推定されたのはQ1，Q4，Q19あたりか．"
]
},
{
"cell_type": "markdown",
"id": "80a3a235-7bdf-405c-b75d-032e8cc670b6",
"metadata": {},
"source": [
"答え合わせのため，下のパラメータを格納した`item_df`を表示．"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "8bb1dd4b-709a-4b2a-be49-ab7646cc8c2e",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>b</th>\n",
" <th>c</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>Q1</th>\n",
" <td>0.530947</td>\n",
" <td>-1.126938</td>\n",
" <td>0.262671</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q2</th>\n",
" <td>0.778774</td>\n",
" <td>0.440435</td>\n",
" <td>0.016508</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q3</th>\n",
" <td>0.851570</td>\n",
" <td>1.997521</td>\n",
" <td>0.078272</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q4</th>\n",
" <td>0.496084</td>\n",
" <td>-1.853463</td>\n",
" <td>0.310627</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q5</th>\n",
" <td>0.592267</td>\n",
" <td>-0.531258</td>\n",
" <td>0.134943</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q6</th>\n",
" <td>0.827886</td>\n",
" <td>1.180413</td>\n",
" <td>0.169449</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q7</th>\n",
" <td>0.686584</td>\n",
" <td>0.559992</td>\n",
" <td>0.363648</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q8</th>\n",
" <td>0.735757</td>\n",
" <td>1.890666</td>\n",
" <td>0.340324</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q9</th>\n",
" <td>0.366514</td>\n",
" <td>0.357466</td>\n",
" <td>0.198492</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q10</th>\n",
" <td>0.504268</td>\n",
" <td>-1.591176</td>\n",
" <td>0.140980</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q11</th>\n",
" <td>0.800319</td>\n",
" <td>1.106090</td>\n",
" <td>0.331420</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q12</th>\n",
" <td>0.704560</td>\n",
" <td>-1.174947</td>\n",
" <td>0.192516</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q13</th>\n",
" <td>0.927804</td>\n",
" <td>1.391537</td>\n",
" <td>0.204997</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q14</th>\n",
" <td>0.744368</td>\n",
" <td>-0.900265</td>\n",
" <td>0.013639</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q15</th>\n",
" <td>0.464219</td>\n",
" <td>1.375144</td>\n",
" <td>0.398396</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q16</th>\n",
" <td>0.433757</td>\n",
" <td>0.358018</td>\n",
" <td>0.399029</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q17</th>\n",
" <td>0.598347</td>\n",
" <td>1.633047</td>\n",
" <td>0.206599</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q18</th>\n",
" <td>0.675967</td>\n",
" <td>0.070730</td>\n",
" <td>0.262606</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q19</th>\n",
" <td>0.431652</td>\n",
" <td>-1.348594</td>\n",
" <td>0.368219</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Q20</th>\n",
" <td>0.537283</td>\n",
" <td>-0.668843</td>\n",
" <td>0.251276</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a b c\n",
"Q1 0.530947 -1.126938 0.262671\n",
"Q2 0.778774 0.440435 0.016508\n",
"Q3 0.851570 1.997521 0.078272\n",
"Q4 0.496084 -1.853463 0.310627\n",
"Q5 0.592267 -0.531258 0.134943\n",
"Q6 0.827886 1.180413 0.169449\n",
"Q7 0.686584 0.559992 0.363648\n",
"Q8 0.735757 1.890666 0.340324\n",
"Q9 0.366514 0.357466 0.198492\n",
"Q10 0.504268 -1.591176 0.140980\n",
"Q11 0.800319 1.106090 0.331420\n",
"Q12 0.704560 -1.174947 0.192516\n",
"Q13 0.927804 1.391537 0.204997\n",
"Q14 0.744368 -0.900265 0.013639\n",
"Q15 0.464219 1.375144 0.398396\n",
"Q16 0.433757 0.358018 0.399029\n",
"Q17 0.598347 1.633047 0.206599\n",
"Q18 0.675967 0.070730 0.262606\n",
"Q19 0.431652 -1.348594 0.368219\n",
"Q20 0.537283 -0.668843 0.251276"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"item_df"
]
},
{
"cell_type": "markdown",
"id": "c14f22cd-be5d-439d-a4a1-98e779f8eb35",
"metadata": {},
"source": [
"推定結果と比較すると，絶対値は外れているものの，大小関係は大体合っている？\n",
"\n",
"推定がまともかどうか調べるために，設問の真の難易度と推定された設問の難易度（期待値），回答者の真の能力値と推定された能力値の相関を調べてみる．"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "dc8198c7-4331-4aa8-b678-5bc0e59e4849",
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import pearsonr"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "b31431b8-f38a-4507-a27b-05759bcccfcf",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"coef for question difficulty: 0.8104923907592895\n"
]
}
],
"source": [
"# 設問の真の難易度\n",
"b_real = item_df.b.values\n",
"\n",
"# 設問の難易度の推定値\n",
"b_predicted = summary['b[Q1]':'b[Q20]']['mean'].values\n",
"\n",
"b_coef, _ = pearsonr(b_real, b_predicted)\n",
"print('coef for question difficulty:', b_coef)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "73ff3349-96c2-4a46-9697-2309c0fd8208",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"coef for participant capability: 0.695506191929421\n"
]
}
],
"source": [
"# 設問の真の難易度\n",
"theta_real = participant_df.theta.values\n",
"\n",
"# 設問の難易度の推定値\n",
"theta_predicted = summary['theta[0]':'theta[1999]']['mean'].values\n",
"\n",
"theta_coef, _ = pearsonr(theta_real, theta_predicted)\n",
"print('coef for participant capability:', theta_coef)"
]
},
{
"cell_type": "markdown",
"id": "89ea3c5f-7807-4800-965d-cf9d797dd99e",
"metadata": {},
"source": [
"問題の難易度も回答者の能力も，真の値と中程度以上の相関があるようだ．まあまあ良いモデリングができているのではないだろうか（PyMCの使い方が合っているかやや不安…）．"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment