% generate data for random walk simulation
clear
 
% n = 300;
n = 30
ns = 20; ng=3;
 
d = [-.75 ; .75];  % noise signal
th =1; z = .2; stdd=.75; sth=.2;  
 
 
 
type = 1;
 
if type ==1
datam = zeros(n,4,ns,ng);
% n x 4 : 
% n trials
% 4 cond's: noise-up noise-down sig-up sig-down
% nstim-Sresp nstim-Nresp sstim-Sresp sstim-Nresp
for age = 1:ng
    age        
    mt = (age/ng);      
        for sub = 1:ns  
        thr = th + sth*(age-1); 
        choice = simulate(n,d,thr,z,stdd,mt);
        datam(:,:,sub,age) = choice;
        end
    thr
end
 
else
    load simdat;
end
 
 
run data_anal;  
run summary;
% run test_diff
save simdat2 datam Ym Mdat 
 
 
 

%%%%%%%%%%%%%%%%%%%%%%%%%%% functions %%%%%%%%%%%%%%%%%%%
function choice = simulate(n,d,th,z,stdd,mt)
% [choice, Rt] = simulate(n,d,th,z,stdd,mt,sth)
% produces a sample of n simulations from a diffusion process
 
h = .0001;
choice = zeros(n,2,2);
for c = 1:2  % pick a d value 
for k = 1:n  % run trials
    
   t = mt;  
   x = z;
   dd = d(c) + stdd*randn(1);
   while abs(x) < th
      t=t+h;
      x = x + dd.*h + sqrt(h).*randn(1);  
   end  % a single trial
   choice(k,:,c) = t.*[(x>0) (x<0)];  
   % [up = resp sig; down = resp noise]
   
end  % all trials
end  % both d's
choice = [choice(:,:,1) choice(:,:,2)];
% trial, cond, d
% cond = noise-up noise-down sig-up sig-down 


% data analysis
ng=3;
ns=20;
 
p = [0; .25; .5; .75; 1];
sz= size(p);
dp = p(2:sz)-p(1:sz-1);
p = p(2:sz-1);
Ym=zeros(ng,ns,4,7);
nt = size(datam(:,:,1,1),1);
for j=1:ng  % 3
    for i = 1:ns   % 20
        X = datam(:,:,i,j);
        fr = sum(X>0);
        Nj = dp*fr ;
        Ij = p*fr + .5; Ijp = ceil(Ij); Ijn = floor(Ij);
        Ijn = (Ijn==0)*1 + (Ijn>0).*Ijn;
        X1 = sort(X(find(X(:,1)),1));
        X2 = sort(X(find(X(:,2)),2));
        X3 = sort(X(find(X(:,3)),3));
        X4 = sort(X(find(X(:,4)),4));
        if size(X4) == 1
            X4 = [X4;X4;X4];
        elseif size(X4) == 2
            X4 = [X4; X4(2)];
        else
        end
        C1 = X1(Ijn(:,1))+(X1(Ijp(:,1))-X1(Ijn(:,1))).*(Ij(:,1)-Ijn(:,1));  %noise up
        C2 = X2(Ijn(:,2))+(X2(Ijp(:,2))-X2(Ijn(:,2))).*(Ij(:,2)-Ijn(:,2));  %noise down
        C3 = X3(Ijn(:,3))+(X3(Ijp(:,3))-X3(Ijn(:,3))).*(Ij(:,3)-Ijn(:,3));  % sig up
        C4 = X4(Ijn(:,4))+(X4(Ijp(:,4))-X4(Ijn(:,4))).*(Ij(:,4)-Ijn(:,4));  % sig down
        
        quart1 = ones(nt,1)*[C1(1) C2(1) C3(1) C4(1)];
        quart2 = ones(nt,1)*[C1(2) C2(2) C3(2) C4(2)];
        quart3 = ones(nt,1)*[C1(3) C2(3) C3(3) C4(3)];
        
        Nj(1,:) = sum((X>0).*(X<=quart1));
        Nj(2,:) = sum((X>quart1).*(X<=quart2));
        Nj(3,:) = sum((X>quart2).*(X<=quart3));
        Nj(4,:) = sum(X>quart3);
        
        Y = [ [Nj(:,1);C1] [Nj(:,2);C2]  [Nj(:,3);C3] [Nj(:,4);C4] ];   
        Ym(j,i,:,:) = Y';
        
    end 
end
 
% Ym( age, subj, 4 , 7)
%  4 rows: (C1=noise up FA, C2=noise down CR, C3=sig up HIT C4=sig down Miss)
%  7 cols: ( freq for 4 categories, 3 quartiles)
 
 
%summary statistics
 % Ym(age,subj, 4 conditions, 7 measures)
 % 4 conditions: (noise up, noise down, sig up, sig down)
 % 7 measures: (4 n's per quartile, 3 quartiles)
 
 Z=squeeze(sum(Ym(:,:,:,1:4),4));  % sum over four n's
C1 = Z(:,:,1)./(Z(:,:,1)+Z(:,:,2));  % prob up for noise 
C2 = Z(:,:,3)./(Z(:,:,3)+Z(:,:,4));  % prob up for sig
mC1 = mean(C1,2);  % mean over subj for each group , noise
mC2 = mean(C2,2);  % mean over subj for each group , sig
' for each group, prob up given noise, prob up given sig'
Pc = [mC1 mC2]
 
' mrt per group for ( up noise, down noise, up sig, down sig)'
mrt = squeeze(mean(Ym(:,:,:,6),2))  % mean over subj of medians 
 
Nq = (n/4)*[Pc  1-Pc];
Nq = cat(3,Nq,Nq,Nq,Nq);
Mq = squeeze(mean(Ym(:,:,:,5:7),2));
Mdat = cat(3,Nq,Mq);
 
 
% (mrt(:,1:2)+mrt(:,3:4))/2  % group by resp
% m1 = mrt(:,1).*Pc(:,1) + mrt(:,2).*(1-Pc(:,1));
% m2 = mrt(:,3).*Pc(:,2) + mrt(:,4).*(1-Pc(:,2));
% (m1+m2)/2