Dateien hochladen nach ''

data_fit_paper.m

@@ -0,0 +1,214 @@
+% fit model to data
+
+load('rep_randwalk.mat')
+nsub=max(allArray.NSub);
+
+secondhalf=false;
+if secondhalf
+    idx=allArray.NT>600 | (allArray.NT>200 & allArray.NT<400);
+    %idx=~idx; % take first instead of second half
+    idx=allArray.RandLevel==2 | idx; % take half, but only for random walk
+    allArray=allArray(idx,:);
+end
+
+% round physical presented duration (100 ms)
+allArray.duration = round(allArray.dur *10)/10;
+durs=unique(allArray.duration(allArray.NSub==1));
+allArray.sim=zeros(size(allArray.duration));
+
+% take only valid data
+% allArray = allArray(allArray.valid==1,:);
+% alternative, leaves temporal structure
+allArray.repDur(allArray.valid==0)=NaN;
+
+cond={'random walk','randomized'};
+parfit=zeros(nsub,2,1);
+linfit=zeros(nsub,2,2);
+simfit=zeros(nsub,2,2);
+
+
+tic
+for j=2:-1:1
+    figure
+    for i=1:nsub
+        idx=(allArray.NSub==i) & (allArray.RandLevel==j);
+        
+        % original model, just 1 parameter
+        [px,ci,resnorm,simres]=fitdata1pv([allArray.dur(idx) allArray.repDur(idx)]);
+        parfit(i,j,1)=px;
+        err=simres-allArray.repDur(idx);
+        idn=isfinite(err);
+        
+        stid=allArray.dur(idx);
+        repd=allArray.repDur(idx);
+        linfit(i,j,:)=polyfit(stid(idn),repd(idn),1);
+        
+        if j==1
+            % predict response from parameter fitted for other condition
+            [~, ~, ~, simres, ~]=kmodel1pv(parfit(i,3-j,1),[allArray.dur(idx) allArray.repDur(idx)]);
+        else
+            % simulate response with fitted parameters for this condition
+            [~, ~, ~, simres, ~]=kmodel1pv(parfit(i,j,1),[allArray.dur(idx) allArray.repDur(idx)]);
+        end
+        
+        allArray.sim(idx)=simres;
+        
+        simfit(i,j,:)=polyfit(stid(idn),simres(idn),1);
+        
+        subplot(4,4,i)
+        prcterr=(allArray.repDur(idx)-allArray.dur(idx))./allArray.dur(idx)*100;
+        prcterrsim=(simres-allArray.dur(idx))./allArray.dur(idx)*100;
+        
+        plot(allArray.dur(idx),allArray.repDur(idx),'.',allArray.dur(idx),simres,'.')
+        ylim([0 3])
+        
+        %plot(allArray.dur(idx),prcterr,'.',allArray.dur(idx),prcterrsim,'.')
+        %ylim([-50 50])
+        
+        meanprct=grpstats(prcterr,allArray.duration(idx));
+        [meanrepdur,stdrepdur]=grpstats(allArray.repDur(idx),allArray.duration(idx),{'mean','std'});
+        durs=unique(allArray.duration(idx));
+        meanprctsim=grpstats(prcterrsim,allArray.duration(idx));
+        [meanrepsim,stdrepsim]=grpstats(simres,allArray.duration(idx),{'mean','std'});
+    end
+end
+toc
+
+
+%%
+i=nsub;
+% i=5;
+% i=8;
+figure('name',['data from subject ' int2str(i)])
+clr=colororder;
+% randomized
+j=2;
+idx=(allArray.NSub==i) & (allArray.RandLevel==j);
+subplot(1,2,1)
+hold on
+plot([0 2],linfit(i,j,1)*[0 2]+linfit(i,j,2),'Color',clr(1,:),'linewidth',2)
+plot([0 2],simfit(i,j,1)*[0 2]+simfit(i,j,2),'--','Color',clr(2,:),'linewidth',2)
+plot(allArray.dur(idx),allArray.repDur(idx),'.','Color',clr(1,:))
+plot(allArray.dur(idx),allArray.sim(idx),'.','Color',clr(2,:))
+plot([0 3],[0 3],'--k')
+hold off
+xlim([0 2])
+ylim([0 2.5])
+xlabel('stimulus duration (s)')
+ylabel('reproduced duration (s)')
+legend('experiment','model simulation')
+title('randomized condition')
+set(gca,'Fontsize',16)
+text('Parent',gca,'FontSize',36,'String','A','Position',[-0.4 2.5 0]);
+
+% calculate error (response-actual)
+errA=allArray.repDur(idx)-allArray.dur(idx);
+presdurA=allArray.dur(idx);
+simerrA=allArray.sim(idx)-allArray.dur(idx);
+
+% now for random walk
+j=1;
+idx=(allArray.NSub==i) & (allArray.RandLevel==j);
+subplot(1,2,2)
+hold on
+plot([0 2],linfit(i,j,1)*[0 2]+linfit(i,j,2),'Color',clr(1,:),'linewidth',2)
+plot([0 2],simfit(i,j,1)*[0 2]+simfit(i,j,2),'--','Color',clr(2,:),'linewidth',2)
+plot(allArray.dur(idx),allArray.repDur(idx),'.','Color',clr(1,:))
+plot(allArray.dur(idx),allArray.sim(idx),'.','Color',clr(2,:))
+plot([0 3],[0 3],'--k')
+hold off
+xlim([0 2])
+ylim([0 2.5])
+xlabel('stimulus duration (s)')
+ylabel('reproduced duration (s)')
+legend('experiment','model prediction')
+title('random walk condition')
+set(gca,'Fontsize',16)
+text('Parent',gca,'FontSize',36,'String','B','Position',[-0.4 2.5 0]);
+set(gcf,'Position',[560   556   836   392])
+
+% calculate error (response-actual)
+errB=allArray.repDur(idx)-allArray.dur(idx);
+presdurB=allArray.dur(idx);
+simerrB=allArray.sim(idx)-allArray.dur(idx);
+
+%%
+figure('name',['data from subject ' int2str(i) ' over trial'])
+subplot(2,1,1)
+j=2;
+idx=(allArray.NSub==i) & (allArray.RandLevel==j);
+plot([allArray.dur(idx),allArray.repDur(idx),allArray.sim(idx)],'.-','Markersize',10)
+xlabel('trial')
+ylabel('duration (s)')
+legend('stimulus','reproduction','simulation')
+set(gca,'Fontsize',16)
+title('randomized')
+subplot(2,1,2)
+j=1;
+idx=(allArray.NSub==i) & (allArray.RandLevel==j);
+plot([allArray.dur(idx),allArray.repDur(idx),allArray.sim(idx)],'.-','Markersize',10)
+xlabel('trial')
+ylabel('duration (s)')
+legend('stimulus','reproduction','prediction')
+set(gca,'Fontsize',16)
+title('random walk')
+
+%%
+allArray.repError = allArray.repDur - allArray.dur;
+allArray.simError = allArray.sim - allArray.dur;
+
+avgdata = grpstats(allArray, {'NSub','RandLevel','duration'},{'mean','sem'}, 'DataVars','repError');
+avgdatadur = grpstats(allArray, {'NSub','RandLevel','duration'},{'mean','sem'}, 'DataVars','repDur');
+
+avgsim = grpstats(allArray, {'NSub','RandLevel','duration'},{'mean','sem'}, 'DataVars','simError');
+avgsimdur = grpstats(allArray, {'NSub','RandLevel','duration'},{'mean','sem'}, 'DataVars','sim');
+
+avgdata.prcterr=avgdata.mean_repError./avgdata.duration*100;
+avgsim.prcterr=avgsim.mean_simError./avgsim.duration*100;
+
+gavgdata = grpstats(avgdata, {'RandLevel','duration'},{'mean','sem'}, 'DataVars','mean_repError');
+gavgdatadur = grpstats(avgdatadur, {'RandLevel','duration'},{'mean','sem'}, 'DataVars','mean_repDur');
+gavgsim = grpstats(avgsim, {'RandLevel','duration'},{'mean','sem'}, 'DataVars','mean_simError');
+gavgsimdur = grpstats(avgsimdur, {'RandLevel','duration'},{'mean','sem'}, 'DataVars','mean_sim');
+
+gavgdata.prcterr=gavgdata.mean_mean_repError./gavgdata.duration*100;
+gavgsim.prcterr=gavgsim.mean_mean_simError./gavgsim.duration*100;
+
+gavgprct=grpstats(avgdata, {'RandLevel','duration'},{'mean','sem'}, 'DataVars','prcterr');
+gavgsimprct=grpstats(avgsim, {'RandLevel','duration'},{'mean','sem'}, 'DataVars','prcterr');
+
+
+%% plot for Vierordt paper
+figure
+
+idx=gavgdata.GroupCount>2;
+subplot(1,2,1)
+hold on
+errorbar(gavgprct.duration(gavgprct.RandLevel == 1 & idx)-0.01, gavgprct.mean_prcterr(gavgprct.RandLevel == 1 & idx),gavgprct.sem_prcterr(gavgprct.RandLevel == 1 & idx),'k-o','MarkerFaceColor','k','MarkerSize',10);
+errorbar(gavgprct.duration(gavgprct.RandLevel == 2 & idx)+0.01, gavgprct.mean_prcterr(gavgprct.RandLevel == 2 & idx),gavgprct.sem_prcterr(gavgprct.RandLevel == 2 & idx),'k-o','MarkerSize',10);
+plot([0 10],[0 0],'--k')
+xlim([0 2])
+ylim([-30 60])
+title('experiment')
+legend('random walk','randomized')
+xlabel('mean target time (s)')
+ylabel('percentage error of reproduction')
+set(gca,'Fontsize',16)
+text('Parent',gca,'FontSize',36,'String','A','Position',[-0.4 63 0]);
+hold off
+
+subplot(1,2,2)
+hold on
+errorbar(gavgsimprct.duration(gavgsimprct.RandLevel == 1 & idx)-0.01, gavgsimprct.mean_prcterr(gavgsimprct.RandLevel == 1 & idx),gavgsimprct.sem_prcterr(gavgsimprct.RandLevel == 1 & idx),'k-o','MarkerFaceColor','k','MarkerSize',10);
+errorbar(gavgsimprct.duration(gavgsimprct.RandLevel == 2 & idx)+0.01, gavgsimprct.mean_prcterr(gavgsimprct.RandLevel == 2 & idx),gavgsimprct.sem_prcterr(gavgsimprct.RandLevel == 2 & idx),'k-o','MarkerSize',10);
+plot([0 10],[0 0],'--k')
+xlim([0 2])
+ylim([-30 60])
+title('simulation')
+xlabel('mean target time (s)')
+ylabel('percentage error of reproduction')
+set(gca,'Fontsize',16)
+text('Parent',gca,'FontSize',36,'String','B','Position',[-0.4 63 0]);
+hold off
+set(gcf,'Position',[560   556   836   392])
+

fitdata1pv.m

@@ -0,0 +1,15 @@
+function [px,ci,resnorm,simres]=fitdata1pv(stimrep)
+% input stimrep : stimulus and response
+% output px : parameters of the model kmodel(p,stimrep)
+% output ci : confidence interval of px
+%
+% fit wrapper for the two-stage model published in
+% Petzschner FH, Glasauer S. Iterative bayesian estimation as an explanation for range and regression effects: a study on human path integration. J Neurosci. 2011 Nov 23;31(47):17220-9.
+% 
+% S.Glasauer 2011/2021
+p0=1;lb=0;
+[px,resnorm,residual,~,~,~,jacobian]=lsqnonlin(@(p)kmodel1pv(p,stimrep),p0,lb);
+ci=nlparci(px,residual,'jacobian',jacobian);
+[~,~,~,simres]=kmodel1pv(px,stimrep);
+% simres=stimrep(:,2)-kmodel(px,stimrep);
+end

kmodel1pv.m

@@ -0,0 +1,63 @@
+function [sres, xest, pest, resp, wp]=kmodel1pv(par,stimrep)
+% input par : current parameter of the model, here only q/r
+% input stimrep : stimulus stimrep(:,1) and response stimrep(:,2)
+% output sres : residual of simulation and response
+% simulation=resp-sres
+
+% par(1) : ratio q/r for kalman filter
+
+% two-stage model published in
+% Petzschner FH, Glasauer S. Iterative bayesian estimation as an explanation 
+% for range and regression effects: a study on human path integration. 
+% J Neurosci. 2011 Nov 23;31(47):17220-9.
+% 
+% S.Glasauer 2011/2021
+
+r=1;
+q=par(1)*r;
+p=r;
+
+a=10;
+off=1;
+
+dist=stimrep(:,1);
+pest=dist*0;
+xest=dist*0;
+resp=dist;
+wp=dist*0;
+for i=1:length(dist) % walk each distance in randomized order
+    d=dist(i);
+    
+    % this is the core estimation routine
+    % measurement
+    z=log(a*d+off); % transform distance to log
+    if i==1 % initial mean prior equals measurement d
+        x=z;
+    end
+    % kalman filter for the mean of the prior
+    km=(p+q)/(p+q+r); % kalman gain
+    p=km*r; % new variance
+    pest(i)=p;
+    %x=(1-km)*x+km*z; % prior
+    x=x+km*(z-x);
+    xest(i)=(exp(x)-off)/a;
+
+    % final fusion
+    w1=1/r/(1/p+1/r); % weight of measurement
+    xe=(1-w1)*x+w1*z; % estimate
+    wp(i)=1-w1; % weight of prior
+    % backtransform
+    resp(i)=(exp(xe)-off)/a; % estimate backtransform    
+end
+
+sres=stimrep(:,2)-resp;
+% analytical steady state solutions
+% k_final=0.5*q/r*(sqrt(1+4*r/q)-1);
+% final prior weight: w_prior=(1-k_final)/(1+k_final);
+
+% remove NaNs ... this treatment of NaNs should lead to better estimation
+% than just setting the residuals to zero, since here the missing values 
+% are really treated as missing
+sres=sres(isfinite(sres));
+
+end