mSignal = 1; %mean of signal distribution
mNoise = 0; %mean of noise distribution
sigma = 1; %standard deviation of both distributions

alpha = log(.9/.1); %decision threshold; corresponds to 90% accuracy

E0 = log(2/1); %starting point; corresponds to prior belief of 2/3 signal, 1/3 noise

tau = .1; %time constant; 100 ms

N = 10000; %number of trials to simulate

H = 1 + (rand(N,1) > 2/3); %1 for signal, 2 for noise; setting 2/3 baserate to match prior (though they don't have to match)

r = zeros(N,1); %response for each trial
RT = zeros(N,1); %RT for each trial

for i = 1:N %loop through trials
    E = E0; %evidence begins at starting point
    n = 0; %tracks number of observations taken on this trial
    while abs(E) < alpha
        if H(i) == 1 %signal trial
            x = mSignal + randn*sigma; %observation
        else %noise trial
            x = mNoise + randn*sigma; %observation
        end

        L = (mSignal-mNoise)/sigma^2 * (x-(mSignal+mNoise)/2); %loglikelihood ratio

        E = E + L; %increment evidence by loglikelihood of current observation
        n = n + 1; %increment count of observations
    end

    r(i) = 1 + (E<0); %response: 1 if E positive (hit upper threshold), 2 if E negative (hit lower threshold)
    RT(i) = n*tau; %RT for this trial
end

rProbs = zeros(2,2); %response probabilities; trial type (signal, noise) by response (signal, noise)
for i=1:2,for j=1:2 %correct answer i, actual answer j
    rProbs(i,j) = mean(r(H==i)==j); %proportion of response on trials of type i that are response j
end,end
disp(' ')
disp('Response probabilities on signal trials')
disp(['Signal response: ' num2str(rProbs(1,1)) '; Noise response: ' num2str(rProbs(1,2))])
disp('Response probabilities on noise trials')
disp(['Signal response: ' num2str(rProbs(2,1)) '; Noise response: ' num2str(rProbs(2,2))])

figure(1),clf
for i=1:2,for j=1:2 %correct answer i, actual answer j
    subplot(2,2,(i-1)*2+j)
    hist(RT(H==i & r==j),tau/2:tau:max(RT)+tau) %RTs for this trial type and this response, with bin size equal to tau
    set(gca,'XLim',[0 max(RT)],'Ylim',[0 max(histc(RT,0:tau:max(RT)))]) %put all graphs on same scale
    if i==1&&j==1,ylabel('Signal Trials'),end
    if i==2&&j==1,xlabel('Signal Responses'),ylabel('Noise Trials'),end
    if i==2&&j==2,xlabel('Noise Responses'),end
end,end

