# How to Implement and Minimize Maximum Likelihood Expression in MATLAB

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));

• Why do you need to iterate? The equation doesn't require it. – Peter K. Aug 27 '16 at 23:15
• @PeterK.:I think in my estimation technique I have to solve using Expectation Maximization which is an iterative optimization algorithm. In my other estimation problem using another model, I have to apply Gradient descent. So, you see that this example would help me to later adapt to other estimation techniques. – SKM Aug 28 '16 at 0:31
• The numerical solution of the ML estimate is quite different from writing the equation. You seem to be confusing the two. – Peter K. Aug 28 '16 at 0:40
• @PeterK.: I have updated the Question based on your suggestions. Please let me know if it is clear or not. Thank you. – SKM Aug 29 '16 at 16:55
• That's somewhat better... but do you know $\mathbf{z}_n$ ? If not, do you know something about it? – Peter K. Aug 29 '16 at 19:37

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)


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?
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).
• @SKM Yes, if your H_est values need to be sorted to get a unimodal plot. – Peter K. Aug 29 '16 at 21:44