Contents
clc; clear; close all;
set(0, 'defaultTextInterpreter', 'latex');
est_fn_ray = @(samples) sqrt(mean(samples.^2)/2);
Q1
Ns = logspace(1, 5, 5);
params = 1:4;
figure('Position', [0 0 750 1000]);
t1 = tiledlayout(length(params), 3);
sgtitle("Exponential R.V.");
figure('Position', [0 0 750 1000]);
t2 = tiledlayout(length(params), 3);
sgtitle("Rayleigh R.V.");
for param=params
exp_res = zeros(length(Ns), 3);
ray_res = zeros(length(Ns), 3);
for Ni=1:length(Ns)
for M=100
N = Ns(Ni);
[exp_bias, exp_var, exp_mse] = ...
run_experiment(N, M, param, @exprnd, @mean);
[ray_bias, ray_var, ray_mse] = ...
run_experiment(N, M, param, @raylrnd, est_fn_ray);
exp_res(Ni, :) = [exp_bias exp_var exp_mse];
ray_res(Ni, :) = [ray_bias ray_var ray_mse];
fprintf("N=%d param=%d: %f %f %f; %f %f %f\n", ...
N, param, ...
exp_bias, exp_var, exp_mse, ...
ray_bias, ray_var, ray_mse);
end
end
nexttile(t1);
semilogx(Ns, exp_res(:,1));
title(sprintf('Bias ($$\\mu=%d$$)', param));
xlabel("Samples");
ylabel("$$E[\hat\mu-\mu]$$");
grid on;
nexttile(t1);
loglog(Ns, exp_res(:,2));
title('Variance');
xlabel("Samples");
ylabel("$$E[(\hat\mu-\bar{\hat\mu})^2]$$");
grid on;
nexttile(t1);
loglog(Ns, exp_res(:,3));
title('MSE');
xlabel("Samples");
ylabel("$$E[(\hat\mu-\mu)^2]$$");
grid on;
nexttile(t2);
semilogx(Ns, ray_res(:,1));
title(sprintf('Bias ($$\\alpha=%d$$)', param));
xlabel("Samples");
ylabel("$$E[\hat\alpha-\alpha]$$");
grid on;
nexttile(t2);
loglog(Ns, ray_res(:,2));
title('Variance');
xlabel("Samples");
ylabel("$$E[(\hat\alpha-\bar{\hat\alpha})^2]$$");
grid on;
nexttile(t2);
loglog(Ns, ray_res(:,3));
title('MSE');
xlabel("Samples");
ylabel("$$E[(\hat\alpha-\alpha)^2]$$");
grid on;
end
Q2
load('data');
data = data.';
x = linspace(min(data), max(data), 1000);
mu_est = mean(data);
alpha_est = est_fn_ray(data);
lambda_est = 1/mu_est;
exp_pdf = @(x) lambda_est * exp(-lambda_est * x);
alpha2_est = alpha_est^2;
ray_pdf = @(x) x/alpha2_est .* exp(-x.^2/(2*alpha2_est));
exp_likelihood = sum(log(exp_pdf(data)));
ray_likelihood = sum(log(ray_pdf(data)));
fprintf("Log-likelihood of exponential distribution: %f\n" ...
+ "Log-likelihood of rayleigh distribution: %f\n", ...
exp_likelihood, ray_likelihood);
Log-likelihood of exponential distribution: 1053.462453
Log-likelihood of rayleigh distribution: 1365.516070
show histogram to visually see which distribution fits better
figure();
hold on;
histogram(data, 'Normalization', 'pdf');
plot(x, exp_pdf(x));
plot(x, ray_pdf(x));
legend(["Sample data (normalized to PDF)", ...
sprintf("Exponential PDF (mu=%f)", mu_est), ...
sprintf("Rayleigh PDF (alpha=%f)", alpha_est)]);
title("Data vs. Exponential and Rayleigh ML-Estimated Distributions");
ylabel("PDF");
xlabel("Values");
helper function to run experiment
function [bias, variance, mse] = run_experiment(N, M, param, randfn, estfn)
samples = randfn(param, N, M);
est = estfn(samples);
bias = mean(est) - param;
variance = var(est);
mse = variance + bias^2;
end
N=10 param=1: 0.002574 0.078209 0.078215; -0.009523 0.028445 0.028536
N=100 param=1: 0.007988 0.011183 0.011247; -0.005721 0.002533 0.002566
N=1000 param=1: 0.000378 0.001286 0.001286; -0.001020 0.000293 0.000294
N=10000 param=1: -0.000920 0.000098 0.000099; -0.000426 0.000024 0.000025
N=100000 param=1: 0.000247 0.000008 0.000008; 0.000111 0.000002 0.000002
N=10 param=2: -0.040950 0.283109 0.284785; 0.011830 0.122604 0.122744
N=100 param=2: -0.007575 0.048410 0.048468; -0.024112 0.009429 0.010010
N=1000 param=2: -0.005067 0.004940 0.004965; 0.003536 0.000870 0.000883
N=10000 param=2: -0.000388 0.000466 0.000466; 0.000211 0.000100 0.000100
N=100000 param=2: 0.000116 0.000036 0.000036; -0.000221 0.000014 0.000014
N=10 param=3: -0.043354 0.686853 0.688733; -0.028280 0.211181 0.211981
N=100 param=3: -0.104062 0.082007 0.092836; 0.005258 0.020529 0.020557
N=1000 param=3: -0.013672 0.009392 0.009579; -0.005131 0.002462 0.002488
N=10000 param=3: -0.005534 0.000959 0.000990; 0.000608 0.000229 0.000229
N=100000 param=3: -0.000150 0.000119 0.000119; -0.000005 0.000019 0.000019
N=10 param=4: 0.127984 1.551503 1.567883; -0.033811 0.389189 0.390333
N=100 param=4: 0.019660 0.182156 0.182542; -0.060038 0.042804 0.046409
N=1000 param=4: 0.011553 0.015798 0.015931; -0.005465 0.003734 0.003764
N=10000 param=4: 0.006851 0.001716 0.001763; 0.001879 0.000401 0.000405
N=100000 param=4: 0.001641 0.000162 0.000164; -0.000130 0.000040 0.000040