<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Testing within-subject models
Issues need to be solved
Step1:  Let's first replicate the tutorial on HDDM website
http
Step2:  Creating simulated data for  the experiment
set number of subjects and number of trials per level for the simulated experiment
Step3:  Next, set up parameters of the DDM for three levels and the first stimulus.
Step4:  Next, set up parameters for the second stimulus, where v is the same as for the first stimulus. This is different for z. In particular
Step5:  Now, let's generate data for both stimuli
Step6:  Setting up the HDDM regression model
The parameter z is bound between 0 and 1, but the standard linear regression does not generate values between 0 and 1. Therefore, we use a link-function, here the inverse logit $1/(1+exp(-x))$, which transform values between plus and minus infinity into values ranging from (just above) 0 to (nearly) 1. [this link function is related to logistic regression].
Next, we need to ensure that bias is $z$ for one stimulus and $1-z$ for the other stimulus. To achive this, we can simply multiply the regression output for one stimulus with $-1$. This is implemented here by dot-multiplying the regression output "x" (which is an array) with equally sized array "stim", which is 1 for all stimulus A trials and -1 for stimulus B trials. We use the patsy command dmatrix to generate such an array from the stimulus column of our stimulated data.
Note
* dot-multiplying is vague here, it should be Hadamard Product, see here.
Step7:  Now, we set up the regression models for z and v and also include the link functions. The relevant string here used by patsy is 1 + C(condition). This will generate a design matrix with an intercept (that's what the 1 is for) and two dummy variables for remaining levels. (The column in which the levels are coded has the default name condition)
Step8:  For v the link function is simply $x = x$, b/c no transformation is needed. (However, you could also analyze this experiment with response coded data. Then you would not stimulus code z but v and you would have to multiply the v for one condition with -1, with a link function like the one for z, but without the additional logit transform)
Step9:  Now, we can finally put the regression description for the hddm model together. The general for this is [{'model'
Step10:  The last step before running the model is to construct the complete hddm regression model by adding data etc.
Step11:  Now, we start the sampling, which takes a long time
Step12:  Comparing generative and recovered model parameters
First, let's print the model stats
Step13:  Working with real data
Now let's move to the real data
<ASSISTANT_TASK:>
Python Code:
# check which python is in use.
import sys
print('Notebook is running:', sys.executable)
# further check your python version
from platform import python_version
print('The current HDDM version is', python_version())
# If you are sure that conda is installed, also check the package that install
#!conda list  # list the conda
import hddm, IPython
import numpy as np
import pandas as pd
print('The current HDDM version is', hddm.__version__) # 0.8.0
# Warning:`IPython.parallel` package has been deprecated since IPython 4.0. 
print('The current IPython version is', IPython.__version__) 
print('The current Numpy version is', np.__version__) 
print('The current Pandas version is', pd.__version__)
%matplotlib inline
# Preparation
import os, time, csv
from datetime import date
import random
import kabuki, hddm
from patsy import dmatrix
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
Explanation: Testing within-subject models
Issues need to be solved:
* Add more than $v$ to the model;
* Stim-coding with HDDMRegression (there was a tutorial on this)
* Flexible design
End of explanation
import sys
# the following will save all the stdout to the txt file, which is not necessary
# sys.stdout = open('ModelRecoveryOutput.txt', 'w')
Explanation: Let's first replicate the tutorial on HDDM website
http://ski.clps.brown.edu/hddm_docs/tutorial_regression_stimcoding.html
End of explanation
n_subjects = 10
trials_per_level = 150
Explanation: Creating simulated data for  the experiment
set number of subjects and number of trials per level for the simulated experiment
End of explanation
level1a = {'v':.3, 'a':2, 't':.3, 'sv':0, 'z':.5, 'sz':0, 'st':0}
level2a = {'v':.4, 'a':2, 't':.3, 'sv':0, 'z':.6, 'sz':0, 'st':0}
level3a = {'v':.5, 'a':2, 't':.3, 'sv':0, 'z':.7, 'sz':0, 'st':0}
Explanation: Next, set up parameters of the DDM for three levels and the first stimulus.
End of explanation
level1b = {'v':.3, 'a':2, 't':.3, 'sv':0, 'z':.5, 'sz':0, 'st':0}
level2b = {'v':.4, 'a':2, 't':.3, 'sv':0, 'z':.4, 'sz':0, 'st':0}
level3b = {'v':.5, 'a':2, 't':.3, 'sv':0, 'z':.3, 'sz':0, 'st':0}
Explanation: Next, set up parameters for the second stimulus, where v is the same as for the first stimulus. This is different for z. In particular: z(stimulus_b) = 1 - z(stimulus_a).
End of explanation
random.seed(123)
np.random.seed(123)
data_a, params_a = hddm.generate.gen_rand_data({'level1': level1a,
                                                'level2': level2a,
                                                'level3': level3a},
                                              size=trials_per_level,
                                              subjs=n_subjects)
data_b, params_b = hddm.generate.gen_rand_data({'level1': level1b,
                                                'level2': level2b,
                                                'level3': level3b},
                                              size=trials_per_level,
                                              subjs=n_subjects)
data_a['stimulus'] = pd.Series(np.ones((len(data_a))), index=data_a.index)
data_b['stimulus'] = pd.Series(np.ones((len(data_b)))*2, index=data_a.index)
mydata = data_a.append(data_b,ignore_index=True)
mydata.head(10)
Explanation: Now, let's generate data for both stimuli
End of explanation
tmpdata = mydata[mydata['subj_idx'] <= 0]
tmpdata
stim = (dmatrix('0 + C(s, [[1], [-1]])', 
                {'s':tmpdata.stimulus.loc[tmpdata.index]},
                return_type="dataframe")
        )
#tmp = np.multiply(tmpdata[['stimulus']], stim[:])
tmp = np.multiply(tmpdata.stimulus.loc[tmpdata.index].to_frame(),stim)
#1/(1+np.exp(-((np.asarray(mydata.stimulus.iloc[mydata.index])) * stim)))
type(stim)
#stim.index
#type(mydata.stimulus.iloc[mydata.index])
tmp
tmp2 = tmpdata.stimulus.loc[tmpdata.index]
tmp2.to_frame().shape
# x is a pandas Series, it should be converted to dataframe to do the multiplication.
def z_link_func(x, data=mydata):
    stim = (dmatrix('0 + C(s, [[1], [-1]])',
                    {'s':data.stimulus.loc[x.index]},
                    return_type="dataframe")
#    stim = (np.asarray(dmatrix('0 + C(s, [[1], [-1]])',
                               #{'s':data.stimulus.ix[x.index]})) # original .ix is deprecated.
           )
    #print(x.shape)
    return 1/(1+np.exp(-np.multiply(x.to_frame(), stim)))
    #return 1 /( 1+np.exp(-(x * stim)))
Explanation: Setting up the HDDM regression model
The parameter z is bound between 0 and 1, but the standard linear regression does not generate values between 0 and 1. Therefore, we use a link-function, here the inverse logit $1/(1+exp(-x))$, which transform values between plus and minus infinity into values ranging from (just above) 0 to (nearly) 1. [this link function is related to logistic regression].
Next, we need to ensure that bias is $z$ for one stimulus and $1-z$ for the other stimulus. To achive this, we can simply multiply the regression output for one stimulus with $-1$. This is implemented here by dot-multiplying the regression output "x" (which is an array) with equally sized array "stim", which is 1 for all stimulus A trials and -1 for stimulus B trials. We use the patsy command dmatrix to generate such an array from the stimulus column of our stimulated data.
Note
* dot-multiplying is vague here, it should be Hadamard Product, see here.
End of explanation
z_reg = {'model': 'z ~ 1 + C(condition)', 'link_func':z_link_func}
Explanation: Now, we set up the regression models for z and v and also include the link functions. The relevant string here used by patsy is 1 + C(condition). This will generate a design matrix with an intercept (that's what the 1 is for) and two dummy variables for remaining levels. (The column in which the levels are coded has the default name condition):
End of explanation
v_reg = {'model': 'v ~ 1 + C(condition)', 'link_func': lambda x:x}
Explanation: For v the link function is simply $x = x$, b/c no transformation is needed. (However, you could also analyze this experiment with response coded data. Then you would not stimulus code z but v and you would have to multiply the v for one condition with -1, with a link function like the one for z, but without the additional logit transform):
End of explanation
reg_descr = [z_reg, v_reg]
Explanation: Now, we can finally put the regression description for the hddm model together. The general for this is [{'model': 'outcome_parameter ~ patsy_design_string', 'link_func': your_link_func}, {...}, ...]
End of explanation
m_reg = hddm.HDDMRegressor(mydata, reg_descr, include='z')
Explanation: The last step before running the model is to construct the complete hddm regression model by adding data etc.
End of explanation
m_reg.sample(5000, burn=200, dbname='within_effect.db', db='pickle') 
m_reg.save('within_effect')
Explanation: Now, we start the sampling, which takes a long time
End of explanation
m_reg.print_stats()
Explanation: Comparing generative and recovered model parameters
First, let's print the model stats
End of explanation
# load data
df1a = hddm.load_csv('df1a.v.hddm_stim.csv')
df = df1a
# randomly chose 10 participants from the dataset
df_subj = df['subj_idx'].unique()
random.seed(10)
df_test_list = []
for i in range(10):
    pos = random.randint(0, (len(df_subj)-1))
    df_test_list.append(df_subj[pos])  
df_test = df[df['subj_idx'].isin(df_test_list)]
def z_link_func(x, data=mydata):
    stim = (dmatrix('0 + C(s, [[1], [-1]])',
                    {'s':data.stimulus.loc[x.index]},
                    return_type="dataframe")
           )
    return 1/(1+np.exp(-np.multiply(x.to_frame(), stim)))
z_reg = {'model': 'z ~ 1 + C(condition)', 'link_func': z_link_func}
v_reg = {'model': 'v ~ 1 + C(condition)', 'link_func': lambda x:x}
a_reg = {'model': 'a ~ 1 + C(condition)', 'link_func': lambda x:x}
reg_descr = [z_reg, v_reg]
m_reg = hddm.HDDMRegressor(mydata, reg_descr, include='z')
m_reg.sample(5000, burn=200, dbname='within_effect.db', db='pickle') 
m_within_subj = hddm.HDDMRegressor(df1a, "v ~ C(match, val, Treatment(’Mismatch.Neutral’))")
Explanation: Working with real data
Now let's move to the real data
End of explanation
<END_TASK>