% This script will visually illustrate how Bayesian inference works within the context of some simple tasks, all of which are related to the flipping of coins.  (Somewhat) Specifically, if we have some collection of coins, each indexed by j, we will let H_j be a variable that indicates the probability that coin j will gives heads on any given flip.
% In this first scenario, suppose that you are told that there are M total coins, all of which are visually indistinguishable, but you are also provided with the value of H_j for each coin (let's suppose they are all different).  You are asked to reach into a bag and grab one coin at random, then you are allowed to flip that coin as many times as you like in an attempt to determine which coin it is that you randomly selected.  How should you do this?
% One method is to use Bayesian inference.  Given the information we have stipulated in this problem, before making any flips of our chosen coin we know what the probability of us having chosen any specific coin k is: P\left(\textrm{coin}\;k\right)=\frac{1}{M}.  We will call this the prior probability of having coin k; "prior" here can be interpreted as the information we have prior to any further information being given to us (such further information would be, for example, how many heads come up in some number of tosses of this coin).  But, this scenario allows us to gather additional, useful data by flipping the coin.  We know that, if we were to toss coin k a number N times, the probability of observing a specific number of heads is related to H_k, and is actually a binomial distribution:
%  P(heads|coin~k)= {N\choose heads}H_k^{heads}(1-H_k)^{N-heads}
% Note that this is a conditional probability: it is the probability of seeing this many heads out of N tosses, conditional on the coing being coin k).
% We can use the prior probability and the conditional probability noted above to determine the posterior probability (the probability after we observe addition data beyond what the prior distribution assumed) of the coin being k through Bayes' rule:
% P(coin~k|heads)=\frac{P(heads|coin~k)P(coin~k)}{P(heads)
% The denominator above is the overall probability of observing that specific number of heads across all coins.  Its role in this equation is really just to serve to normalize things, so that in the end the sum
% \sum_k P(coin~k|heads)=1.
% As is often the case when using Bayes' rule, then, the denominator does not require too much thought to calculate: just compute the numerators for all possible prior cases k, then the sum of all of those numerators is the denominator.
% Operationally, then, we can flip our coin N times, record the number of heads, and use the above equation to specify the probability that our coin is coin k, for all coins.  We might suspect that, as the number of flips goes to infinity, then eventually we will end up with a scenario where all of the posterior probabilities approaches 0 except for that of the true coin, which will approach 1.
% Let's set up some code that will simulate all of the above:

M=10; %how many coins are there total
%H=sort(rand(M,1)); %the H values for each coin, just chosen randomly for this example

%truecoin=randi(M); %this is the coin you selected from the bag, and you don't know this value
N=1000; %how many times you will flip your coin
heads=binornd(N,H(truecoin)) %the number of heads you observe

prior=ones(M,1)/M; %prior probability of coin k, the same for all k in this example
conditional=(H.^heads).*((1-H).^(N-heads)); %the conditional probability for each coin; note the binomial coefficient is missing, as it will cancel out with the denominator term so is unneeded
posterior=conditional.*prior;
posterior=posterior/sum(posterior);
plot(H,posterior,'ro')
ylim([0 1]);
xlabel("Possible H values (coins)");
ylabel("Probability of given H value (coin)");
[~,maxposterior]=max(posterior)
truecoin

% Try playing around with the above for different values of M and N.  Note that the coin that has highest posterior probability will not always be the true coin (especially at small N, and when some the true coin has a very similar H to another coin), though it should be if you let N be very large.
%
% One neat thing about Bayes' rule is that you can either lump all of your observations into one overall observation (like we did above, with just the total number of heads), or you can perform sequential updates after each observation one at a time (or in batches of whatever number you like).  If you are performing these sequential updates, the posterior from the previous observation just becomes the prior for the next one!  We can visualize this for our problem: instead of doing a single Bayes' rule with N observations and a total number of heads, let's do N sequential Bayes' rule updates, each of which corresponds to the results of a single flip and uses the posterior from the last observation as the prior for this one.  The only (minor) difference in this case is that, since each observation is a single coin flip, the probability P(heads|H) is Bernoulli distributed:
% P(heads|coin~k)=H_k^{heads}(1-H_k)^{1-heads}
% In the equation above, we should just note that, since this is only a single flip, the value of heads will either be 1 or 0.

M=10; %how many coins are there total
H=sort(rand(M,1)); %the H values for each coin, just chosen randomly for this example

truecoin=randi(M); %this is the coin you selected from the bag, and you don't know this value
N=100; %how many times you will flip your coin
prior=ones(M,1)/M; %prior probability of coin k, the same for all k in this example
gprior=prior; %the original prior
totheads=0; %to track the total number of heads across flips

for flips=1:N
    heads=binornd(1,H(truecoin)); %heads or not on this flip
    totheads=totheads+heads;
    conditional=(H.^heads).*((1-H).^(1-heads)); %the conditional probability for each coin; note the binomial coefficient is missing, as it will cancel out with the denominator term so is unneeded
    posterior=conditional.*prior;
    posterior=posterior/sum(posterior);
    plot(H,gprior,"b",H,posterior,"g"); % visualizing
    ylim([0 1]);
    legend("Prior","Current Posterior");
    xlabel("Possible H values (coins)");
    ylabel("Probability of given H value (coin)");
    prior=posterior;
    pause(.00001)
end
%now compute as above with a single Bayes' rule and totheads
conditional=(H.^totheads).*((1-H).^(N-totheads)); %the conditional probability for each coin; note the binomial coefficient is missing, as it will cancel out with the denominator term so is unneeded
gposterior=conditional.*gprior;
gposterior=gposterior/sum(gposterior);
hold on
plot(H,gposterior,"ro");
legend("Prior","Current Posterior","Single Update Posterior");
hold off

% Again, try messing around with the values above.
%
% A related example to the one above is to determine whether or not a particular coin is "fair".  More precisely, here we will attempt to determine what the probability H is of this coin to land on heads, given some number of observations we have made on N individual flips.
% At first, let's just compute this probability for some discrete set of possible H values:

M=100; % the number of discrete H values we will allow for
H=linspace(0,1,M)'; % the actual discrete H values we will allow

% For this example, we will first choose an uninformative prior distribution for H. That is, an initial probability distribution that really doesn't provide any indication of what values for H are truly likely to be correct.  In this problem, the clear choice for an uninformative prior is just a uniform distribution: as far as we know going in, all values of H are equally likely.

gprior=ones(1,M)'/M; % a completely uninformative, uniform prior
plot(H,gprior,"b"); % visualizing
ylim([0 1]);
legend("Prior");
xlabel("Possible H values");
ylabel("Probability of given H value");

%
% Now, we may not know what the true value of H is, but the computer needs to know this to proceed, so we'll set that here, randomly choosing some value between 0 and 1, uniformly.

secretH=rand();

% Alternatively, you could enter a value that you would like to use for H here (if you want to use a random H, be sure this slider is set to 0):

userH=0.25;

%
% Suppose that this coin is now flipped N times.  Pick a (reasonable!) N value here:

N=100;

%

if userH==0 % set the true value of H to use for this coin
    trueH=secretH;
else
    trueH=userH;
end
flips=randsrc(N,1,[1,0;trueH,1-trueH]); % generate many random "flips"
heads=sum(flips); % the number of "heads" observed on these flips
fprintf('Your virtual coin came up heads %g times out of %g flips.',heads,N)

%
% Now, we will use the information on the total number of heads observed to determine our posterior distribution for the value of H via Bayes' rule, just like above:
% P(H|heads)=\frac{P(heads|H)P(H)}{P(heads)
% Finally, we note again that P(heads|H) is a binomial distribution, with N total trials and probability of success on any given trial being H:
% P(heads|H)= {N\choose heads}H^{heads}(1-H)^{N-heads}
% Let's use this to compute the posterior, then also say some interesting things about it:

gposterior=H.^heads.*(1-H).^(N-heads).*gprior; % unnormalized Bayes' rule
gposterior=gposterior/sum(gposterior); % now normalize
plot(H,gprior,"b",H,gposterior,"r"); % visualizing
ylim([0 1]);
legend("Prior","Posterior");
xlabel("Possible H values");
ylabel("Probability of given H value");
meanpostH=gposterior'*H; % the mean of the posterior
cumu=cumsum(gposterior); % finding a 90% confidence interval
lower=max(H(cumu<.05));
upper=min(H(cumu>.95));
fprintf('The mean value of your posterior distribution for H is %g.',meanpostH)
fprintf('There is a 90 percent chance that H is between %g and %g.',lower,upper)

%
% Mess around with different values for N, and see how the 90% confidence interval varies, as well as how close the mean gets to the true value.
%
%
% As above, we can also compute the posterior after each individual flip observation.  Let's do this, and animate the resulting sequence of posterior distributions we obtain as we run through the sequence of all N observations:

prior=gprior; % starting from our original, uninformative prior
for iter=1:N
    flip=flips(iter); % consider only a single flip
    posterior=H.^flip.*(1-H).^(1-flip).*prior; % the posterior after that single flip
    posterior=posterior/sum(posterior); % normalizing the posterior
    prior=posterior; % set up this posterior as the prior for the next observation
    plot(H,gprior,"b",H,gposterior,"r",H,posterior,"g"); % visualizing
    ylim([0 1]);
    legend("Prior","Full Posterior","Current Posterior");
    xlabel("Possible H values");
    ylabel("Probability of given H value");
    pause(.00001)
end

%
% So far in this example, we have been computing everything numerically using a "large" number of discrete possible H values to evaluate the prior and posterior at.  However, in general we would like to allow H to be any real number on [0,1].  From an analytical point of view, this is no problem, as we just use Bayes' rule as usual, and our prior will just be some function of the continuous variable H instead of the discrete values we've used above.  The denominator from Bayes' rule will now generally be an integral instead of a sum, but still just really serves the purpose to normalize everything.
% When doing things this way, we can of course use whatever prior distribution P(H) we would like.  But, in some cases of Bayes' rule we can leverage the idea of using a conjugate prior to make things a little simpler for ourselves.  In order to define and discuss conjugate priors, let us first rewrite Bayes' rule in a more generic way that doesn't necessarily reference coins:
% P(A|B)=\frac{P(B|A)P(A)}{P(B)}
% Suppose we are working on a Bayes' rule problem with a particular functional form for the conditional probability P\left(B|A\right); in our coin example, this is binomial.  A conjugate prior for that conditional probability is any prior P\left(A\right) such that, upon multiplying by the conditional probability, the result (which is the unnormalized posterior) is of the same functional form as P\left(A\right); albeit probably with different values for any parameters (called hyperparameters) that P\left(A\right) might contain.  To hopefully make things more clear, let's again consider our coin example where the conditional proability is binomial.  It turns out that a conjugate prior to the binomial distribution is the Beta distribution:
% Beta(H;\alpha,\beta)\propto H^{\alpha-1}(1-H)^{\beta-1},~\alpha>0,~\beta>0 ;
% this distribution has two hyperparameters, \alpha and \beta, which dictate the overall shape of the distribution, along with things like its mean, variance, etc.  Note that for \alpha =\beta =1, the Beta distribution is just the uniform distribution, so our example above is already using a prior of this form!  We can see why this is a conjugate prior to the binomial distribution by multiplying the two:
% Binom(heads|H;N)Beta(H;\alpha,\beta)\propto H^{heads}(1-H)^{N-heads}H^{\alpha-1}(1-H)^{\beta-1}=H^{\alpha+heads-1}(1-H)^{\beta+N-heads-1}=Beta(H;\alpha+heads,\beta+N-heads)
% So, we see that the posterior is also a Beta distribution, but with new hyperparameters that depend on the prior hyperparameters we use as well as the observations that are made.
% Conjugate priors are very useful, as they vastly simplify Bayesian calculations.  This is especially true when making many sequential Bayesian updates, as using a conjugate prior to start out with will guarantee that at every step along the way we will end up with our prior/posteriors also being of that same form, and all we need do is update the hyperparameters correctly.
% Let's explore this in some code:

alpha=1; %initial hyperparameter alpha for the Beta conjugate prior
beta=3; %initial hyperparameter beta for the Beta conjugate prior
x=0:.001:1;
conjprior=betapdf(x,alpha,beta); %the conjugate prior
postalpha=alpha+heads; %updated hyperparameter alpha after observations
postbeta=beta+N-heads; %updated hyperparameter beta after observations
conjposterior=betapdf(x,postalpha,postbeta); %the posterior, which is also beta distributed

%let's now compute things the less efficient way that we used above
gprior=betapdf(H,alpha,beta);
gposterior=H.^heads.*(1-H).^(N-heads).*gprior; % unnormalized Bayes' rule
gposterior=gposterior/sum(gposterior); %normalize
dH=H(2)-H(1); %this is used to turn our gposterior above into an approximate probability density function
gposterior=gposterior/dH;

plot(x,conjprior,x,conjposterior,H,gposterior,'ro')
legend("Beta Prior","Beta Posterior","Discrete Posterior as PDF");
xlabel("Possible H values");
ylabel("Probability density of given H value");

% Play around with different values of the hyperparameters (just be sure they are positive).  For reference, the mean of the Beta distribution is \frac{\alpha}{\alpha+\beta} and the variance is \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}.  If you choose initial hyperparameters that put very little probability in the region of trueH, your end result might not be very good.