%% Problem 1
fprintf('\nProblem 1 Results\n\n')

% Markov Chain Monte Carlo (MCMC) in 2D

N=10000;  %number of steps in  chain

%Intitialize chain vector
Chain=zeros(N,2);

%Setup initial value
X0=5;
Y0=5;
Chain(1,1)=X0;
Chain(1,2)=Y0;

%Setup probability function to explore
Mu_X=0;
Sigma_X=4;
Mu_Y=0;
Sigma_Y=10;
p_xy=@(X) normpdf(X(1),Mu_X,Sigma_X)*normpdf(X(2),Mu_Y,Sigma_Y);

%Setup test jump size
TestSigma=2;

%run chain

for nn=1:N-1
    %Make test jump
    CurrentX=Chain(nn,:);
    TestX(1)=CurrentX(1)+TestSigma*randn();
    TestX(2)=CurrentX(2)+TestSigma*randn();
    
    %Metropolis-Hastings rule
    if rand() < (p_xy(TestX)/p_xy(CurrentX))
        Chain(nn+1,:)=TestX;
    else
        Chain(nn+1,:)=CurrentX;
    end
end

% Histogram of x,y
subplot(2,1,1)
hist(Chain(:,1),50)
xlabel('x')
ylabel('p(x)')
subplot(2,1,2)
hist(Chain(:,2),50)
xlabel('y')
ylabel('p(y)')

% Visualize Chain

Xleft=-20
Xright=20;

Yleft=-20;
Yright=20;

Xvec=(Xleft:.1:Xright);
Yvec=(Yleft:.1:Yright);

Zvec=1:N;
FigH=figure

for nn=1:100:N
    figure(FigH)
    %plot chain
    subplot(2,1,1)
    plot3(Chain(1:nn,1),Chain(1:nn,2),Zvec(1:nn),'b.')
    hold on
    plot3(Chain(nn,1),Chain(nn,2),Zvec(nn),'ro','Linewidth',3)
    % axis([Xleft Xright 0 N])
    hold off
    xlabel('x')
    ylabel('y')
    zlabel('Chain Index')
    
    %plot histogram X
    subplot(2,1,2)
    hist3(Chain(1:nn,:),[20 20])
    xlabel('x')
    ylabel('y')
    zlabel('p(x,y)')
    pause(.01)
end


% Look at the histogram in 2D and its shadow
n=hist3(Chain,[20 20]);
% change the size of figure
fh=figure('units','normalized','outerposition',[0,0,1,1]);
figure(fh)

hold on
% Plot histogram in 3D
hist3(Chain);

% Make a shadow from histogram
n1=n';
n1(size(n,1) + 1, size(n,2) + 1) = 0;
xb = linspace(min(Chain(:,1)),max(Chain(:,1)),size(n,1)+1);
yb = linspace(min(Chain(:,2)),max(Chain(:,2)),size(n,1)+1);
h = pcolor(xb,yb,n1);
h.ZData = ones(size(n1)) * -500;
colormap(hot) % heat map
grid on
view(3);


%% Problem 2
fprintf('\nProblem 2 Results\n\n')
clear all
close all

%load Data
load HW11Data

%Plot given data
figure;
plot(T,Volts,'.b');
xlabel('Time');
ylabel('Volts')
legend('Expeimental data for T=A*V')

%The goal is to estimate the slope, 'A' along with the standard error on 'A'. 
Sigma_Noise=3;
Volts=T*3+Sigma_Noise*randn(size(Volts));

%find likelihood function of Mean(A)
Avec=linspace(0,5,1000);
LogL=zeros(size(Avec));

for ll=1:length(Avec)
   %See HW11_Solution_SlopeLogL
   [ LL ] = HW11_Solution_SlopeLogL(Avec(ll),Sigma_Noise,Volts,T); 
   LogL(ll)=LL; 
end

%Show the -LogL plot
figure
plot(Avec,LogL)
xlabel('A')
ylabel('Neg Log Likelihood')

%Show the likelihood plot
figure
plot(Avec,exp(-LogL))
xlabel('Mean(A)')
ylabel('Likelihood')

% Bayes rule:
% P(theta|data)=P(data|theta)*P(theta)/P(data)

% P(data) is integral of P(data|theta)*P(theta) over theta

% By Bayes rule, using a uniform prior(it is given in problem), P(theta)=1/L 
% where L is the size of interval of integration. 

L_A=exp(-LogL); %This is P(data|theta)

%Step size for numeric integral
d_Avec=Avec(2)-Avec(1);

%This is denominator in Bayes rule.  
NormInt=sum(d_Avec*L_A);

%This is the probability distribution of Mean_A given our data,
%P(theta|data)
P_AGData=L_A/NormInt;

figure
plot(Avec,P_AGData)
xlabel('A')
ylabel('P(A|Data)')

% We can approximate this shape using a 1D Gaussian
% the second derivative of a the log of a 1D Gaussian is 1/sigma^2

%The 'fminunc' function works like fminsearch, but also returns the 
% 'Hessian' - the matrix of second derivates.  For a scalar paremeter 
% like we have here, then Hessian is 1/sigma^2

%Setup fminunc and find optimum values. 
Options=optimset();
Theta0=2.5;
[Theta_Found,~,~,Output,Grad,Hess]=fminunc(@HW11_Solution_SlopeLogL,Theta0,Options,Sigma_Noise,Volts,T);

%We can give the error on our estimate as the standard error, which is the 
%sigma of our Gaussian distribution
StandardError=sqrt(inv(Hess))

%Plot the distrubution and the Gaussian approximation together
LHess=normpdf(Theta_Found,Avec,StandardError)
hold on
plot(Avec,LHess,'--','linewidth',2)
legend('P(A|Data)','Gaussian From Hessian')

fprintf('The parameter A would be \n\n A=%g +/-%g\n\n', Theta_Found,StandardError)
% Now look at this A as a slope of the given data
A_Found=Theta_Found;
Voltsfit=A_Found*T;
fprintf('Look at the figure to compare found A with data !\n\n')
figure 
plot(T,Volts,'b.');
hold on
plot(T,Voltsfit,'r-');
xlabel('Time')
ylabel('Volts')
legend('Experimental data','Curve fit')



%% Problem 3
fprintf('\nProblem 3 Results\n\n')


HW11_Solution_GUI();



%% Problem 4
syms a b c d
A = [a,b;c,d]; %A doesn't need to be symbolic. It's just a matrix of symbolics
%Print nicely
fprintf('Inverse of:\n');
pretty(A);
fprintf('is\n');
pretty(inv(A))

%% Problem 5
syms x mu
%Print nicely
fprintf('Integral of e^(x^2) over all x:\n');
f=int(exp(-mu^2),-inf,inf)
fprintf('Pretty Version: \n');
pretty(f)


