clear all
GAL1pmodstruct;
bestfitparam;

k(1) = 1.769292045467269;
k(2) = 20.589996419308118;

%interpolate as chebfun for speedup
prior2 = chebfun(@(x)nctpdf(x,k(1),k(2)),domain(0,2e3));


param={'mu1','k1','alpha','g1','k2','g2','sigma1','sigma2'};

fitDataE = [ dE(:,1); dE(:,2) ].';
fitSTDE  = [ stdE(:,1); stdE(:,2) ].';

mini                        = 1e-5*ones(1,length(p));
maxi                        = 1e5*ones(1,length(p));

% Include priors?
%mStruct.p.mu1=1;
%mStruct.p.k1=2;
%mStruct.p.alpha=3;
%mStruct.p.g1=4;
%mStruct.p.Kg=5;
%mStruct.p.k2=6;
%mStruct.p.g2=7;
priors = 0;

fun = @(par)residGAL1p( t, x0, par, fitDataE, fitSTDE, input, prior2 );


initPars                    = p;

opts = optimset( 'TolX', 1e-12, 'TolFun', 1e-6, 'MaxFunEvals', 1e4, 'MaxIter', 1e4, 'Jacobian', 'off', 'Display', 'off' );


[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );

initPars

kprev = initPars(1:6);
sigma = initPars(7:8);

CovMat = [
 24000,      0,      0,      0,      0,      0;
     0, 0.0058, 0.0005,      0,      0,      0;
     0, 0.0005, 0.0022,      0,      0,      0;  
     0,      0,      0,     50,      0,      0;
     0,      0,      0,      0, 160000,      0;
     0,      0,      0,      0,      0,     20;
     ];

kc=0.0035;
acc=0;
nacc=0;
n=1000000;
sample=zeros(n,6);
prob=zeros(n,1);
lh=zeros(n,1);
accepted=zeros(n,1);

            p = [kprev sigma];
            plotGAL1p
            drawnow

%h=waitbar(0,'MCMC');

k=1;
while (k<n+1)
    if(mod(k,500)==0)
        (k)/n
        display(['Acceptance Rate: ' num2str(acc/(nacc+acc)*100) '%'])
    end
    try
    % the step is sampled from a multivariate normal distribution and
    % scaled with kc
    kcur = kprev + kc*mvnrnd(zeros(1,6),CovMat);
    sample(k,:) = kcur;
    % calculate the a-posteriori probability of the new sample
    prob(k) = -sum(fun([kcur sigma]).^2);
    if k>1
        % check whether we reject the sample or not. mind that the
        % probability is log-scaled
        if log(rand(1)) < +prob(k)-prob(k-1)  
            % update current sample
            kprev = kcur;
            % count accepted samples
            acc=acc+1;
            % update waitbar
            %waitbar((k-1)/n)
            % save current acceptance rate for post-processing
            lh(k)=acc/(nacc+acc)*100;
            k=k+1;
            %p = [kcur sigma];
            %plotETOHPROM
            %drawnow
        else
            % count not accepted samples
            nacc=nacc+1;
        end
    else
        % we need to do something else for the first sample.
        if log(rand(1)) < prob(k)+RN
            kprev = kcur;
            acc=acc+1;
            %waitbar((k-1)/n)
            lh(k)=acc/(nacc+acc)*100;
            k=k+1;
        else
            nacc=nacc+1;
        end  
    end

    catch ME
        disp(ME)
    end
end
%close(h)

%%
% Thinning
thinstep=100;

tsample=sample(1:thinstep:end,:);


%% Post Processing 
%

%%
% Heat Map of the joint sample density of two parameters
%
figure(7)
for m = 1:length(kprev)^2
    subplot(length(kprev),length(kprev),m)
        smoothhist2D(tsample(:,[mod(m-1,length(kprev))+1, floor((m-1)/length(kprev))+1]),1,[50 50])
        xlabel(param{mod(m-1,length(kprev))+1})
        ylabel(param{floor((m-1)/length(kprev))+1})
end

%%
% Scatter plot of the joint probability of two parameters
%

figure(8)
for m = 1:length(kprev)^2
    subplot(length(kprev),length(kprev),m)
    if((mod(m-1,length(kprev))+1)==floor((m-1)/length(kprev))+1)
        [minquantile,maxquantile] = mininterval(tsample(:,(mod(m-1,length(kprev))+1)),0.95);
        kde(tsample(:,(mod(m-1,length(kprev))+1)))
        gridxy(minquantile)
        gridxy(maxquantile)
        xlabel(param{floor((m-1)/length(kprev))+1})
    else
        scatter(tsample(:,mod(m-1,length(kprev))+1),tsample(:,floor((m-1)/length(kprev))+1),4,prob(1:thinstep:end))
        xlabel(param{mod(m-1,length(kprev))+1})
        ylabel(param{floor((m-1)/length(kprev))+1})
    end
end

%%
% iterative time-course of the parameters
%

figure(9)
for m = 1:length(kprev)
    subplot(3,3,m)
    plot(sample(:,m))
    xlabel('iteration')
    ylabel(param{m})
end


%%
%

figure(11)
clf

plot(prob,'.-r')
%hold on
%
xlabel('sample')
legend('log a-posteriori')

figure(12)
clf 
plot(lh,'.-g')
xlabel('sample')
legend('acceptance rate [%]')





%% Credibility Intervals
% 

%%
% 95% 
%
    display(['----------------------------------------------'])
for m = 1:length(kprev)
    [minquantile,maxquantile] = mininterval(tsample(:,m),0.95);
    display(['95% Credibility Interval for ' param{m} ': [' num2str(minquantile) ' , ' num2str(maxquantile) ']']);
end

%%
% 75%
%
    display(['----------------------------------------------'])
for m = 1:length(kprev)
    [minquantile,maxquantile] = mininterval(tsample(:,m),0.75);
    display(['75% Credibility Interval for ' param{m} ': [' num2str(minquantile) ' , ' num2str(maxquantile) ']']);
end

%%
% 50%
%
    display(['----------------------------------------------'])
for m = 1:length(kprev)
    [minquantile,maxquantile] = mininterval(tsample(:,m),0.50);
    display(['50% Credibility Interval for ' param{m} ': [' num2str(minquantile) ' , ' num2str(maxquantile) ']']);
end



%% Covariance Matrix
%

cov(sample)

display(['Acceptance Rate: ' num2str(acc/(nacc+acc)*100) '%'])

clear vxi
clear vxu
tn =1000;
vxi = zeros(size(tsample,1),tn);
vxu = zeros(size(tsample,1),tn);
tt=linspace(t(1),t(end),tn);
for kk = 1:size(tsample,1)
    x0(1)       = tsample(kk,3)/tsample(kk,4);
    x0(2) 		= tsample(kk,5)/tsample(kk,6)*x0(1);
    x0(3)       = tsample(kk,3)/tsample(kk,4);
    x0(4)       = tsample(kk,5)/tsample(kk,6)*x0(3);
    [tx,xx]=odeGAL1p( tt, x0, tsample(kk,1:6), input, [reltols abstols 10] );
    vxi(kk,:) = tsample(kk,1)*xx(2,:);
    vxu(kk,:) = tsample(kk,1)*xx(4,:);
end
ttt=repmat(tt,size(tsample,1),1);
[bandwidth,density,X,Y]=kde2d([vxi(:),ttt(:)]);
surf(X,Y,density,'LineStyle','none'), view([0,70])
colormap hot, hold on, alpha(.8)
zlim([0, 10e-3])


save mcmc_metropolis.mat
