-
Notifications
You must be signed in to change notification settings - Fork 4
/
driftdiff.m
95 lines (76 loc) · 3.02 KB
/
driftdiff.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
function [RT,decision,evidence,pm] = driftdiff(pm)
%function [RT,decision,evidence] = driftdiff(pm.['trials','driftrate','noise','zeropoint','bias','biasrange','upperbound','nondectime','ndcrange','deadline','leak','lambda'])
%
% Simulate drift diffusion process (+ optional plot).
%
% Input:
% - pm
% - trials: number of trials [1000]
% - driftrate: stimulus quality [0]
% - noise: SD of driftrate [10]
% - zeropoint: drift criterion [0]
% - bias: prior bias[.5]
% - biasrange: range of bias [0]
% - upperbound: upper bound a [500]
% - nondectime: non-decision time [0]
% - ndcrange: range of non-decision time [0]
% - deadline: length of epoch in ms [1000]
% - leak: decay of information, reasonably between 0 - 1 [0]
% - lambda: Ornstein-Uhlbeck decay in drift (mean DV*l) [0]
%
% Output:
% - RT: reaction times for each trial
% - decision: choice (-1,0,1)
% - evidence: evidence trace for each trial
%
% Fabrice Luyckx, 23/3/2017 (adapted from scripts by Chris Summerfield)
%% DEFAULT VALUES
paramnames = {'trials','driftrate','noise','zeropoint','bias','biasrange','upperbound','nondectime','ndcrange','deadline','leak','lambda'};
paramvals = {1000,0,10,0,.5,0,500,0,0,1000,0,0};
structnames = fields(pm);
for i = 1:length(paramnames)
if sum(strcmp(paramnames{i},structnames)) == 0 % field doesn't exist
pm = setfield(pm,paramnames{i},paramvals{i});
else
if isempty(getfield(pm,paramnames{i})) % field has no content
pm = rmfield(pm,paramnames{i});
pm = setfield(pm,paramnames{i},paramvals{i});
end
end
end
%% Initialise
RT = nan(pm.trials,1);
decision = zeros(pm.trials,1);
evidence = zeros(pm.trials,pm.deadline);
% Set drift rate
if length(pm.driftrate) == 1
pm.driftrate = pm.driftrate.*ones(pm.trials,1);
end
% Set zeropoint
if length(pm.zeropoint) == 1
pm.zeropoint = pm.zeropoint.*ones(pm.trials,1);
end
% Draw nondecision time from uniform distribution
pm.nondectime = pm.nondectime.*ones(pm.trials,1) + (pm.ndcrange*rand(pm.trials,1));
% Determine bias from uniform distribution
pm.bias = pm.upperbound.*(pm.bias.*ones(pm.trials,1) + (pm.biasrange*rand(pm.trials,1)));
% Evidence starts as bias
DV = pm.bias;
%% Simulate drift process
for t = 1:pm.deadline
% Add increment
inc = pm.driftrate + pm.zeropoint + randn(pm.trials,1)*pm.noise; % increment
inc = inc.*(t>pm.nondectime); % only add increment after nondecision time has passed
DV = DV + inc - (pm.leak*sign(DV)) + (mean(DV)*pm.lambda); % add leak and increase rate
% Log trials where bound is reached (a or 0)
RT(DV > pm.upperbound & decision == 0) = t;
RT(DV < 0 & decision == 0) = t;
% Log decision
decision(DV > pm.upperbound & decision == 0) = 1;
decision(DV < 0 & decision == 0) = -1;
% Log evidence
evidence(:,t) = DV;
evidence(decision == 1,t) = pm.upperbound;
evidence(decision == -1,t) = 0;
end
end