rng(0); %Include seed so it will reproduce exactly
NumTrials = 240;
grain = 201;
PF = @PAL_Gumbel; %assumed psychometric function

stimRange = (linspace(PF([0 1 0 0],.1,'inverse'),PF([0 1 0 0],.9999,'inverse'),21));
priorAlphaRange = linspace(PF([0 1 0 0],.1,'inverse'),PF([0 1 0 0],.9999,'inverse'),grain);
priorBetaRange =  linspace(log10(.0625),log10(16),grain);
priorGammaRange = 0.5;
priorLambdaRange = 0:.01:.2;

%Initialize PM structure
PM = PAL_AMPM_setupPM('priorAlphaRange',priorAlphaRange,...
                      'priorBetaRange',priorBetaRange,...
                      'priorGammaRange',priorGammaRange,...
                      'priorLambdaRange',priorLambdaRange,...
                      'numtrials',NumTrials,...
                      'PF' , PF,...
                      'stimRange',stimRange);                  

paramsGen = [0, 1, .5, .02];   %generating parameter values [alpha, beta, gamma, lambda] (or [threshold, slope, guess, lapse])
                                
while PM.stop ~= 1

    response = rand(1) < PF(paramsGen, PM.xCurrent);    %simulate observer
    %update PM based on response:
    PM = PAL_AMPM_updatePM(PM,response);

end

data.x = PM.x(1:end-1);
data.y = PM. response;
data.n = ones(size(data.x));
pfhb = PAL_PFHB_fitModel(data,'PF','Gumbel');

plot(priorAlphaRange,PF([PM.threshold(end) 10.^PM.slope(end) PM.guess(end) PM.lapse(end)],priorAlphaRange))
hold on
plot(priorAlphaRange,PF([pfhb.summStats.a.mean 10.^pfhb.summStats.b.mean .5 pfhb.summStats.l.mean],priorAlphaRange))