clear all
close all
%% this is a simple bit of code to explore the very basics of Bayesianism.

%% Simulate draws from linear model (for different X values)
numberofx=30; %% this is the number of measurements we have - the distinct X values
hyperparameterbeta=100; %% this captures how strong our prior belief is that alpha=bias
hyperparameterbiasonalpha=30; %% the a-priori expected value for alpha

dataY=cumsum(ones(numberofx,1))+randn(numberofx,1); % these are the simulated draws Y
wheresampleX=cumsum(ones(numberofx,1)); % these are the points of measurement X
figure; scatter(dataY,wheresampleX);

%% Write down likelihood and maximize it using a simple numerical optimization (perhaps fix sigma=1 and mu=0)

%% P(dataY|X,sig,mu,alpha)=\Pi_i (1/sig^2)* exp(-(dataY_i-\alpha X_i)^2/sig^2)


%% let's just minimize the log-likelihood.
logliknoconstants=@(alpha)sum((dataY-alpha*wheresampleX).^2);

initializeparametersearch=2;
[maximumlikelihoodalpha]  = fminsearch(logliknoconstants,initializeparametersearch)
%% (of course we could just to this analytically but the generic need for numerical optimization is a problem for both max likelihood and maximum aposteriori methods)

%% let's consider the posterior P(alpha|Y) \propto P(Y|X,sig,mu,alpha)P(alpha|beta,meanalpha) where P(alpha|beta,meanalpha) is the prior on alpha with the hyperparameters (parameter of the prior) beta and meanalpha
%% this adds an extra term to my likelihood.

logposteriornoconstants=@(alpha)(sum((dataY-alpha*wheresampleX).^2))+hyperparameterbeta*(hyperparameterbiasonalpha-alpha)^2;
[maximumaposteriorialpha]  = fminsearch(logposteriornoconstants,initializeparametersearch)

