% script m-file GottliebAnal.m computes Bayes factors for Model that says
% there is no difference vs a model that says that the mean off the Many
% group is on the order of 10% greater than the mean of the Few group.
% Assumes that the sigma of the Gaussian likelihood function is known and is
% equal to the sigma of the best fit to the distribution of the combined
% data.

FewData=[1.3300;1.67;4.33;7.33;7.67;9.67;15.33;16.33];

ManyData=[4.3333;4;21.3333;7.6667;19.3333;0;7.6667;1.3333;...
    6;7.3333;8.6667;15.3333;11.3333;15.6667;7;0.3333];

%%    
CmbdData=[FewData;ManyData];

stairs(sort(CmbdData),(1:length(CmbdData))/length(CmbdData),'k')
% plot of cumulative distribution of combined data
hold on

stairs(sort(FewData),(1:length(FewData))/length(FewData),'g')
% plot of the cumulative distribution of the data from the Few Group

stairs(sort(ManyData),(1:length(ManyData))/length(ManyData),'r')
% plot of the cumulative distribution of the data from the Many Group

xlabel('Mean Differential Responses Per Minute')
ylabel('Cumulative Fraction of Subjects')
title('Gottlieb, Experiment 1, Cumulative Distributions')

MxLp=mle(CmbdData) % maximum likelihood parameters for the normal distribution
% when fitted to the combined data
% mean = 8.3747  sigma = 5.9644

% Although the normal distribution extends below 0, which is anlytically
% impossible, it has the advantage that it does allow observations of
% 0, and there is one in the data. Moreover, considerations of the
% step-like form of the learning curve for individual subjects implies that
% if there were a difference between the groups, it would take the form of
% there being more 0s in the Few Group than in the Many Group (that is,
% more subjects who had not yet learned).

RperMin=[0:.5:30]'; % x values for plotting

plot(RperMin,normcdf(RperMin,MxLp(1),MxLp(2)),'k--')
% superposing best-fitting Gaussian

legend('Combined','Many','Few','Best-Fit')


%% Strategy
% The Bayes factor for Model 1 (no diff) vs Model 2 (vaguely specified positive
% diff) is determined entirely by the relative likelihoods of the two
% models given the data from the Many group, because the likelihoods of the
% two models given the data for the Few Group are the same. The likelihood
% of the Many data given Model 1 can be computed using as a prior the
% marginal posterior density of the mean of the Few Group. This embodies
% the hypothesis that the Many data are drawn from the same distribution as
% the Few data and that our prior uncertainty about that distribution is
% specified by the 8 draws we have already made. We get this prior (that is,
% the posterior of the Few Group) using a uniform prior distribution on the
% interval between 0 and 20. To get the prior for the
% Many group, we convolve a prior for the increments with the just
% specified prior derived from the the Few data.  By convolving the increment prior with the prior
% derived from the Few data, we create a prior distribution for Model 2
% that takes into account our uncertainty about what the true value of the
% mean of the distribution from which the Few data were drawn really is.
%
% Computing the posterior marginal probability distribution for the mean
% of the Few Group, using a uniform prior. Because we will normalize the
% likelihood function to make it integrage to 1, we do not have to worry
% about the probability density of this uniform prior

M=linspace(0,20,200); % 200 points evenly spaced between 0 and 20

for i=1:length(FewData) % stepping through the Few data
    
    LgLkFew(i,:)=log10(normpdf(FewData(i),M,6));
    % log10(normpdf(FewData(i),M,6)) gives the log of the probability
    % density of the likelihood function at the datum, FewData(i), for each
    % of the 100 possible positions of the mean given in the row vector M.
    % log10(NormFactors) applies the appropriate normalizing factor for
    % each of these likelihoods
    
end % of stepping through the Few data

LgLkFnctnFew=sum(LgLkFew); % summing across the data

LkFnctnFew=10.^LgLkFnctnFew; % antilogs

LkFnctnFew=LkFnctnFew/sum(LkFnctnFew); % normalizing to get the posterior
% probability density function for the mean of the Few Group. When we
% normalize by summing and then dividing through by the sum, we treat the
% probability density function that we have computed as a set of rectangles
% (bars in a histogram) whose area must sum to 1. This normalization
% assumes that each of the points along the x axis at which we have
% computed a probability density has a unit width. Recall that in
% probability density functions when you rescale the x axis, you also must
% rescale the y axis, because probability density is probability mass PER
% UNIT ON THE X AXIS! In the present case, the units with respect to which
% we have normalized the probability density function are 1/10th of our
% responses/minute units on the x axis of our graphs (because the
% M = linspace(0,20,200) command created 10 points per unit

figure % posterior probability density for mean of Few Group

plot(M,LkFnctnFew*10) % the *10 rescales the probability densities to make
% them appropriate to the units on the x-axis of our graph
xlabel('Possible Values of Few Mean')
ylabel('Posterior Probability Density')
title('Posterior Density Function for Mean of Few Group')

hold on

[Mx,I]=max(LkFnctnFew);

FewMode=M(I) % maximum likelihood value for the mean of the Few group

plot([FewMode FewMode],ylim,'--') % dashed line at mode

plot([mean(FewData) mean(FewData)],ylim,'-.')

legend('Post PD','Mode','Sample Mean')

% The posterior probability density function for the mean of the Few Group
% is Model 1's prior probability distribution for the mean of the Many
% Group

M1prior=LkFnctnFew;

%%

figure % Cum distribution of Few data w maximum likelihood Gaussian
% likelihood function

stairs(sort(FewData),(1:length(FewData))/length(FewData),'g')

hold on

plot(M,normcdf(M,FewMode,6),'g--')

xlabel('Responses per Minute')
ylabel('Cumulative Fraction of Subjects')
title('Cumulative Distribution of Few Data w Max Like Gaussian')
legend('Data','Max Lkhd Gaussian')

%% the increment prior

% Instead of what I did in class, where I used a truncated normal function
% for the increment prior (which had substantial probability density at 0),
% I now use a Weibull distribution, which has 0 probability density at 0
% (hence, makes increments of 0 impossible). It also has a long tail, so it
% allows much bigger increments.

figure

IncrPrior=wblpdf(M,2,2);

IncrPrior=IncrPrior/sum(IncrPrior); % normalizing

plot(M,IncrPrior*10) % This distribution expresses in probabilistic terms
% our prior beliefs (or at least conceivable assumptions we could make)
% about the size of the effect of the 8-fold increase in number of trials
% on the mean strength of the CR
xlabel('Mean Increment in CR Strength')
ylabel('Probability Density')
title('Prior Distribution of the Increments')
%%
% To get the prior distribution for the mean of the Many Group, we have to
% integrate our prior belief (or possible assumption) about the size of the
% effect of the ADDITIONAL trials (the increment distribution) with our
% uncertainty about the size of the effect of the initial trials (the
% trials that the Few and Many Groups both got). Our uncertainty about the
% size of the effect of those initial trials is specified by the
% distribution for the mean of the Few Group (NB, not the distribution from
% which the Few sample was drawn, but rather the distribution of the means
% of samples of that size). We put them together by convolving them

M2prior=conv(IncrPrior,M1prior); 

M2prior=M2prior(1:200); % The number of points in the function
% obtained by convolving two vectors is the sum of their lengths minus 1.
% Only the first 100 entries fall on the same time stretch as the two
% functions that were convolved. This trims the M2prior
%
figure % M1 prior, Increment prior and M2 prior
plot(M,M1prior*10,'k'); hold on
plot(M,IncrPrior*10,'k--')
plot(M,M2prior*10,'r')
xlabel('Possible Values of the Mean')
ylabel('Prior Probability Density')
legend('Model 1 Prior','Increment Prior','Model 2 Prior')
title('Prior Distributions')
text(4,.3,{['The prior for Model 2 is the convolution'];...
    ['of the Increment prior and the Model 1 prior']})

%% Computing the integrated likelihoods
% For Model 1, we assume that we are drawing the Many data from the same
% distribution that we drew the Few data from. Thus, our prior for Model 1
% is the posterior probability distribution for the mean of the Few Group.
% For Model 2, we assume that we are drawing from a distribution t

for i=1:length(ManyData) % stepping through data, computing values of the
    % posterior likelihood of each datum (row variable is the datum) as a
    % function of the assumed value of the mean (column variable)
    
    LogLkFunction(i,:)=log10(normpdf(ManyData(i),M,6));
    % log likelihood of the datum, ManyData(i), for each of the possible
    % values of the mean (in the vector M), with the log of the normalizing
    % factor added
    
    Z1many(i,:)=log10(M1prior)+LogLkFunction(i,:); % Adding the log
    % likelihood to the log of the prior probability is equivalent to
    % multiplying the likelihood of a datum (for a given value of the mean
    % of the likelihood function) by the prior probability of that value
    % for the mean. Doing this for all the different values of the mean
    % gives the posterior log likelihood function for that datum. Thus, the
    % resulting array has a row for each datum and the columns give the
    % posterior log likelihood for that datum, assuming the Model 1 prior
    % probability distribution
    Z2many(i,:)=log10(M2prior)+LogLkFunction(i,:);
    % The posterior likelihood of each datum in the Many Group, using the
    % Model 2 prior
    
end % Of computing the likelihood array


Z1manySum=sum(Z1many); % sum across all the data in Many. This is the
% posterior log likelihood function for the mean of the Many group given Model
% 1's prior distribution function

L1many=10.^Z1manySum; % Taking antilogs to get the likelihoods given Model 1

Z2manySum=sum(Z2many); % This is the posterior log likelihood function for the
% mean of the Many group given Model 2

L2many=10.^Z2manySum; % Taking antilogs to get the likelihoods given Model 2

figure % plotting the posterior likelihood functions

plot(M,L1many,'k')

hold on

plot(M,L2many,'r')

plot([mean(FewData) mean(FewData)],ylim,'k--') % plotting the mean of the
% Few data

plot([mean(ManyData) mean(ManyData)],ylim,'r--') % plotting the mean of the
% Many data

BF=sum(L1many)/sum(L2many); % Bayes factor (ratio of the two intergrals)

xlabel('Possible Values of Many Mean')
ylabel('Posterior Likelihood')
title('Posterior Likelihood Functions for Mean of Many Group')
legend('Model 1','Model 2','Few Mean','Many Mean')
Top=ylim;
text(9,.95*Top(2),['Bayes factor = ' num2str(BF)])

% One can see at a glance that--with this prior for Model 2  the red (Model
% 2) likelihood function is much less than the area under the black
% likelihood function, meaning that the likelihood ratio (Bayes factor)
% favors Model 1. And, indeed, the Bayes factor is 4.33

% The modes (locations of the peaks) of the likelihood
% functions give the maximum likelihood estimates of the location of the
% means for the two models. The dashed lines show the sample means. If
% Model 1 is the right one, then the mean of the Many group lies to the
% right of the most likely "true" value for this mean, whereas if Model 2
% is true, then the it lies to the left of its "true" value. The fact that
% the sample means are not at the maxima of the posterior likelihood
% functions reflects the effect of the informative prior distributions on
% our estimates of where the mean value for all possible subjects in the
% Many condition really is. The posterior likelihood function for Model 2
% has less area than the posterior likelihood function for Model 1, because
% its prior is both more vague (diffuse, spread out) and located somewhat
% farther beyond where the data suggest that the mean really is. Thus, we
% can improve the posterior likelihood of Model 2 relative to Model 1 by
% making it less vague and by tailoring its prediction of the increment
% factor more closely to what the data imply

%% Loop to compute Bayes factor as a function of increment factor and CV

for i=1:10 % Stepping through possible values for the increment factor
    
    IncrFactor=1+.05*i; % Increment factors ranging from 5% to 50%
    
    IncrMean=IncrFactor*FewMode-FewMode;
 
    for s=1:10 % Stepping through possible values for the sharpness of
        % the weibull distribution, the higher s (the shape or sharpness
        % parameter
        % 
        
        Sharpness=10/s; % The higher the sharpness parameter, the narrower
        % and more sharply peaked the weibull distribution is. When
        % sharpness= 1, it has the shape of the exponential; when sharpness
        % equals 3.4, it approximates a normal distribution
        
        IncrPDF=wblpdf(M,IncrMean,Sharpness);
        
        IncrPrior=IncrPDF/sum(IncrPDF); % Normalizing the increment distribution
        
        M2prior=conv(IncrPrior,M1prior); 

        M2prior=M2prior(1:200);
        
        Z2many=repmat(log10(M2prior),16,1)+LogLkFunction; % By replicationg the
        % vector for the prior probability distribution 16 times, we make
        % an array of the same size as the array for the datum-by-datum
        % likelihood function
        
        Z2manySum=sum(Z2many); % summing across the data gives log likelihoods
        % for Model 2
        
        L2many=10.^Z2manySum; % Likelihoods, Model 2
        
        BF2(i,s)=sum(L2many)/sum(L1many); % Bayes factor
        
    end % of stepping throguh vagueness
    
end % of stepping through increments
%        

MxBF2=max(max(BF2)) % The maximum value of BF2 (the odds in favor of M2)

[MxBF2row,MxBF2col]=find(BF2==max(max(BF2))); % location of maximum

ContLvls=[.5*MxBF2 .1*MxBF2 .05*MxBF2 .01*MxBF2]; % contour levels at .5,.1
% .05, and .01 of the maximum odds in 

figure % contour plot of BF2

contour(BF2,ContLvls)

hold on

plot(MxBF2row,MxBF2col,'ok') % point of maximum relative likelihood

xlabel('Weibull S')

ylabel('Increment Factor (fraction of Few Mean)')

title({['Contour Plot of Bayes Factor for Model 2 vs 1'];...
    ['as Function of Parameters of Weibull Increment Prior']})

text(1.1,1.1,['Max BF = ' num2str(MxBF2)])

%%
figure % the increment prior that gives even odds for Model 2 vs Model 1

IncrMag=1.05*FewMode-FewMode; % The magnitude of the prior
% increment function that achieves equal likelihood with Model 1 is the
% magnitude that produces the highers BF times the best estimate
% of the mean of the Few Group minus the mean of the Frew Group

plot(M,wblpdf(M,IncrMag,10))

hold on
xlabel('Difference Between Many & Few Means')
ylabel('Prior Probability')
title('Increment Prior Required for Parity with No-Diff Model')

plot([IncrMag IncrMag],ylim,'--')

text(2,2.7,{['This is the increment in mean CR strength'];...
    ['that an alternative model must predict'];...
    ['to achieve odds parity with the No-Increment Model']})



% The most important thing to note is that essentially all of the
% probability mass of the parity increment prior lies within an increment
% of less than 1 response per minute on a base of about 8 responses per
% minute, which means that a viable model for the additional effect of the
% 8-fold increase in the number of trials must predict an increment no
% larger than 12%

%% Looking at how well the models fit the data

figure % Cumulative distribution of Many data and model fits

stairs(sort(ManyData),(1:length(ManyData))/length(ManyData),'r')
% cumulative distribution of Many data

hold on

[MxL1many,I_L1many]=max(L1many); % finding the column in the
% L1many vector that has the maximum value

ModeMany=M(I_L1many); % maximum likelihood value for the mean of
% the Many Group given Model 1. (This is also the new ml value for the mean
% of the Few Group under this model, because the model asserts that the
% data in both groups are draws from the same distribution, so, given that
% model, we have now made 24 draws from that distribution

Mm=linspace(0,30);

plot(Mm,normcdf(Mm,ModeMany,6),'r--') % Best fitting likelihood function
% given Model 1 (No-Diff)

plot(Mm,normcdf(Mm,1.05*FewMode,6),'r-.') % Best fitting likelihood function
% for the mode of the prior distribution for the Mean of the Many Group in
% the "yes-there-is-a-difference" model that achives odds parity with the
% No-Diff model. It would be a truly remarkable model that predicted a
% distribution that close to--but nonetheless distinguished from--the
% distribution expected under the no-difference model

xlabel('Responses Per Minute')
ylabel('Cumulative Fraction of Subjects')
title({['Comparison of Model 2''s Prior Mode'];['with Model 1''s Max Posterior Likelihood']})
legend('Many Data','Model 1 ML Fit','Fit Using Modal M2 prior')


%% Gamma increment prior

A=2; % Shape parameter of gamma

B=.5; % Scale parameter. Product AxB is the mean increment

IncrPDF=gampdf(M,A,B);

IncrPrior=IncrPDF./sum(IncrPDF);

 %figure
plot(M,IncrPrior,'k--') % This distribution expresses in probabilistic terms
% our prior beliefs (or at least conceivable assumptions we could make)
% about the size of the effect of the 8-fold increase in number of trials
% on the mean strength of the CR
xlabel('Mean Increment in CR Strength')
ylabel('Probability Density')
title('Gamma Increment Prior')
%%
% To get the prior distribution for the mean of the Many Group, we have to
% integrate our prior belief (or possible assumption) about the size of the
% effect of the ADDITIONAL trials (the increment distribution) with our
% uncertainty about the size of the effect of the initial trials (the
% trials that the Few and Many Groups both got). Our uncertainty about the
% size of the effect of those initial trials is specified by the
% distribution for the mean of the Few Group (NB, not the distribution from
% which the Few sample was drawn, but rather the distribution of the means
% of samples of that size). We put them together by convolving them

M2prior=conv(IncrPrior,M1prior); 

M2prior=M2prior(1:100); % The number of points in the function
% obtained by convolving two vectors is the sum of their lengths minus 1.
% Only the first 100 entries fall on the same time stretch as the two
% functions that were convolved. This trims the M2prior
%
figure % M1 prior, Increment prior and M2 prior
plot(M,M1prior,'k'); hold on
plot(M,IncrPrior,'k--')
plot(M,M2prior,'r')
xlabel('Possible Values of the Mean')
ylabel('Prior Probability Density')
legend('Model 1 Prior','Gamma Increment Prior','Model 2 Prior')
title('Prior Distributions')
text(4,.073,{['The prior for Model 2 is the convolution'];...
    ['of the Increment prior and the Model 1 prior']})
%%
clear A B Mean Sigma BF2
for i=1:20 % Stepping through possible values for the shape parameter of
    % the gamma 
    
    A(i)=.7+.3*i; % (1= exponential; the higher the value of A, the less the
    % skew
    
 
    for j=1:20 % Stepping through possible values for the shape.

        B(j)=.3*j/A(i);  % The product of A&B is the mean; thus, to keep the same
        % scale as A increases, one has to normalize by A.
        
        Mean(i,j)=A(i)*B(j);
        
        Sigma(i,j)=sqrt(A(i))/B(j);
                
        IncrPDF=gampdf(M,A(i),B(j));
        
        IncrPrior=IncrPDF/sum(IncrPDF); % Normalizing the increment distribution
        
        M2prior=conv(IncrPrior,M1prior); 

        M2prior=M2prior(1:100);
        
        Z2many=repmat(log10(M2prior),16,1)+LogLkFunction; % By replicating the
        % vector for the prior probability distribution 16 times, we make
        % an array of the same size as the array for the datum-by-datum
        % likelihood function
        
        Z2manySum=sum(Z2many); % summing across the data gives log likelihoods
        % for Model 2
        
        L2many=10.^Z2manySum; % Likelihoods, Model 2
        
        BF2(i,j)=sum(L2many)/sum(L1many); % Bayes factor
        
    end % of stepping through scale parameter
    
end % of stepping through shape parameter
%        

MxBF2=max(max(BF2)) % The maximum value of BF2 (the odds in favor of M2)

[MxBF2row,MxBF2col]=find(BF2==max(max(BF2))); % location of maximum

ContLvls=[.5*MxBF2 .1*MxBF2 .05*MxBF2 .01*MxBF2]; % contour levels at .5,.1
% .05, and .01 of the maximum odds in 

figure % contour plot of BF2

contour(BF2,ContLvls)

xlabel('B (scale)')

ylabel('A (shape)')

title('Contour Plot for BF2 With Gamma Increment Prior')

text(1.3,19.3,'1.03')

text(12,8,'Levels are .5, .1,.05, & .01 of maximum')

%%
figure % the increment prior that gives even odds for Model 2 vs Model 1

plot(M,gampdf(M,A(20),B(1)))

hold on
xlabel('Difference Between Many & Few Means')
ylabel('Prior Probability')
title('Gamma Increment Prior Required for Parity with No-Diff Model')

[Mx,R]=max(gampdf(M,A(20),B(1)));

plot([M(R) M(R)],ylim,'--')

text(2,2.7,{['This is the increment in mean CR strength'];...
    ['that an alternative model must predict'];...
    ['to achieve odds parity with the No-Increment Model']})

% The most important thing to note is that essentially all of the
% probability mass of the parity increment prior lies within an increment
% of less than 1 response per minute on a base of about 8 responses per
% minute, which means that a viable model for the additional effect of the
% 8-fold increase in the number of trials must predict an increment no
% larger than 12%