
% Download 'Data.mat' or 'Data.csv' from: 
% http://panda3.phys.unm.edu/~klidke/TestData/

load('Data')
imshow(Data,[],'InitialMagnification','fit')

%% Gain offset correction
Gain=20
Offset=100

D=(double(Data)-Offset)/Gain

%% Create a Model Generating function for 2D Gaussian PSF model

% Our fit parameters
% Theta = [X,Y,Photons,Bg]

%Microscope PSF width (pixels)
PSFSigma=1.3

% In case we don't have 'normpdf' from stats toolbox
my_normpdf=@(X,MU,Sigma)1/sqrt(2*pi*Sigma.^2)*exp( - (X-MU).^2/(2*Sigma.^2))

% Make anonymous function for model
[XGrid,YGrid]=meshgrid(1:8,1:8)  %Note 1-based grid
my_Model=@(Theta)Theta(4)+Theta(3)*my_normpdf(XGrid,Theta(1),PSFSigma).*my_normpdf(YGrid,Theta(2),PSFSigma)

% Test
Theta=[4,4,1000,5]
M=my_Model(Theta)
imshow(M,[],'InitialMagnification','fit')

%% Making a function to calculate the likelihood ratio used in MH 

% Direct approach
% LR = @(D,MC,MT) prod(exp(MC-MT).*(MT./MC).^D); 

% Easier to use log
% LLR = @(D,MC,MT) sum(sum(MC-MT+D.*log(MT./MC)));

% And then exponentiate
LR = @(D,MC,MT) exp(sum(sum(MC-MT+D.*log(MT./MC))));

%% The MCMC method using Metropolis-Hastings (MH)

% Setup Chain

NTrials=10000;
JumpSigma=[.1,.1,10,.1] % size scale of jumps matched to parameters

Chain=zeros(NTrials,4);
Chain(1,:)=[4,4,1000,5] %Starting Guess


for nn=2:NTrials
    ThetaC=Chain(nn-1,:);  %Theta current
    ThetaT=ThetaC+JumpSigma.*randn(1,4); %Theta test
    ModelC=my_Model(ThetaC);
    ModelT=my_Model(ThetaT);
    Alpha=LR(D,ModelC,ModelT); %accept based on likelihood ratio
    if rand()<Alpha %accept
        Chain(nn,:)=ThetaT;
    else
        Chain(nn,:)=ThetaC;
    end
end

% Show chain on top of data
figure
imshow(D,[],'InitialMagnification','fit')

hold on
X=Chain(:,1);
Y=Chain(:,2);
plot(X,Y,'g.')

%% Look at Chain 

NBurnin=1000;

XNB=Chain(NBurnin:end,1);
YNB=Chain(NBurnin:end,2);
PNB=Chain(NBurnin:end,3);
BGNB=Chain(NBurnin:end,4);

% Histogram of 'X'.  These can seen as samples from the posterior
% probabability distribution P(X|Data)
figure
hist(XNB,50)

%% Look at all paramters and standard errors. 

X_Found=mean(XNB)
X_SE=std(XNB)

Y_Found=mean(XNB)
Y_SE=std(XNB)

P_Found=mean(PNB)
P_SE=std(PNB)

BG_Found=mean(BGNB)
BG_SE=std(BGNB)

%% Diffusing particles 2 Frames

N=3
X=.5*randn(1,3)
Y=.5*randn(1,3)

figure
plot(X,Y,'ro')

D=.1 %um^2/s
DT=.1;

Sigma_D=sqrt(2*D*DT);

X2=X+randn(1,3).*Sigma_D;
Y2=Y+randn(1,3).*Sigma_D;

hold on
plot(X2,Y2,'go')
xlabel('X (micron)')
ylabel('Y (micron)')

%%  Simple Brownian Motion Generator

Traj=cumsum(Sigma_D*randn(1000,2))
figure
plot(Traj(:,1),Traj(:,2))














