function w8c_laplace_question_code()
% MLPR, University of Edinburgh
% Iain Murray, 2014, 2016
%
% The script runs in Octave, but the plots only come out exactly as intended in Matlab.
% (And Matlab seems buggy, I sometimes have to run the script twice for the
% plots to export properly.)
%
% I'm not porting this script to Python, which seems pointless grunt work.
% There's no interesting computation here, and I've provided other plotting demos.

set(0,'defaulttextinterpreter','latex');

normal = @(x,m,v) exp(-0.5*(x-m).^2/v)/sqrt(2*pi*v);
hh = 0.01; % discretization for plots

for rr = [2 20]
    % Posterior over lambda
    figure(1); clf; hold on;
    lmu = rr-1;
    lstd = sqrt(rr-1);
    lambda = (lmu-4*lstd):hh:(lmu+4*lstd); % Points to plot at
    pstar_lambda = exp(-lambda).*(lambda.^(rr-1)).*(lambda>0);
    % numerically normalizing pstar because lazy. Not accurate if range of plot
    % doesn't actually capture all the mass (when approximation is poor), so one
    % of the plots in the answers document is slightly wrong:
    p_lambda = pstar_lambda / (sum(pstar_lambda)*hh);
    laplace_lambda = normal(lambda, lmu, lstd^2);
    plot(lambda, p_lambda, 'b');
    plot(lambda, laplace_lambda, 'r--');
    xlabel('$\lambda$');
    ylabel('$P(\lambda\,|\,r)$');
    xtix = get(gca, 'xtick');
    q3_saveplot(sprintf('q3lambda_r%d', rr));
    clear lambda lmu lstd
    
    % Posterior over log(lambda)
    figure(2); clf; hold on;
    lmu = log(rr);
    lstd = sqrt(1/rr);
    log_lambda = (lmu-4*lstd):hh:(lmu+4*lstd);
    explog_lambda = exp(log_lambda);
    pstar_loglambda = exp(-explog_lambda).*(explog_lambda.^rr);
    p_loglambda = pstar_loglambda / (sum(pstar_loglambda)*hh);
    laplace_loglambda = normal(log_lambda, lmu, lstd^2);
    plot(explog_lambda, p_loglambda, 'b');
    plot(explog_lambda, laplace_loglambda, 'r--');
    set(gca, 'xscale', 'log');
    % I've done the slightly weird thing of plotting lambda on the x-axis and
    % then using a log scale. So the points are spaced just as they would be if
    % I'd plotted log-lambda. This way the numbers are more easily comparable to
    % the other plot.
    xlabel('$\lambda$');
    ylabel('$P(\log\lambda\,|\,r)$');
    set(gca, 'xtick', xtix);
    q3_saveplot(sprintf('q3loglambda_r%d', rr));
end


function q3_saveplot(figname)

% Use the whole plot area:
axis tight
% Once again, Matlab is annoying. The mode gets clipped off a little when
% plotting with thick lines. Tweak y-axis height:
ax = axis;
ax(4) = 1.01*ax(4);
axis(ax);

% Make line widths and fonts reasonable for inclusion as a small figure:
set(gcf, 'PaperPosition', [0 0 4 2]*2);

% Matlab bug? The axis labels get clipped off for me with the above
% PaperPosition. I fix that by nudging the plot position and size:
pos = get(gca, 'position'); % [left bottom width height]
pos(1:2) = pos(1:2) + 0.2;
pos(3:4) = pos(3:4) - 0.2;
set(gca, 'position', pos);

% For me, -dpdf doesn't work well, so I export to .eps and convert:
print([figname, '.eps'], '-depsc');
system(['epstopdf ' figname, '.eps']);