load('fig4_histogram_data.mat')
R = data.R;
llik = data.llik;
kappa = data.kappa;

%%%%%%%%%%%%%%%%%

x_vals = 0:0.01:4;  
y_vals = 2 * x_vals .* dlrlim(x_vals.^2);
dolphinBlue = 0.3 * [0 0.4470 0.7410] + 0.7 * [0.65 0.65 0.65];
figure()
histogram(sqrt(2*llik./kappa), 50, 'Normalization', 'pdf', 'FaceColor', dolphinBlue, 'EdgeColor', 'black');
hold on;
plot(x_vals, y_vals, 'LineWidth', 2, 'Color', 'blue');
hold on
fplot(@(x) 2 * normpdf(x), [0, 3.1], 'LineStyle', '--', 'LineWidth', 2, 'Color', 'black');

ylim([0,1.05])

xlabel('$(2\ell(\hat{\theta}_n)/\hat{\kappa}_n)^{1/2}$', 'Interpreter', 'latex');
ylabel('Density', 'Interpreter', 'latex');

legend({'Empirical', 'Limit density', 'Half-normal'}, ...
       'Interpreter', 'latex', 'Location', 'northeast');

box on;
set(gca, 'FontSize', 16);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x_vals2 = 0:0.01:1;
y_vals2 = ones(size(x_vals2));
figure()
histogram(R, 50, 'Normalization', 'pdf', 'FaceColor', dolphinBlue, 'EdgeColor', 'black');
hold on;
plot(x_vals2, y_vals2, 'LineWidth', 1, 'Color', 'blue', 'LineWidth', 2);

xlabel('$R$', 'Interpreter', 'latex');
ylabel('Density');

legend({'Empirical', 'Uniform limit'}, 'Location', 'northeast');

box on;
set(gca, 'FontSize', 16);

