%% Lab 7 - Bayesian Models
%% Exercise 1: numerical calculation of the posterior
% The stimulus distribution P(S)is a Gaussian with a mean of 20 and standard deviation of 4. The measurement distribution P(X|S) is Gaussian with standard deviation of 5. A Bayesian observer infers S from an observed measurement X_{obs} = 30. We are now going to calulate the poserior pdf explictly.
% The vector of hypothesized stimulus values S_{hyp}: 0, .2, .4,...40

sHyp = 0:.2:40;

xObs = 30;
sdMeas = 5;
meanSig = 20;
sdSig = 4;

% Calculate the likelihood function and the prior on this vector of S_{hyp} values, remember the formula for a Gaussian is f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}

%--------------your code here--------------------
prior = [];
likelihood = [];

% Calculate the posterior, normalize and convert to a probability density:

protoPos = prior .* likelihood;     %multiply to get posterior
normPos = protoPos/sum(protoPos);   %normalize by the sum
probDens = normPos/.2;              %divide by step size to get the density

% Calculate the mean of the posterior \mu_{pos} = \frac{\frac{\mu_1}{\sigma^2_1}+\frac{\mu_2}{\sigma^2_2}}{\frac{1}{\sigma^2_1}+\frac{1}{\sigma^2_2}}, \sigma_{pos} =\sqrt{ \frac{1}{\frac{1}{\sigma_1^2}+\frac{1}{\sigma_2^2}}}

%--------------your code here--------------------
meanPos = [];
sdPos = [];

% Plot the likelihood, prior and posteior in the same plot. Is the posterior wider or narrower than likelihood and prior?

figure
plot(sHyp,prior,'r','LineWidth',2)
hold on
xline(meanSig,'r--','LineWidth',2)
plot(sHyp,likelihood,'b','LineWidth',2)
xline(xObs,'b--','LineWidth',2)
plot(sHyp,probDens,'m','LineWidth',2)
xline(meanPos,'m--','LineWidth',2)
text(5,.1,['\mu Pos: ' num2str(meanPos)])
text(5,.08,['\sigma Pos: ' num2str(sdPos)])
box off
set(gca,'TickDir','out')
ylabel('Probability')
xlabel('Stimulus Intensity (a.u.)')
xlim([.5 40.5])
xticks([0:10:40])
legend('Prior','\mu prior','Likelihood','\mu like','Posterior','\mu pos')
title('Original plot)')

% Change the standard deviation of the measurement distribution to a very large value. What happens to the posterior?

%--------------your code here--------------------
sdMeas = []; %choose a number
likelihood = []; %same as above!
meanPos = []; %same as above!
sdPos = []; %same as above!

protoPos = prior .* likelihood;     %multiply to get posterior
normPos = protoPos/sum(protoPos);   %normalize by the sum
probDens = normPos/.2;              %divide by step size to get the density

figure
plot(sHyp,prior,'r','LineWidth',2)
hold on
xline(meanSig,'r--','LineWidth',2)
plot(sHyp,likelihood,'b','LineWidth',2)
xline(xObs,'b--','LineWidth',2)
plot(sHyp,probDens,'m','LineWidth',2)
xline(meanPos,'m--','LineWidth',2)
text(5,.1,['\mu Pos: ' num2str(meanPos)])
text(5,.08,['\sigma Pos: ' num2str(sdPos)])
box off
set(gca,'TickDir','out')
ylabel('Probability')
xlabel('Stimulus Intensity (a.u.)')
xlim([.5 40.5])
xticks([0:10:40])
legend('Prior','\mu prior','Likelihood','\mu like','Posterior','\mu pos')
title('Large measurement SD')

% Change the standard deviation of the measurment distirbution to a very small number. What happens to the poserior now?

%--------------your code here--------------------
sdMeas = []; %choose a number
likelihood = []; %same as above!
meanPos = []; %same as above!
sdPos = []; %same as above!

protoPos = prior .* likelihood;     %multiply to get posterior
normPos = protoPos/sum(protoPos);   %normalize by the sum
probDens = normPos/.2;              %divide by step size to get the density


figure
plot(sHyp,prior,'r','LineWidth',2)
hold on
xline(meanSig,'r--','LineWidth',2)
plot(sHyp,likelihood,'b','LineWidth',2)
xline(xObs,'b--','LineWidth',2)
plot(sHyp,probDens,'m','LineWidth',2)
xline(meanPos,'m--','LineWidth',2)
text(5,.1,['\mu Pos: ' num2str(meanPos)])
text(5,.08,['\sigma Pos: ' num2str(sdPos)])
box off
set(gca,'TickDir','out')
ylabel('Probability')
xlabel('Stimulus Intensity (a.u.)')
xlim([.5 40.5])
xticks([0:10:40])
legend('Prior','\mu prior','Likelihood','\mu like','Posterior','\mu pos')
title('Small measurement SD')

%
%% Exercise 2: conditional and marginal probability (prosecutor fallacy)
% You are receiving a test for a relatively rare disease. The test is fairly accurate--if you have the disease, the test will most likely (sensitivity = 0.9) report a positive result, and if you don't have it, the test will most likely (specificity = 0.8) report a negative result. However, in the absence of a test, you are unlikely to have the disease, as only 2% of the population has it. Now suppose you take a test, and the result is positive. How likely are you to actually have the disease? Solve this in two ways:
% 1. Simulate this scenario for many trials, and count the number of times the positive subjects have the disease relative to total positive counts.
% 2. Calculate the probability theoretically using Bayes' rule.

%Variables for the problem:
sensitivity = 0.9;  % probability of positive test given disease
specificity = 0.8;  % probability of negative test given no disease
prevalence  = 0.02; % probability of having the disease 
n           = 1e5;  % number of test subjects

%
%% Exercise 2: method 1 - Simulation
% Simulate this scenario for n trials, and count the number of times the positive subjects have the disease relative to total positive counts.

outcome_table = zeros(2,2); %counts for having the disease.

%--------------your code here--------------------

%
%% Exercise 2: method 2 -  Calculate the probability theoretically using Bayes' rule.
% Below is the formula for Bayes' rule:
% P(A|B) = \frac{P(B|A)P(A)}{P(B)
% and in the case of this problem:
% P(disease|positive) = \frac{P(positive|disease)P(disease)}{P(positive)}


%Assign the variables to caclulate the outcome:

%--------------your code here--------------------

% Are these values bigger or smaller than you would have thought? Can you explain why? Are the conditions of having the disease and not having the disease symmetric? What is the primary asymmetry about these outcomes?
%
% Bonus exercise: What about the probability you *don't* have the disease if you test negative? Is this also small? Is the test just worse than advertised?

%--------------your code here--------------------

%
%% Exercise 3: binomidal
% Suppose we collected N flips of a coin with true probability p of coming up heads. We count the number of heads we actually observe, n. Compute the mean and standard deviation of n by 3 different methods:
% 1. The "Look it up" method--use "binom" in the scipy.stats module and create a distribution with given parameter values, then access the means and standard deviations computed for you.

%--------------your code here--------------------

% 2. Manually compute the probability distribution P(n) over the number of heads n ∈ {0, 1, ..., N} and compute the mean and standard deviation from that.

%--------------your code here--------------------

% 3. Compute the mean and standard deviation theoretically using the mean and std for an individual coin flip and considering the binomial as a sum over many independent coin flips.

%--------------your code here--------------------

%
%% Exercise 4: Bayesian updating (Beta and Bernoulli distribution)
% What if we don't know p? Can we infer p based on coin flip outcomes? Can we be certain based on finite data? What if we simply know that our coin is fair, before observing any flips? (E.g. maybe we can just tell by looking at it the coin is perfectly symmetric, and no amount of data will change our minds.) Does our interpretation of the data change? What if, on the other hand, the "coin" being tossed is some other physical process that we could easily imagine is not 50-50?
% The Bayesian approach is to treat p as itself a random variable. One might argue philosophically that p is fixed but simply unknown, and it's a fallacy to model p as the outcome of some random variable. The Bayesians disagree, and they're fine with treating "uncertainty" and "randomness" with the same math---probability distributions over possibilities---and even having them interact! The coin flips are obviously "random," whereas the parameter p may simply be "unknown," but we can still model our uncertainty about p using a probability distribution over the interval [0, 1].
% The point of the first paragraph is that our prior probability of p matters. If we were 100% confident the coin is 50-50 (say based on physical measurements), we could flip 10 heads in a row and be amused at how unlikely that outcome was---yet still believe as certainly as before that p = 0.5. What that if, on the other hand, we went in with an open mind? What if our prior probability was that p could be any value between 0 and 1 with equal probability? Then how does observing some sequence of heads and tails update our initial belief? Let's start with one coin flip:
% Prior : P(p) = 1
% Likelihood : P(H|p) = p
% Posterior : P(p|H) = \frac{P(H|p)P(p)}{P(H)}
%                  = p\frac{1}{\int P(H|p)P(p) \,dp.\
%                    =\frac{p}{\int p \,dp} = 2p
% The Beta distribution was designed to handle this specific inference problem. The Beta distribution has 2 parameters, α and β. By definition, the distribution Beta_{α,β}(p) gives the probability of p given the observation of α heads and β tails, assuming you started with a flat prior distribution on p.
% For this exercise, generate iteratively individual Bernoulli outcomes from a distribution with some (ground-truth) p parameter. As you get each outcome, update the posterior distribution over p incorporating each new outcome, by two methods and verify they get the same answer:
% 1. Taking a Beta distribution out of the box using betapdf() and changing the α and β values.
% 2. Using Bayes rule on the previous posterior in light of the new observation.
%
%% Exercise 4: Method 1 - using betapdf()

p_axis = linspace(0,1,101);  % calculate the distribution over x = 0,.01,.02,...,1
Beta   = betapdf(p_axis, 1, 1);   % compute the beta "prior" distribution
figure
plot(p_axis, Beta); hold on

%ground truth p
p         = 0.3; %you can change this value!
N_samples = 8; %play with this value as well
alpha     = 1;
beta      = 1;
lgd       = {'Prior'};
for i = 1:N_samples
    outcome = binornd(1, p);
    
    %Update alpha and beta, get new Beta pdf
    
   %--------------your code here--------------------
    
    lgd{end+1} = ['\alpha = ', num2str(alpha), ', \beta = ', num2str(beta)];
    gray_scale = 1- i/N_samples;
    plot(p_axis, Beta, 'Color', ones(1,3).*gray_scale); hold on;
end
Lgd = {lgd{:}};
hold off; legend(Lgd); box off;

%
%% Exercise 4: Method 2 - updating posterior using Bayes rule

dp = p_axis(2) - p_axis(1);

%really this is the "prior" but we'll be recursively updating it
posterior = ones(1,length(p_axis));

figure
plot(p_axis, posterior); hold on;

%ground truth p
p = 0.3; %you can change this value!
N_samples = 8; %play with this value as well
heads = 0;
tails = 0;
lgd  = {'Prior'};
for i = 1:N_samples
    outcome = binornd(1, p);
    
    %Update posterior using Bayes rule
    
    %--------------your code here--------------------
    
    lgd{end+1} = ['Heads = ', num2str(heads), ', Tails = ', num2str(tails)];
    gray_scale = 1- i/N_samples;
    plot(p_axis, posterior, 'Color', ones(1,3).*gray_scale); hold on;
end
Lgd = {lgd{:}};
hold off; legend(Lgd); box off;

%
%% Exercise 5: Bayesian updating (Normal distribution)
% The above exercise shows that the "conjugate" distribution to the binomial is the beta, since our posterior of the bionimal parameters given binomial observation follows a beta distribution. In this exercise we will show the conjugate of the normal distribution is... the normal distribution!
% EXERCISE: Play with these values, examining what happens as the sample deviates more and more from the mean of the prior, and as the sample size increase

% Pick a mean and standard deviation for your prior distribution
priormean = 0;
priorvar = 1;

% Pick a set of parameters from which you will draw your samples (number % 
%of samples and mean). You might as well use the common variance, since 
%that's what we are assuming
N = 20;
datamean = 3;

% Draw a sample
priorsd = sqrt(priorvar);
sample = datamean + priorsd*randn(N,1);

% Compute the prior over a fine, discrete set of values that is wide
% enough to encompass all samples and most of the probability mass
x = -5:.01:7;
nsample = length(x);
prior = (1/(sqrt(2*pi)*priorsd))*exp(-(x-priormean).^2/(2*priorvar));

% Compute the posterior using Bayes Rule, normalizing both prior and 
% posterior to sum to one over the discrete points
deltaxs = (sample * ones(1,nsample)) - (ones(N,1)*x);
probvalues = (1/(sqrt(2*pi)*priorsd))*exp(-deltaxs.^2/(2*priorvar));
likelihood = prod(probvalues);
posterior = likelihood.*prior;

prior = prior/sum(prior);
posterior = posterior/sum(posterior);

% Plot prior and posterior
figure;
plot(x,prior,'LineWidth',2);
hold on;
plot(x,posterior,'LineWidth',2);
box off;
legend('Prior','Posterior')

%