%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code estimates the mean and the variance from a sample
% of i.i.d. normal random variables using maximum likelihood
% estimation.
Copyright By PowCoder代写 加微信 powcoder
% You can modify:
% T: the sample size
% mu0: the mean of the population
% sigma0sq: the variance of the population
% The standard errors of the MLE estimates are computed in 3 different ways:
% 1) We compute the first derivative by hand
% 2) We compute the second derivative by hand
% 3) We compute the first derivatives using the gradient function in matlab (in the symbolic math
% toolbox)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear environment
clear all; clear variables; clc;
% fix seed of random number generator
% number of observations
T = 1000; %sample size
% true parameters in the population
sigma0 = 2;
sigma0sq = sigma0^2;
% sample is iid from normal(mu0,sigma0^2)
X = normrnd(mu0,sigma0,T,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ESTIMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ESTIMATION METHOD 1: Naive estimates
%Estimate with sample mean and sample variance
meanX = mean(X);
stdX = std(X);
varX = var(X);
%ESTIMATION METHOD 2: Maximum Likelihood estimation
% initial guess
parameter_guess = [1,1];
% maximization of the likelihood
fprintf(‘Maximum Likelihood Estimation\n’)
options = optimset(‘Display’,’off’,’TolFun’,1e-5,’TolX’,1e-5);
% compare results between methods
opt_actual_table = table([mu0 sigma0sq]’, [meanX varX]’, para_opt’,’VariableNames’,{‘True’ ‘Naive’ ‘MLE’},…
‘RowNames’,{‘Mean’ ‘Variance’});
disp(‘Actual vs. Naive vs. MLE estimates’);
disp(opt_actual_table);
disp(‘Note: The naive estimates are the sample mean and the sample variance (the sample variance is divided by T-1 for unbiasedness)’)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% INFERENCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(‘——————————————————————————-\n’)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Method 1: Compute first derivative (gradient) by hand
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(‘Compute standard errors: Method 1 – first derivative (gradient) by hand\n’)
muhat = para_opt(1); % estimate for mean
sigmasqhat=para_opt(2); % estimate for variance
varcov_mle = zeros(2,2); % initialize matrix
% compute gradient at time t as a vector
gra_t = [(X(t)-muhat)./(sigmasqhat) , -(1/(2*sigmasqhat)) + (1/(2*(sigmasqhat^2)))*((X(t)-muhat)^2) ]’;
% add to the var/cov matrix
varcov_mle = varcov_mle + gra_t * gra_t’;
% invert the matrix and divide by T to estimate varcov of parameters
varcov = (1/T)*inv((1/T)*varcov_mle);
var_mle = diag(varcov); % variances are on the diagonal
stderror_mle = sqrt(var_mle); % std errors are square roots of the variances
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Method 2: Compute second derivative (Hessian) by hand
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(‘Compute standard errors: Method 2 – second derivative (Hessian) by hand\n’)
varcov_mle2 = zeros(2,2); %initialize matrix
% compute Hessian at time t as a matrix
hes_t = [-1/(sigmasqhat) , -(X(t)-muhat)./(sigmasqhat^2) ; …
-(X(t)-muhat)./(sigmasqhat^2), (1/(2*sigmasqhat^2)) – (1/(sigmasqhat^3))*((X(t)-muhat)^2)];
% add to the var/cov matrix
varcov_mle2 = varcov_mle2 + hes_t;
varcov2 = -(1/T)*inv((1/T)*varcov_mle2);
var_mle2 = diag(varcov2); % variances are on the diagonal
stderror_mle2 = sqrt(var_mle2); % std errors are square roots of the variances
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Method 3: Compute gradient using the symbolic tools in Matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(‘Compute standard errors: Method 3 – gradient using matlab symbolic functions\n’)
% set up symbolic values for the parameters
syms mu sigmasq
% vector of symbolic likelihood
likelihood = zeros(T,1,’sym’);
likelihood = -(1./2).*log(sigmasq) – (1./(2.*sigmasq))*((X-mu).^2);
% gradient initialized with zeros
gra_lik = zeros(T,2,’sym’) ;
% matrix of var-cov initialized with zeros
varcov_g = zeros(2,2,’sym’);
% take derivative of the likelihood for each observation
gra_lik(t,:) = gradient(likelihood(t,1), [mu sigmasq]);
% compute var/cov for that observation and add to matrix var/cov
varcov_g = varcov_g + gra_lik(t,:)’*gra_lik(t,:);
% substitute the value of the parameters estimated with MLE
varcov_g_mle = subs(varcov_g, [mu sigmasq], [para_opt(1), para_opt(2)]);
% var/cov should now become numerical
varcov_g_mle = double(varcov_g_mle);
% compute var/cov of the parameters
varcov3 = (1/T)*inv((1/T)*varcov_g_mle);
var_mle3 = diag(varcov3); % the variances are on the diagonal
stderror_mle3 = sqrt(var_mle3); % the std errors are square roots of variances
% display all results
% compare results between methods
estimate_table = table([mu0 sigma0sq]’, [meanX varX]’, para_opt’, stderror_mle, stderror_mle2, stderror_mle3, ‘VariableNames’,{‘True’ ‘Naive’ ‘MLE’ ‘Std_Error_1’ ‘Std_Error_2’ ‘Std_Error_3’},…
‘RowNames’,{‘Mean’ ‘Variance’});
disp(‘——————————————————————————‘);
disp(estimate_table);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Log-likelihood function for iid normal (mu,sigma^2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = loglik(parameters, sample)
mu = parameters(1); % mean
sigmasq = parameters(2); % variance
T = length(sample); % sample size
% compute likelihood of iid normal sample
f = -(T/2).*log(sigmasq) – (1/(2*sigmasq))*sum((sample – mu).^2);
% use likelihood/T and make sure to change sign to feed into fminsearch
% we need to change sign because fminsearch looks for a minimum, not a max
f = -f/T; % for fminsearch
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com