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']);