-
Notifications
You must be signed in to change notification settings - Fork 2
/
dmcSim.m
executable file
·200 lines (180 loc) · 7.17 KB
/
dmcSim.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
function res = dmcSim(varargin)
% res = dmcSim(varargin)
%
% DMC model simulation
% Ulrich, R., Schröter, H., Leuthold, H., & Birngruber, T. (2015).
% Automatic and controlled stimulus processing in conflict
% tasks: Superimposed diffusion processes and delta functions.
% Cognitive psychology, 78, 148-174.
% Code adapted from Appendix C. Basic Matlab Code.
%
% Tested using Matlab 2017a+ (implicit broadcasting required from Matlab 2016b+)
%
% varargin:
% 'amp', amplitude of automatic activation
% 'tau', time to peak automatic activation
% 'aaShape', shape parameter of automatic activation
% 'mu', drift rate of controlled processes
% 'sigma', diffusion constant
% 'bnds', +- response barrier
% 'resMean', mean of non-decisional component
% 'resSD', standard deviation of non-decisional component
% 'tmax', integer
% 'varSP', true/false variable starting point
% 'spShape', shape parameter of starting point distribution
% 'varDR', true/false variable drift rate
% 'drLim', limit range of distribution of drift rate
% 'drShape', shape parameter of drift rate
% 'nTrl', integer between ~1000 and 200000
% 'nTrlData', integer beteen ~5 and 10
% 'cafBins', bins to calculate conditional accuracy functions
% 'deltaBins', bins to calculate incomp-comp delta plots
% 'makePlots', true/false
%
% Outputs:
% struct with the following fields:
% activation: drift, comp, incomp
% trials: comp, incomp
% rts: comp, incomp with n rows (trials) and 3 columns (rt, error, bin)
% caf: 2 rows (comp, incomp * n columns (bins)
% rtDist: 4 rows (comp, incomp, mean, effect) * n columns (bins)
% prms: model input parameters
% summary: table with comp, incomp with mean/sd/percentage errors
%
% Reproduces Figures 3 (30), 4 (150), & 5 (90) with changes in tau parameter (see Table 1)
% Reproduces Figure 6 with varSP (variable start point)
% Reproduces Figure 7 with varDR (variable drift rate)
% Reproduces Figure 8 (CAF) varies with changes in bounds (speed/accuracy)
%
% Examples:
% res = dmcSim(); Fig 3
% res = dmcSim('tau', 150); Fig 4
% res = dmcSim('tau', 90); Fig 5
% res = dmcSim('varSP', true); Fig 6
% res = dmcSim('varDR', true); Fig 7
% res = dmcSim('amp', 25, 'tau', 15, 'resMean', 350, ...)
%% setup
prms.amp = 20; % amplitude of automatic activation
prms.tau = 30; % time to peak automatic activation
prms.aaShape = 2; % shape parameter of automatic activation
prms.mu = 0.5; % drift rate of controlled processes
prms.sigma = 4; % diffusion constant
prms.bnds = 75; % +- response barrier
prms.resMean = 300; % mean of non-decisional component
prms.resSD = 30; % standard deviation of non-decisional component
prms.tmax = 1000; % vector length
prms.varSP = false; % variable start point
prms.spShape = 3; % shape parameter of variable start point
prms.varDR = false; % variable drift rate
prms.drLim = [0.1 0.7]; % range of beta distribution
prms.drShape = 3; % shape parameter of variable drift rate
prms.nTrl = 100000; % number of trials within each comp/incomp conditions
prms.nTrlData = 5; % number of individual trials to return
prms.cafBins = 0:20:100; % bins to calculate conditional accuracy functions
prms.deltaBins = 5:10:95; % bins to calculate incomp-comp delta plots
prms.makePlots = true;
for i = 1:2:length(varargin)
switch varargin{i}
case 'amp'
prms.amp = varargin{i+1};
case 'tau'
prms.tau = varargin{i+1};
case 'aaShape'
prms.aaShape = varargin{i+1};
case 'mu'
prms.mu = varargin{i+1};
case 'sigma'
prms.sigma = varargin{i+1};
case 'bnds'
prms.bnds = varargin{i+1};
case 'resMean'
prms.resMean = varargin{i+1};
case 'resSD'
prms.resSD = varargin{i+1};
case 'tmax'
prms.tmax = varargin{i+1};
case 'varSP'
prms.varSP = varargin{i+1};
case 'spShape'
prms.spShape = varargin{i+1};
case 'varDR'
prms.varDR = varargin{i+1};
case 'drLim'
prms.drLim = varargin{i+1};
case 'drShape'
prms.drShape = varargin{i+1};
case'nTrl'
prms.nTrl = varargin{i+1};
case'nTrlData'
prms.nTrlData = varargin{i+1};
case 'cafBins'
prms.cafBins = varargin{i+1};
case 'deltaBins'
prms.deltaBins = varargin{i+1};
case'makePlots'
prms.makePlots = varargin{i+1};
otherwise
error('varargin not recognised');
end
end
%% create results structure
res = resStruct(prms);
%% simulation
drift = prms.amp .* exp(-(1:prms.tmax) ./ prms.tau) .* ((exp(1) .* (1:prms.tmax) ./ (prms.aaShape-1) ./prms.tau) .^ (prms.aaShape-1));
for comp = {'comp', 'incomp'}
if strcmp(comp, 'comp'); sign = 1; else, sign = -1; end
if ~prms.varDR % constant drift rate across trials
muVec = (sign * drift) .* ((prms.aaShape-1) ./ (1:prms.tmax)-1/prms.tau) + prms.mu; % eq7/eq8
else % variable drift rate (vdr): beta distribution
muVec = (sign * drift) .* ((prms.aaShape-1) ./ (1:prms.tmax)-1/prms.tau) + rand_beta(prms.nTrl, prms.drShape, prms.drLim); % eq7/eq8
end
activation = muVec + (prms.sigma*randn(prms.nTrl, prms.tmax, 'single'));
if prms.varSP % variable starting point: beta distribution (+- bounds)
activation(:, 1) = activation(:, 1) + rand_beta(prms.nTrl, prms.spShape, [-prms.bnds prms.bnds]);
end
% accumulate activation
activation = cumsum(activation, 2);
%% find RTs for correct/incorrect trials (point where X exceeds crit +-bounds)
[~, rt] = max(abs(activation) > prms.bnds, [], 2);
rt(rt == 1, 1) = prms.tmax; % does not reach boundary classified as error
rt_idx = sub2ind(size(activation), 1:size(activation, 1), rt');
rt = [rt + normrnd(prms.resMean, prms.resSD, prms.nTrl, 1), (activation(rt_idx) < prms.bnds)'];
%% calculate conditional accuracy functions (CAF)
[~, ~, rt(:, 3)] = histcounts(rt(:, 1), prctile(rt(:, 1), prms.cafBins));
res.caf = [res.caf; transpose(calcCAF(rt))];
%% calculate percentiles
res.rtDist = [res.rtDist; prctile(rt(rt(:, 2) == 0, 1), prms.deltaBins)];
%% store required results
res.activation.drift = drift;
res.activation.(comp{:}) = mean(activation);
res.trials.(comp{:}) = activation(1:prms.nTrlData, :);
res.rts.(comp{:}) = rt;
res.summary = [res.summary; [comp(:)', ...
{round(mean(rt(rt(:, 2) == 0)))}, ...
{round(std(rt(rt(:, 2) == 0)))}, ...
{round(sum(rt(:,2)/prms.nTrl) * 100, 1)}, ...
{round(mean(rt(rt(:, 2) == 1)))}, ...
{round(std(rt(rt(:, 2) == 1)))}]];
end
%% calculate comp effect and delta
res.rtDist(3,:) = (res.rtDist(1, :) + res.rtDist(2, :)) / 2;
res.rtDist(4,:) = (res.rtDist(2, :) - res.rtDist(1, :));
%% plots
if prms.makePlots
plotDMC(res);
end
end
%%
function res = resStruct(prms)
% res = resStruct(prms)
%
% Creates output structure for dmcSim.
res.prms = prms;
res.activation = [];
res.trials = [];
res.rts = [];
res.caf = [];
res.rtDist = [];
res.summary = array2table(nan(0, 6), ...
'VariableNames', {'Comp', 'rtCorr', 'sdRtCorr', 'perErr', 'rtErr', 'sdRtErr'});
end