How to Implement and Minimize Maximum Likelihood Expression in MATLAB

Question

I need to check if the estimation algorithm has converged or not. I am using the Maximum Likelihood estimation method. For convergence check, we see if the log-likelihood has reached its maximum value or not. But I am having difficulty in implementing the log-likelihood expression.

As an example, I am estimating the model parameters of a Moving Average model of order d =3 expressed in Eq(1). The known coefficients are h = [1 0.45 -0.2].The pdf is given in Eq(2) and the log-likelihood in Eq(3).

I have tried the following steps. For 9 different values of the channel coefficients, I calculate the likelihood expression. I then iterate 9 times to get 9 different values. After the vector of likelihood values are obtained, I should take the log of the values of likelihood and plot. However, I am stuck in the plot as I cannot understand how to show the maximum for the coefficients of the channel which is in a vector.

Question 1: I am stuck in implementation in the step. How do I plot the curve and show at which value of $H$ the curve reaches its maximum value?

Question 2: I want to add noise of 5dB using the awgn() function as follows

 y=awgn(y,5,'measured');

But in that case I cannot understand what value I should take for the standard deviation of the noise, s?

 clear all


N=256; %Number of Samples to collect

H_est = [1.07648398119767,0.469220009437168,-0.0245459792881367;
1.01925235521823,0.449715366555421,-0.0287160909270392;
1.05903405523260,0.568022131010140,0.0261264158909992;
1.03693902020750,0.421708336037157,-0.102062737844706;
1.07969493608754,0.481505366768731,-0.0505480628598221;
1.05948950849694,0.450984209695615,-0.110305164694739;
0.933641686721610,0.310664477120195,-0.0803322291311877;
1.13980384345926,0.530800706049278,-0.0314494051827266;
0.949926376465499,0.498110774466918,0.0473466083388137];
 %Generating source information signal 
z = rand(1,N);
h = [1    0.45   -0.2]; %true known channel impulse response

y = filter(h,1,z);
y=y+randn(1,N);
s=1; %Assume standard deviation s=1
tol = 1e-4;
maxiter = 9; %because there are 9 values of the channel coefficients
llh = -inf(1,maxiter);

 llh=zeros(1,9); %Place holder for likelihoods

for iter=1:9
%Calculate Likelihoods for each parameter value in the range


 H = H_est(iter,:);


     Hz=  filter(H,1,z);
        llh(iter) = exp(-sum((y-Hz).^2)/(2*s^2));  %Neglect the constant term (1/(sqrt(2*pi)*sigma))^N as it will pull %down the likelihood value to zero for 

      %  increasing value of N
        if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter-1)); break; end   % check likelihood for convergence
  end
        [maxL,index]=max(llh); %Select the parameter value with Maximum Likelihood
        display('Maximum Likelihood of H');
display(H_est(index,:));

plot(1:9,log(llh));

Answer 1

OK, let's see if I can answer this now.

The original log likelihood expression is $$ L(\mathbf{y} | \mathbf{\Theta}) = -\frac{N}{2} \ln (2\pi \sigma^2_w) - \frac{1}{2\sigma^2_w} \sum_{n=0}^{N-1} ( y_n - \mathbf{h}^T \mathbf{z}_n)^2 $$ where $y_n, n=0,1,\ldots,N-1$ are the known noisy measurements, $\mathbf{h}$ are your (unknown) filter coefficients and $\mathbf{z}_n$ are the known random inputs to your filter.Here $\mathbf{\Theta}$ is just your unknown $\mathbf{h}$ coefficients if that is all that is to be estimated.

So, the way to calculate it is just (in R) as:

N <- 256
h_true <- c(1, 0.45, -0.2)    
z <- rnorm(N,0,1)
hz <- filter(z,h_true)
sigma_w <- 1

y <- hz + rnorm(N,0,sigma_w)

llh <- -N/2*log(2*pi*sigma_w*sigma_w) - 1/(2*sigma_w*sigma_w)*sum((y-hz)^2)

Now, to your question:

Question 1: I am stuck in implementation in the step. How do I plot the curve and show at which value of $H$ the curve reaches its maximum value?

Well, that's a little hard because your log likelihood is a multi-dimensional plot (four dimensional in your specific case -- one for each coefficient and one for the value of the log likelihood).

So you need to plot llh on three dimensions H_est(:,1), H_est(:,2), and H_est(:,3) and then the value. You might want to just look at plot3(H_est(:,1), H_est(:,2), llh) (or something that approximate this).