
paretopdf = @(x, a, b) (x > b) .* (a * b^a * x.^(-a-1));
    
alphas = [0.1, 2.0, 1.0];
betas = [0.1, 0.1, 1.0];

X = 0:0.01:4;
    
for i = 1:3

    alpha = alphas(i);
    beta = betas(i);
    
    alpha_post = alpha + 3;
    beta_post = max([0.7, 1.3, 1.9, beta]);
    
    figure()   
    subplot(2, 1, 1)
    plot(X, paretopdf(X, alpha, beta))
    title(sprintf('Prior: alpha = %f, beta = %f', alpha, beta))
    
    subplot(2, 1, 2)
    plot(X, paretopdf(X, alpha_post, beta_post))
    title(sprintf('Posterior: alpha = %f, beta = %f', alpha_post, beta_post))
    
    saveas(gcf, sprintf('q4_%d.eps', i), 'eps')
    
end