Hierarchical Bayesian Modeling in logistic regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import mstats\n",
    "import pystan\n",
    "import os\n",
    "import data_import as di\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved sessions loaded from: sessions.pkl\n",
      "Loading new data files..\n",
      "Importing data file: WT2-2020-01-07-110105.txt\n",
      "Importing data file: WT2-2020-01-08-110007.txt\n",
      "Importing data file: WT2-2020-01-09-110001.txt\n",
      "Importing data file: WT2-2020-01-10-110845.txt\n",
      "Importing data file: WT2-2020-01-11-110007.txt\n",
      "Importing data file: WT2-2020-01-13-110447.txt\n",
      "Importing data file: WT2-2020-01-14-110003.txt\n",
      "Importing data file: WT2-2020-01-15-105112.txt\n",
      "Importing data file: WT2-2020-01-16-104853.txt\n",
      "Importing data file: WT2-2020-01-17-110702.txt\n",
      "Importing data file: WT2-2020-01-18-110001.txt\n",
      "Loading new data files..\n",
      "Importing data file: WT3-2020-01-02-130009.txt\n",
      "Importing data file: WT3-2020-01-03-130002.txt\n",
      "Importing data file: WT3-2020-01-04-130001.txt\n",
      "Importing data file: WT3-2020-01-06-130611.txt\n",
      "Importing data file: WT3-2020-01-07-130009.txt\n",
      "Importing data file: WT3-2020-01-08-140707.txt\n",
      "Importing data file: WT3-2020-01-09-130829.txt\n",
      "Importing data file: WT3-2020-01-10-131002.txt\n",
      "Importing data file: WT3-2020-01-11-130000.txt\n",
      "Importing data file: WT3-2020-01-13-130015.txt\n",
      "Importing data file: WT3-2020-01-14-130800.txt\n",
      "Importing data file: WT3-2020-01-15-124334.txt\n",
      "Importing data file: WT3-2020-01-16-124533.txt\n",
      "Importing data file: WT3-2020-01-17-125659.txt\n",
      "Importing data file: WT3-2020-01-18-125733.txt\n",
      "Loading new data files..\n",
      "Importing data file: WT4-2020-03-31-091700.txt\n",
      "Importing data file: WT4-2020-04-01-090623.txt\n",
      "Importing data file: WT4-2020-04-02-093656.txt\n",
      "Importing data file: WT4-2020-04-03-091040.txt\n",
      "Importing data file: WT4-2020-04-04-092112.txt\n",
      "Importing data file: WT4-2020-04-06-092056.txt\n",
      "Importing data file: WT4-2020-04-07-090829.txt\n",
      "Importing data file: WT4-2020-04-08-091022.txt\n",
      "Importing data file: WT4-2020-04-10-090658.txt\n",
      "Importing data file: WT4-2020-04-11-090952.txt\n",
      "Importing data file: WT4-2020-04-13-090722.txt\n",
      "Importing data file: WT4-2020-04-14-091043.txt\n",
      "Importing data file: WT4-2020-04-15-085847.txt\n",
      "Importing data file: WT4-2020-04-16-090043.txt\n",
      "Importing data file: WT4-2020-04-17-092026.txt\n",
      "Loading new data files..\n",
      "Importing data file: WT5-2020-03-31-105409.txt\n",
      "Importing data file: WT5-2020-04-01-105007.txt\n",
      "Importing data file: WT5-2020-04-02-111352.txt\n",
      "Importing data file: WT5-2020-04-03-110055.txt\n",
      "Importing data file: WT5-2020-04-04-110006.txt\n",
      "Importing data file: WT5-2020-04-06-110000.txt\n",
      "Importing data file: WT5-2020-04-07-105814.txt\n",
      "Importing data file: WT5-2020-04-08-110001.txt\n",
      "Importing data file: WT5-2020-04-09-110254.txt\n",
      "Importing data file: WT5-2020-04-10-110011.txt\n",
      "Importing data file: WT5-2020-04-11-110605.txt\n",
      "Importing data file: WT5-2020-04-13-110000.txt\n",
      "Importing data file: WT5-2020-04-14-110708.txt\n",
      "Importing data file: WT5-2020-04-15-110001.txt\n",
      "Importing data file: WT5-2020-04-16-110250.txt\n",
      "Importing data file: WT5-2020-04-17-105703.txt\n",
      "Loading new data files..\n",
      "Importing data file: WT6-2020-03-31-131013.txt\n",
      "Importing data file: WT6-2020-04-01-130046.txt\n",
      "Importing data file: WT6-2020-04-02-125501.txt\n",
      "Importing data file: WT6-2020-04-03-130948.txt\n",
      "Importing data file: WT6-2020-04-04-130837.txt\n",
      "Importing data file: WT6-2020-04-06-131703.txt\n",
      "Importing data file: WT6-2020-04-07-131833.txt\n",
      "Importing data file: WT6-2020-04-08-130831.txt\n",
      "Importing data file: WT6-2020-04-09-130538.txt\n",
      "Importing data file: WT6-2020-04-10-130343.txt\n",
      "Importing data file: WT6-2020-04-11-130645.txt\n",
      "Importing data file: WT6-2020-04-13-131341.txt\n",
      "Importing data file: WT6-2020-04-14-130847.txt\n",
      "Importing data file: WT6-2020-04-15-130731.txt\n",
      "Importing data file: WT6-2020-04-16-132836.txt\n",
      "Importing data file: WT6-2020-04-17-133722.txt\n"
     ]
    }
   ],
   "source": [
    "data_dir = 'C:\\\\Users\\\\Masakazu TAIRA\\\\Desktop\\\\Two_step_data\\\\'\n",
    "experiment_WT1 = di.Experiment(data_dir+'WT1') \n",
    "experiment_WT1.save()\n",
    "WT1_sessions = experiment_WT1.get_sessions(subject_IDs='all', when=['2020-01-7',...,'2020-01-18'])\n",
    "\n",
    "\n",
    "WT2experiment = di.Experiment(data_dir+'WT2')\n",
    "WT2experiment.save()\n",
    "WT2_sessions = WT2experiment.get_sessions(subject_IDs='all', when=['2020-01-07',...,'2020-01-18'])\n",
    "\n",
    "WT3experiment = di.Experiment(data_dir+'WT3')\n",
    "WT3experiment.save()\n",
    "WT3_sessions = WT3experiment.get_sessions(subject_IDs='all', when=['2020-01-02',...,'2020-01-18'])\n",
    "\n",
    "WT4experiment = di.Experiment(data_dir+'WT4')\n",
    "WT4experiment.save()\n",
    "WT4_sessions = WT4experiment.get_sessions(subject_IDs='all', when=['2020-03-31',...,'2020-04-17'])\n",
    "\n",
    "WT5experiment = di.Experiment(data_dir+'WT5')\n",
    "WT5experiment.save()\n",
    "WT5_sessions = WT5experiment.get_sessions(subject_IDs='all', when=['2020-03-31',...,'2020-04-17'])\n",
    "\n",
    "WT6experiment = di.Experiment(data_dir+'WT6')\n",
    "WT6experiment.save()\n",
    "WT6_sessions = WT6experiment.get_sessions(subject_IDs='all', when=['2020-03-31',...,'2020-04-17'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[11, 11, 15, 15, 16, 16]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_sessions = [len(WT1_sessions),len(WT2_sessions),len(WT3_sessions),len(WT4_sessions),len(WT5_sessions),len(WT6_sessions)]\n",
    "n_sessions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sessioninfo(session, session_type = 0):\n",
    "    # if session_type = 1, the session did not give photostimulation\n",
    "    b = session.print_lines\n",
    "    trial = b[0].split(\" \")\n",
    "    n_trial    = len(b)\n",
    "\n",
    "    trans_state = int(trial[8][3] == 'A') # 0 = B, 1 = A\n",
    "    choice     = np.zeros(n_trial, dtype = int) \n",
    "    second_state = np.zeros(n_trial, dtype = int) \n",
    "    outcome = np.zeros(n_trial, dtype = int) \n",
    "    exp_mov    = np.zeros(n_trial, dtype = int)\n",
    "    trial_type = np.zeros(n_trial, dtype = int)\n",
    "    stim       = np.zeros(n_trial, dtype = int) \n",
    "    reward_state = np.zeros(n_trial, dtype = int)\n",
    "\n",
    "\n",
    "    for i in range(n_trial):\n",
    "        d = b[i].split(\" \")\n",
    "        if d[4] == 'C:1': # 0 = right, 1 = left\n",
    "            choice[i] = 1\n",
    "        if d[5] == 'S:1': # 0 = down, 1 = up \n",
    "            second_state[i] = 1\n",
    "        if d[6] == 'O:1': # 0 = unrewarded, 1 = rewarded\n",
    "            outcome[i] = 1        \n",
    "        exp_mov[i] = float(d[7][3:8])\n",
    "        if d[8][2] == 'D':\n",
    "            reward_state[i] = 0\n",
    "        elif d[8][2] == 'U':\n",
    "            reward_state[i] = 1\n",
    "        else:\n",
    "            reward_state[i] = 2\n",
    "        if d[9] == 'CT:FC': # 1 = free choice, 0 = forced choice\n",
    "            trial_type[i] = 1\n",
    "        if session_type  == 1:\n",
    "            if d[11] == 'STIM:Y': # 0 = blue light, 1 = yellow light\n",
    "                stim[i] = 1\n",
    "    \n",
    "    block_info = {'block_type':[], 'start':[0], 'end':[], 'block_length':[]}\n",
    "    block_counter = 0\n",
    "\n",
    "    for i in range (n_trial-1):\n",
    "        rew_state = reward_state[i]\n",
    "        if rew_state != reward_state[i+1]:\n",
    "            if rew_state == 1:\n",
    "                block_info['block_type'].append('Up')\n",
    "            elif rew_state == 2:\n",
    "                block_info['block_type'].append('Neutral')\n",
    "            else:    \n",
    "                block_info['block_type'].append('Down')\n",
    "\n",
    "            block_info['end'].append(i)\n",
    "            block_info['start'].append(i+1)\n",
    "            length = (block_info['end'][block_counter] - block_info['start'][block_counter]) + 1         \n",
    "            block_info['block_length'].append(length)\n",
    "            block_counter += 1\n",
    "        if i == n_trial-2:\n",
    "            block_info['end'].append(n_trial-1)\n",
    "            length = (block_info['end'][block_counter] - block_info['start'][block_counter]) + 1\n",
    "            block_info['block_length'].append(length)\n",
    "            if rew_state == 1:\n",
    "                block_info['block_type'].append('Up')\n",
    "            elif rew_state == 2:\n",
    "                block_info['block_type'].append('Neutral')\n",
    "            else:    \n",
    "                block_info['block_type'].append('Down')    \n",
    "    \n",
    "    \n",
    "    \n",
    "    return n_trial, choice, second_state, outcome, trial_type, trans_state, reward_state, stim, block_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GetPredictors(session, session_type = 0):\n",
    "    n_trial, choice, second_state, outcome, trial_type, trans_state, reward_state, stim, block_info = sessioninfo(session, session_type)\n",
    "    predictors = {'p_correct':[], 'p_choice':[], 'p_outcome':[], 'p_trans':[], 'p_transxoutcome':[]}\n",
    "    for i in range (n_trial):\n",
    "        if reward_state[i] == 2: #neutral block\n",
    "            rew_state = 1\n",
    "        elif reward_state[i] == 1: # up reward block       \n",
    "            rew_state = 2\n",
    "        else: rew_state = reward_state[i]\n",
    "        correct = 0.5 * (rew_state - 1) * (2 * trans_state - 1)\n",
    "        p_outcome = (outcome[i] == choice[i]) - 0.5\n",
    "        #p_trialtype = (trial_type[i]  == choice[i]) - 0.5\n",
    "        if trans_state == 1:\n",
    "            trans = int(choice[i] == second_state[i])             \n",
    "        else:\n",
    "            trans = 1 - int(choice[i] == second_state[i])\n",
    "            \n",
    "        interaction = (trans == outcome[i]).astype(int)\n",
    "        \n",
    "        predictors['p_choice'].append(choice[i]-0.5)\n",
    "\n",
    "        predictors['p_outcome'].append(p_outcome)\n",
    "        #predictors['p_trialtype'].append(p_trialtype)\n",
    "        predictors['p_correct'].append(correct)\n",
    "        predictors['p_trans'].append((trans == choice[i]) - 0.5)\n",
    "        predictors['p_transxoutcome'].append((interaction== choice[i]) - 0.5)       \n",
    "\n",
    "    a = pd.DataFrame.from_dict(predictors)\n",
    "    #print(predictors)\n",
    "    trial_type = pd.DataFrame(trial_type, columns = ['trial_type'])\n",
    "    target = (choice[1:]).astype(int)\n",
    "    b = pd.DataFrame(target,columns =['t_choices'])\n",
    "\n",
    "    data = pd.concat([a, b], axis=1, join = 'inner')\n",
    "    \n",
    "    return data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 10):\n",
    "    \n",
    "    data_selected = data\n",
    "    \n",
    "    df = pd.DataFrame(trial_type, columns = ['trial_type'])\n",
    "    \n",
    "    \n",
    "    if type == 'exr':\n",
    "        for value in block_info['start']:\n",
    "            for j in range(n_selection):\n",
    "                if value+j < len(data):\n",
    "                    data_selected = data_selected.drop(value+j)\n",
    "                    df = df.drop(value+j)\n",
    "                elif value+j == len(data):\n",
    "                    df = df.drop(value+j)\n",
    "            \n",
    "        index = df[df['trial_type'] == 0].index\n",
    "#         print(data_selected)\n",
    "#         print(index)\n",
    "\n",
    "        for i in (index):\n",
    "            if ((i-n_selection) in block_info['start']) == False: # if a trial before forced choice trials are already excluded, nothing will happen\n",
    "                data_selected = data_selected.drop(i-1)    \n",
    "\n",
    "    if type == 'All':\n",
    "        index = df[df['trial_type'] == 0].index     \n",
    "        for i in (index):\n",
    "            if i != 0:\n",
    "                data_selected = data_selected.drop(i-1)\n",
    "        \n",
    "    return data_selected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.DataFrame([])\n",
    "df2 = pd.DataFrame([])\n",
    "df3 = pd.DataFrame([])\n",
    "df4 = pd.DataFrame([])\n",
    "df5 = pd.DataFrame([])\n",
    "df6 = pd.DataFrame([])\n",
    "for i in range (n_sessions[0]):\n",
    "    n_trial, _,_ ,_ ,trial_type ,_ ,_ ,_ , block_info = sessioninfo(WT1_sessions[i],session_type=0)\n",
    "    data = GetPredictors(WT1_sessions[i], session_type = 0)\n",
    "    data_analyzed = selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 5)\n",
    "    data_analyzed[\"session_id\"] = i+1\n",
    "    df1 = df1.append(data_analyzed, ignore_index=True)\n",
    "for i in range (n_sessions[1]):\n",
    "    n_trial, _,_ ,_ ,trial_type ,_ ,_ ,_ , block_info = sessioninfo(WT2_sessions[i],session_type=0)\n",
    "    data = GetPredictors(WT2_sessions[i], session_type = 0)\n",
    "    data_analyzed = selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 5)\n",
    "    data_analyzed[\"session_id\"] = i+1\n",
    "    df2 = df2.append(data_analyzed, ignore_index=True)\n",
    "for i in range (n_sessions[2]):\n",
    "    n_trial, _,_ ,_ ,trial_type ,_ ,_ ,_ , block_info = sessioninfo(WT3_sessions[i],session_type=0)\n",
    "    data = GetPredictors(WT3_sessions[i], session_type = 0)\n",
    "    data_analyzed = selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 5)\n",
    "    data_analyzed[\"session_id\"] = i+1\n",
    "    df3 = df3.append(data_analyzed, ignore_index=True)    \n",
    "for i in range (n_sessions[3]):\n",
    "    n_trial, _,_ ,_ ,trial_type ,_ ,_ ,_ , block_info = sessioninfo(WT4_sessions[i],session_type=0)\n",
    "    data = GetPredictors(WT4_sessions[i], session_type = 0)\n",
    "    data_analyzed = selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 5)\n",
    "    data_analyzed[\"session_id\"] = i+1\n",
    "    df4 = df4.append(data_analyzed, ignore_index=True)    \n",
    "for i in range (n_sessions[4]):\n",
    "    n_trial, _,_ ,_ ,trial_type ,_ ,_ ,_ , block_info = sessioninfo(WT5_sessions[i],session_type=0)\n",
    "    data = GetPredictors(WT5_sessions[i], session_type = 0)\n",
    "    data_analyzed = selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 5)\n",
    "    data_analyzed[\"session_id\"] = i+1\n",
    "    df5 = df5.append(data_analyzed, ignore_index=True)\n",
    "for i in range (n_sessions[5]):\n",
    "    n_trial, _,_ ,_ ,trial_type ,_ ,_ ,_ , block_info = sessioninfo(WT6_sessions[i],session_type=0)\n",
    "    data = GetPredictors(WT6_sessions[i], session_type = 0)\n",
    "    data_analyzed = selectionoftrials(data, block_info, trial_type, type = 'exr',n_selection = 5)\n",
    "    data_analyzed[\"session_id\"] = i+1\n",
    "    df6 = df6.append(data_analyzed, ignore_index=True)    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "stan_model_CoChOTI = \"\"\"\n",
    "data {\n",
    "  int N;\n",
    "  int N_id;\n",
    "  real Ch[N];\n",
    "  real Co[N];\n",
    "  real O[N];\n",
    "  real T[N];\n",
    "  real I[N];\n",
    "  int Y[N];\n",
    "  int<lower=1, upper=N_id> s_id[N];\n",
    "}\n",
    "\n",
    "parameters{\n",
    "  real a10;\n",
    "  real a20;\n",
    "  real a30;\n",
    "  real a40;\n",
    "  real a50;\n",
    "  real b0;\n",
    "  real a1_id[N_id];\n",
    "  real a2_id[N_id];\n",
    "  real a3_id[N_id];\n",
    "  real a4_id[N_id];\n",
    "  real a5_id[N_id];\n",
    "  real b_id[N_id];\n",
    "  real<lower=0> s_a1;\n",
    "  real<lower=0> s_a2;\n",
    "  real<lower=0> s_a3;\n",
    "  real<lower=0> s_a4;\n",
    "  real<lower=0> s_a5;\n",
    "  \n",
    "  real<lower=0> s_b;\n",
    "}\n",
    "\n",
    "transformed parameters{\n",
    "  real a1[N_id];\n",
    "  real a2[N_id];\n",
    "  real a3[N_id];\n",
    "  real a4[N_id];\n",
    "  real a5[N_id];\n",
    "  real b[N_id];\n",
    "  for (n in 1:N_id){\n",
    "      a1[n] = a10 + a1_id[n];\n",
    "      a2[n] = a20 + a2_id[n];\n",
    "      a3[n] = a30 + a3_id[n];\n",
    "      a4[n] = a40 + a4_id[n];\n",
    "      a5[n] = a50 + a5_id[n];\n",
    "      b[n] = b0 + b_id[n];\n",
    "  }\n",
    "\n",
    "}\n",
    "\n",
    "model {\n",
    "    for (id in 1:N_id){\n",
    "      a1_id[id] ~ normal(0,s_a1);\n",
    "      a1_id[id] ~ normal(0,s_a2);\n",
    "      a1_id[id] ~ normal(0,s_a3);      \n",
    "      a1_id[id] ~ normal(0,s_a4);\n",
    "      a1_id[id] ~ normal(0,s_a5);\n",
    "      b_id[id] ~ normal(0,s_b);\n",
    "    }\n",
    "    \n",
    "    for (n in 1:N){\n",
    "      Y[n] ~ bernoulli_logit(a1[s_id[n]] * Ch[n] + a2[s_id[n]] * Co[n] + a3[s_id[n]] * O[n] + a4[s_id[n]] * T[n] + a5[s_id[n]] * I[n] + b[s_id[n]]);\n",
    "    }\n",
    "\n",
    "}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "stan_model_CoChOTINoH = \"\"\"\n",
    "data {\n",
    "  int N;\n",
    "  real Ch[N];\n",
    "  real Co[N];\n",
    "  real O[N];\n",
    "  real T[N];\n",
    "  real I[N];\n",
    "  int Y[N];\n",
    "}\n",
    "\n",
    "parameters{\n",
    "  real a1;\n",
    "  real a2;\n",
    "  real a3;\n",
    "  real a4;\n",
    "  real a5;\n",
    "  real b;\n",
    "\n",
    "}\n",
    "\n",
    "model {\n",
    "    for (n in 1:N){\n",
    "      Y[n] ~ bernoulli_logit(a1 * Ch[n] + a2 * Co[n] + a3* O[n] + a4 * T[n] + a5 * I[n] + b);\n",
    "    }\n",
    "\n",
    "}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_24c01e7229a959f75d6909c37e60dea8 NOW.\n"
     ]
    }
   ],
   "source": [
    "sm = pystan.StanModel(model_code = stan_model_CoChOTI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1266a92c64ec34c25eba155f9758eea3 NOW.\n"
     ]
    }
   ],
   "source": [
    "smnoH = pystan.StanModel(model_code = stan_model_CoChOTINoH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "stan_data = {\"N\":df.shape[0], \"N_id\":n_sessions, \"Ch\":df[\"p_choice\"], \"Co\":df[\"p_correct\"], \"O\":df[\"p_outcome\"],\"T\":df[\"p_trans\"],\"I\":df[\"p_transxoutcome\"],\"Y\": df[\"t_choices\"], \"s_id\":df[\"session_id\"]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_24c01e7229a959f75d6909c37e60dea8.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "             mean se_mean     sd    2.5%     25%      50%     75%  97.5%  n_eff   Rhat\n",
      "a10          0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a20         -0.44    0.26   0.32   -0.79   -0.78    -0.53   -0.02  -0.02      2 114.87\n",
      "a30          0.54    1.34   1.64    -0.8   -0.79    -0.44    2.85   2.85      2 235.86\n",
      "a40          2.45    0.71   0.87    1.45    1.46     2.33    3.56   3.57      2 252.02\n",
      "a50         -0.19    2.58   3.16   -4.65   -4.64     1.64    2.41   2.43      2 652.49\n",
      "b0           0.04    0.05   0.07   -0.02   -0.02   8.7e-3    0.12   0.15      2  12.89\n",
      "a1_id[1]   2.2e-9  5.5e-9 4.2e-7 -8.8e-7 -2.4e-7   7.7e-9  2.4e-7 8.6e-7   5797    1.0\n",
      "a1_id[2]  -7.8e-9  5.3e-9 4.1e-7 -8.7e-7 -2.6e-7   1.1e-9  2.4e-7 8.1e-7   5881    1.0\n",
      "a1_id[3]   2.3e-9  5.3e-9 4.2e-7 -8.5e-7 -2.5e-7   4.8e-9  2.5e-7 8.7e-7   6197    1.0\n",
      "a1_id[4]  -8.9e-9  5.2e-9 4.1e-7 -8.4e-7 -2.6e-7  -5.6e-9  2.4e-7 8.2e-7   6127    1.0\n",
      "a1_id[5]  -8.7e-9  5.6e-9 4.3e-7 -8.9e-7 -2.6e-7  -1.0e-8  2.4e-7 8.9e-7   5960    1.0\n",
      "a1_id[6]  9.6e-10  5.4e-9 4.1e-7 -8.2e-7 -2.5e-7  -1.1e-9  2.4e-7 8.4e-7   5776    1.0\n",
      "a1_id[7]   4.7e-9  5.2e-9 4.1e-7 -8.1e-7 -2.4e-7   4.2e-9  2.5e-7 8.2e-7   6106    1.0\n",
      "a1_id[8]   1.0e-9  5.4e-9 4.2e-7 -8.5e-7 -2.4e-7 -5.4e-10  2.5e-7 8.6e-7   6005    1.0\n",
      "a1_id[9]  -6.0e-9  5.5e-9 4.2e-7 -8.8e-7 -2.5e-7  -7.6e-9  2.3e-7 8.7e-7   5774    1.0\n",
      "a1_id[10]  2.9e-9  5.3e-9 4.2e-7 -8.2e-7 -2.5e-7  -4.2e-9  2.5e-7 8.6e-7   6182    1.0\n",
      "a1_id[11]  4.6e-9  5.7e-9 4.2e-7 -8.5e-7 -2.5e-7   9.7e-9  2.5e-7 8.7e-7   5488    1.0\n",
      "a2_id[1]     0.66    0.37   0.45    0.17    0.17     0.58    1.25   1.26      2  50.12\n",
      "a2_id[2]     0.79    0.51   0.62   -0.05   -0.04     0.97    1.44   1.45      2  87.66\n",
      "a2_id[3]     1.55    0.32   0.39    0.98    1.01     1.72     1.9   1.92      2  54.45\n",
      "a2_id[4]      1.7    0.41   0.51    1.27    1.27     1.41     2.4   2.42      2   79.4\n",
      "a2_id[5]     2.13    0.38   0.47    1.54    1.55     2.15    2.69    2.7      2  93.27\n",
      "a2_id[6]     1.97    0.11   0.13    1.78     1.8     2.03    2.09    2.1      2  40.04\n",
      "a2_id[7]      1.0    0.06   0.07    0.92    0.94     0.99    1.09   1.11      2  16.84\n",
      "a2_id[8]     1.67    0.26   0.32    1.23    1.23     1.81    1.96   1.99      2  45.17\n",
      "a2_id[9]     1.42    0.13   0.16    1.21    1.22     1.44    1.59    1.6      2  12.02\n",
      "a2_id[10]    1.47    0.36   0.44    0.84    0.85     1.78    1.79    1.8      2   89.5\n",
      "a2_id[11]    0.21    0.13   0.15   -0.02    0.01     0.29    0.34   0.39      2  16.84\n",
      "a3_id[1]     0.57    1.61   1.98   -2.22   -2.19     1.67    2.24   2.29      2 244.29\n",
      "a3_id[2]    -0.91    1.35   1.66    -3.1    -3.1    -0.54    0.91   0.92      2 618.67\n",
      "a3_id[3]    -0.73    1.32   1.61   -3.01   -3.01     0.37    0.44   0.45      2 347.01\n",
      "a3_id[4]     -0.7    1.32   1.61   -2.95   -2.94     0.06    0.76   0.78      2  80.09\n",
      "a3_id[5]    -0.96    1.21   1.48   -3.08   -3.03     0.04    0.13   0.14      2 148.34\n",
      "a3_id[6]    -0.43    1.38    1.7   -2.76   -2.76     0.23    1.22   1.24      2 272.05\n",
      "a3_id[7]     -0.3    1.41   1.73   -2.74   -2.72     0.67    1.16   1.18      2 229.08\n",
      "a3_id[8]     0.17    1.44   1.76   -2.23   -2.21     0.75    1.97    2.0      2 355.62\n",
      "a3_id[9]    -0.73    1.31    1.6   -2.97   -2.96     0.03    0.72   0.75      2 327.96\n",
      "a3_id[10]   -0.63    1.72   2.11   -3.62   -3.59     0.57    1.12   1.13      2  210.9\n",
      "a3_id[11]   -0.58    0.71   0.87   -1.73   -1.71    -0.41    0.34   0.44      2   58.1\n",
      "a4_id[1]    -1.67    0.75   0.92   -2.74   -2.74    -1.79   -0.49  -0.49      2 262.14\n",
      "a4_id[2]    -0.91    0.53   0.64   -1.73   -1.71    -0.86   -0.15  -0.13      2 122.59\n",
      "a4_id[3]    -1.54     0.9    1.1   -2.76   -2.75    -1.79   -0.09  -0.07      2 192.85\n",
      "a4_id[4]    -0.68    0.71   0.87   -1.73   -1.72    -0.74    0.42   0.42      2  178.8\n",
      "a4_id[5]    -1.79     0.6   0.73   -2.74   -2.71    -1.73   -0.94  -0.92      2  70.83\n",
      "a4_id[6]    -0.46    0.57    0.7   -1.16   -1.15    -0.72     0.5   0.52      2  79.13\n",
      "a4_id[7]    -0.59    0.86   1.06    -2.1   -2.08     0.12    0.17    0.2      2 158.54\n",
      "a4_id[8]    -0.13    0.75   0.92   -1.44   -1.42     0.42    0.62   0.63      2  256.7\n",
      "a4_id[9]     -1.6    0.37   0.46   -2.22   -2.21    -1.47   -1.12   -1.1      2  63.25\n",
      "a4_id[10]   -0.76    0.66   0.81   -1.87   -1.85     -0.5    0.06   0.07      2 175.85\n",
      "a4_id[11]   -0.95    0.41   0.51   -1.51    -1.5    -1.07   -0.28  -0.26      2  83.74\n",
      "a5_id[1]     1.47    2.78    3.4   -1.12   -1.11    -0.75    6.28   6.29      2 256.43\n",
      "a5_id[2]     1.79    2.54   3.11    -1.0   -0.99     0.24    6.13   6.13      2 439.73\n",
      "a5_id[3]     1.29    2.61   3.19   -1.36   -1.33    -0.57    5.79   5.79      2 537.81\n",
      "a5_id[4]     1.56    2.63   3.23   -1.29   -1.28    -0.11    6.07   6.07      2 1228.9\n",
      "a5_id[5]     1.54    2.45   3.01   -1.39   -1.36     0.32    5.68   5.69      2 390.24\n",
      "a5_id[6]     1.71    2.56   3.14   -1.07   -1.04     0.07    6.08   6.11      2 243.49\n",
      "a5_id[7]     1.77    2.43   2.98   -0.61   -0.61    -0.04    5.97   5.98      2 588.33\n",
      "a5_id[8]     1.42    2.32   2.84   -0.63   -0.62    -0.55    5.44   5.45      2 757.65\n",
      "a5_id[9]     1.44    2.67   3.27   -1.63   -1.54    -0.09    5.99    6.0      2 203.06\n",
      "a5_id[10]    2.12    2.75   3.37   -0.74   -0.71     0.23    6.85   6.86      2 523.75\n",
      "a5_id[11]    1.93    2.06   2.53   -0.07   -0.07     0.37    5.49   5.51      2 627.05\n",
      "b_id[1]     -0.32    0.04   0.04   -0.38   -0.37     -0.3   -0.28  -0.26      2  11.97\n",
      "b_id[2]     -0.05  2.8e-3 7.8e-3   -0.06   -0.05    -0.05   -0.04  -0.03      8   1.52\n",
      "b_id[3]   -9.1e-3     0.1   0.12   -0.11    -0.1    -0.09    0.16   0.17      2  32.01\n",
      "b_id[4]      0.53    0.05   0.06    0.44    0.49     0.51    0.61   0.62      2   9.04\n",
      "b_id[5]      0.06    0.12   0.14   -0.06   -0.05    -0.03    0.26   0.27      2  23.04\n",
      "b_id[6]      0.03    0.05   0.07   -0.06   -0.06     0.06    0.08   0.12      2  10.25\n",
      "b_id[7]     -0.29    0.12   0.15   -0.51    -0.5    -0.22   -0.16  -0.15      2  35.79\n",
      "b_id[8]     -0.38     0.1   0.12   -0.51   -0.49    -0.42   -0.22  -0.21      2  15.95\n",
      "b_id[9]      0.32    0.07   0.08    0.22    0.26     0.29    0.43   0.44      2  11.01\n",
      "b_id[10]     0.19    0.02   0.02    0.15    0.16     0.19    0.21   0.21      2   5.04\n",
      "b_id[11]     0.65    0.16    0.2     0.4    0.41     0.65    0.87   0.92      2  21.86\n",
      "s_a1       2.5e-5  2.7e-5 3.4e-5  7.7e-7  9.1e-7   9.9e-7  7.1e-5 7.4e-5      2  66.87\n",
      "s_a2       2.2e-5  1.9e-5 2.3e-5  4.9e-6  5.6e-6   6.1e-6  5.4e-5 5.4e-5      2  95.07\n",
      "s_a3       4.5e-5  4.6e-5 5.7e-5  8.1e-7  1.4e-6   8.7e-6  1.2e-4 1.3e-4      2 301.42\n",
      "s_a4       6.1e-6  6.1e-6 7.5e-6  2.4e-7  4.2e-7   1.3e-6  1.7e-5 1.7e-5      2  65.26\n",
      "s_a5       7.1e-7  2.4e-7 3.0e-7  2.6e-7  4.2e-7   6.7e-7  1.0e-6 1.2e-6      2    4.1\n",
      "s_b          0.37    0.05   0.06    0.28    0.29     0.37    0.44   0.45      2  18.49\n",
      "a1[1]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[2]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[3]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[4]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[5]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[6]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[7]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[8]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[9]        0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[10]       0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a1[11]       0.84    0.12   0.15    0.62    0.63     0.92    0.97   0.98      2  45.11\n",
      "a2[1]        0.22    0.15   0.18 -2.0e-4    0.06     0.15    0.46   0.49      2  20.98\n",
      "a2[2]        0.34    0.25    0.3   -0.07   -0.06     0.44    0.66   0.67      2  35.99\n",
      "a2[3]         1.1    0.07   0.09    0.96     1.0     1.12    1.19   1.22      2  10.93\n",
      "a2[4]        1.25    0.42   0.51     0.6    0.63     1.25    1.87    1.9      2  61.55\n",
      "a2[5]        1.69    0.28   0.35    1.34    1.39     1.53    2.16   2.18      2  58.78\n",
      "a2[6]        1.53    0.36   0.44     1.0    1.01      1.5    2.07   2.09      2 102.72\n",
      "a2[7]        0.56    0.31   0.38    0.19     0.2      0.4    1.07   1.09      2 109.26\n",
      "a2[8]        1.23    0.03   0.04    1.16    1.19     1.21    1.27   1.31      2   5.46\n",
      "a2[9]        0.98    0.18   0.22    0.63     0.7     1.07    1.19   1.21      2  15.55\n",
      "a2[10]       1.03    0.14   0.17    0.82    0.83      1.0    1.24   1.26      2  32.38\n",
      "a2[11]      -0.23    0.29   0.36   -0.55   -0.51    -0.44    0.27   0.28      2   34.1\n",
      "a3[1]        1.11    0.28   0.34    0.62    0.66     1.23    1.45   1.49      2  34.11\n",
      "a3[2]       -0.37    0.37   0.46   -1.01   -0.96    -0.25    0.11   0.12      2  73.05\n",
      "a3[3]       -0.19    0.14   0.18   -0.44   -0.42    -0.16 -9.3e-3   0.03      2  20.84\n",
      "a3[4]       -0.17    0.36   0.44   -0.79    -0.7     -0.1    0.32   0.37      2  19.86\n",
      "a3[5]       -0.42    0.16   0.19   -0.67   -0.66     -0.4   -0.21  -0.17      2  15.61\n",
      "a3[6]        0.11    0.21   0.26   -0.25   -0.19     0.09    0.42   0.45      2  23.41\n",
      "a3[7]        0.24    0.08    0.1    0.11    0.13     0.23    0.36   0.39      2   10.8\n",
      "a3[8]        0.71    0.53   0.65   -0.06   -0.05     0.63    1.53   1.56      2  82.95\n",
      "a3[9]        -0.2    0.12   0.15   -0.44   -0.38    -0.11   -0.08  -0.05      2  15.52\n",
      "a3[10]      -0.09    0.39   0.47   -0.77   -0.74     0.14    0.33   0.34      2  59.24\n",
      "a3[11]      -0.04    0.69   0.85   -0.88   -0.82    -0.41    1.12   1.14      2  47.61\n",
      "a4[1]        0.78    0.14   0.17    0.53    0.55     0.83    0.96   0.97      2  34.06\n",
      "a4[2]        1.54    0.18   0.22     1.3    1.32     1.47    1.84   1.85      2  44.07\n",
      "a4[3]        0.91    0.28   0.34    0.53    0.55     0.81    1.36   1.39      2   53.8\n",
      "a4[4]        1.77     0.1   0.12    1.58     1.6     1.84    1.87   1.88      2  22.06\n",
      "a4[5]        0.66    0.11   0.14    0.51    0.53      0.6    0.83   0.87      2  12.55\n",
      "a4[6]        1.99    0.27   0.33    1.58    1.62     1.95    2.41   2.42      2  35.85\n",
      "a4[7]        1.86    0.38   0.47    1.47    1.49     1.58    2.51   2.53      2  66.11\n",
      "a4[8]        2.33    0.38   0.46    1.87    1.88     2.14    2.95   2.97      2 102.95\n",
      "a4[9]        0.85    0.34   0.41    0.33    0.35     0.86    1.35   1.36      2  52.51\n",
      "a4[10]       1.69    0.11   0.13    1.51    1.53     1.71    1.83   1.84      2  22.39\n",
      "a4[11]        1.5    0.32    0.4    1.15    1.19     1.26    2.06   2.07      2  55.58\n",
      "a5[1]        1.28    0.25    0.3    0.87    0.92     1.31    1.63   1.65      2  24.92\n",
      "a5[2]         1.6    0.17   0.21    1.41    1.44     1.49    1.88   1.91      2  25.95\n",
      "a5[3]         1.1    0.03   0.03    1.06    1.07     1.08    1.14   1.16      2   5.04\n",
      "a5[4]        1.37    0.14   0.17    1.12    1.14     1.43    1.53   1.55      2  30.85\n",
      "a5[5]        1.35    0.36   0.44    1.03    1.04     1.05    1.94   2.01      2  44.72\n",
      "a5[6]        1.51    0.12   0.15    1.35    1.37     1.45     1.7   1.78      2  11.06\n",
      "a5[7]        1.58    0.16   0.19    1.32    1.34      1.6     1.8   1.82      2  25.64\n",
      "a5[8]        1.23    0.38   0.46    0.79    0.81     1.02    1.87   1.88      2  77.14\n",
      "a5[9]        1.25    0.24    0.3    0.79    0.88     1.35    1.55   1.59      2   19.4\n",
      "a5[10]       1.93    0.18   0.21    1.68    1.71     1.87     2.2   2.23      2   24.7\n",
      "a5[11]       1.74    0.52   0.64    0.84    0.86     2.01    2.35   2.36      2  96.54\n",
      "b[1]        -0.28    0.09    0.1    -0.4   -0.39     -0.3   -0.15  -0.13      2  18.09\n",
      "b[2]      -9.3e-3    0.05   0.07   -0.08   -0.07    -0.03    0.06   0.11      2   8.92\n",
      "b[3]         0.03    0.08   0.09    -0.1   -0.08     0.03    0.14   0.15      2   15.3\n",
      "b[4]         0.57    0.06   0.07    0.45     0.5     0.59    0.63   0.66      2   7.57\n",
      "b[5]          0.1    0.09   0.11   -0.04   -0.01     0.07    0.24   0.25      2  14.67\n",
      "b[6]         0.07    0.02   0.03    0.02    0.05     0.06     0.1   0.13      2   3.73\n",
      "b[7]        -0.25    0.07   0.09   -0.39   -0.36     -0.2   -0.18  -0.18      2  15.39\n",
      "b[8]        -0.34    0.06   0.08   -0.45   -0.39    -0.36   -0.24  -0.23      2   7.23\n",
      "b[9]         0.36    0.06   0.08    0.22    0.28      0.4    0.41   0.44      2   7.27\n",
      "b[10]        0.23    0.07   0.08    0.13    0.14      0.2    0.33   0.35      2  13.57\n",
      "b[11]        0.69    0.16    0.2     0.4    0.43     0.77    0.85    0.9      2  17.39\n",
      "lp__       -634.2   27.91  34.41  -684.6  -677.1   -625.4  -599.5 -591.3      2  11.26\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 14:17:08 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n",
      "WARNING:pystan:6000 of 6000 iterations saturated the maximum tree depth of 10 (100 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n"
     ]
    }
   ],
   "source": [
    "fit =sm.sampling(data = stan_data, iter = 3000, warmup = 1000, chains = 3, seed = 123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "stan_datanoH1 = {\"N\":df1.shape[0],\"Ch\":df1[\"p_choice\"], \"Co\":df1[\"p_correct\"], \"O\":df1[\"p_outcome\"],\"T\":df1[\"p_trans\"],\"I\":df1[\"p_transxoutcome\"],\"Y\": df1[\"t_choices\"]}\n",
    "stan_datanoH2 = {\"N\":df2.shape[0],\"Ch\":df2[\"p_choice\"], \"Co\":df2[\"p_correct\"], \"O\":df2[\"p_outcome\"],\"T\":df2[\"p_trans\"],\"I\":df2[\"p_transxoutcome\"],\"Y\": df2[\"t_choices\"]}\n",
    "stan_datanoH3 = {\"N\":df3.shape[0],\"Ch\":df3[\"p_choice\"], \"Co\":df3[\"p_correct\"], \"O\":df3[\"p_outcome\"],\"T\":df3[\"p_trans\"],\"I\":df3[\"p_transxoutcome\"],\"Y\": df3[\"t_choices\"]}\n",
    "stan_datanoH4 = {\"N\":df4.shape[0],\"Ch\":df4[\"p_choice\"], \"Co\":df4[\"p_correct\"], \"O\":df4[\"p_outcome\"],\"T\":df4[\"p_trans\"],\"I\":df4[\"p_transxoutcome\"],\"Y\": df4[\"t_choices\"]}\n",
    "stan_datanoH5 = {\"N\":df5.shape[0],\"Ch\":df5[\"p_choice\"], \"Co\":df5[\"p_correct\"], \"O\":df5[\"p_outcome\"],\"T\":df5[\"p_trans\"],\"I\":df5[\"p_transxoutcome\"],\"Y\": df5[\"t_choices\"]}\n",
    "stan_datanoH6 = {\"N\":df6.shape[0],\"Ch\":df6[\"p_choice\"], \"Co\":df6[\"p_correct\"], \"O\":df6[\"p_outcome\"],\"T\":df6[\"p_trans\"],\"I\":df6[\"p_transxoutcome\"],\"Y\": df6[\"t_choices\"]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_noH1 = smnoH.sampling(data = stan_datanoH1, iter = 3000, warmup = 1000, chains =3, seed = 123)\n",
    "fit_noH2 = smnoH.sampling(data = stan_datanoH2, iter = 3000, warmup = 1000, chains =3, seed = 123)\n",
    "fit_noH3 = smnoH.sampling(data = stan_datanoH3, iter = 3000, warmup = 1000, chains =3, seed = 123)\n",
    "fit_noH4 = smnoH.sampling(data = stan_datanoH4, iter = 3000, warmup = 1000, chains =3, seed = 123)\n",
    "fit_noH5 = smnoH.sampling(data = stan_datanoH5, iter = 3000, warmup = 1000, chains =3, seed = 123)\n",
    "fit_noH6 = smnoH.sampling(data = stan_datanoH6, iter = 3000, warmup = 1000, chains =3, seed = 123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_1266a92c64ec34c25eba155f9758eea3.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
      "a1     0.78  1.6e-3   0.12   0.55    0.7   0.78   0.86   1.02   5621    1.0\n",
      "a2     0.68  1.6e-3   0.13   0.41   0.59   0.68   0.77   0.94   6985    1.0\n",
      "a3     0.01  1.6e-3   0.12  -0.23  -0.07   0.02    0.1   0.25   5553    1.0\n",
      "a4     1.25  1.6e-3   0.12   1.01   1.17   1.25   1.33   1.49   5904    1.0\n",
      "a5     1.26  1.8e-3   0.13   1.01   1.17   1.26   1.35   1.52   5630    1.0\n",
      "b      0.13  5.5e-4   0.05   0.03   0.09   0.13   0.16   0.22   7943    1.0\n",
      "lp__  -1288    0.03   1.76  -1292  -1289  -1288  -1287  -1286   2728    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 15:20:05 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit_noH1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_1266a92c64ec34c25eba155f9758eea3.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
      "a1     0.66  1.5e-3   0.11   0.44   0.58   0.66   0.73   0.87   5544    1.0\n",
      "a2     1.23  1.5e-3   0.12   0.99   1.15   1.23   1.31   1.47   6421    1.0\n",
      "a3     -0.1  1.5e-3   0.11  -0.32  -0.18   -0.1  -0.03   0.11   5596    1.0\n",
      "a4     1.48  1.4e-3   0.11   1.27   1.41   1.48   1.55   1.69   5646    1.0\n",
      "a5     1.16  1.5e-3   0.11   0.94   1.09   1.16   1.24   1.38   5466    1.0\n",
      "b      0.24  5.2e-4   0.04   0.15   0.21   0.24   0.27   0.33   7290    1.0\n",
      "lp__  -1606    0.03   1.78  -1610  -1607  -1605  -1604  -1603   2859    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 15:20:32 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit_noH2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_1266a92c64ec34c25eba155f9758eea3.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
      "a1     0.84  1.1e-3   0.09   0.66   0.78   0.84    0.9   1.01   6256    1.0\n",
      "a2     0.81  1.2e-3    0.1   0.62   0.74   0.81   0.88    1.0   6605    1.0\n",
      "a3    -0.12  1.1e-3   0.09   -0.3  -0.18  -0.12  -0.06   0.05   6289    1.0\n",
      "a4     0.87  1.1e-3   0.09    0.7   0.81   0.87   0.93   1.04   6178    1.0\n",
      "a5     0.47  1.4e-3    0.1   0.28    0.4   0.47   0.54   0.67   5467    1.0\n",
      "b      0.03  4.3e-4   0.04  -0.04 7.9e-3   0.03   0.06    0.1   7232    1.0\n",
      "lp__  -2131    0.03   1.74  -2136  -2132  -2131  -2130  -2129   3128    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 15:21:00 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit_noH3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_1266a92c64ec34c25eba155f9758eea3.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
      "a1     0.98  1.1e-3   0.09   0.81   0.92   0.98   1.04   1.14   5861    1.0\n",
      "a2     1.23  1.2e-3   0.09   1.05   1.16   1.23   1.29   1.41   6174    1.0\n",
      "a3     0.18  1.0e-3   0.08   0.01   0.12   0.18   0.23   0.34   6544    1.0\n",
      "a4     0.98  1.0e-3   0.08   0.81   0.92   0.98   1.04   1.15   6995    1.0\n",
      "a5     0.71  1.2e-3   0.09   0.53   0.64    0.7   0.77   0.89   6291    1.0\n",
      "b      -0.1  4.0e-4   0.04  -0.18  -0.13   -0.1  -0.08  -0.03   8597    1.0\n",
      "lp__  -2267    0.03   1.71  -2272  -2268  -2267  -2266  -2265   2557    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 15:21:31 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit_noH4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_1266a92c64ec34c25eba155f9758eea3.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
      "a1     0.97  1.0e-3   0.09   0.81   0.92   0.97   1.03   1.14   6782    1.0\n",
      "a2     0.76  1.1e-3   0.09   0.58    0.7   0.76   0.83   0.95   7595    1.0\n",
      "a3     0.26  1.1e-3   0.09    0.1   0.21   0.26   0.32   0.44   6455    1.0\n",
      "a4      0.9  1.1e-3   0.09   0.73   0.84    0.9   0.96   1.07   6101    1.0\n",
      "a5     0.64  1.2e-3   0.09   0.46   0.58   0.64    0.7   0.82   6213    1.0\n",
      "b     -0.02  4.0e-4   0.04  -0.09  -0.04  -0.02 6.8e-3   0.05   8415    1.0\n",
      "lp__  -2190    0.03   1.68  -2194  -2191  -2190  -2189  -2188   2976    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 15:22:01 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit_noH5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_1266a92c64ec34c25eba155f9758eea3.\n",
      "3 chains, each with iter=3000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=2000, total post-warmup draws=6000.\n",
      "\n",
      "       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
      "a1      1.8  1.1e-3   0.09   1.63   1.74    1.8   1.86   1.98   7070    1.0\n",
      "a2     1.34  1.4e-3   0.11   1.12   1.26   1.34   1.41   1.56   6186    1.0\n",
      "a3     0.62  1.1e-3   0.09   0.44   0.55   0.62   0.68    0.8   6473    1.0\n",
      "a4      1.0  1.1e-3   0.09   0.83   0.94    1.0   1.06   1.18   6373    1.0\n",
      "a5     0.41  1.4e-3    0.1   0.21   0.34   0.41   0.48   0.62   5464    1.0\n",
      "b      0.29  4.4e-4   0.04   0.21   0.26   0.29   0.32   0.37   8519    1.0\n",
      "lp__  -1920    0.03   1.76  -1924  -1921  -1920  -1919  -1918   2797    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Sun Mar 28 15:22:33 2021.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(fit_noH6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}