0
$\begingroup$

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

image

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));
$\endgroup$
  • $\begingroup$ Why do you need to iterate? The equation doesn't require it. $\endgroup$ – Peter K. Aug 27 '16 at 23:15
  • $\begingroup$ @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. $\endgroup$ – SKM Aug 28 '16 at 0:31
  • $\begingroup$ The numerical solution of the ML estimate is quite different from writing the equation. You seem to be confusing the two. $\endgroup$ – Peter K. Aug 28 '16 at 0:40
  • $\begingroup$ @PeterK.: I have updated the Question based on your suggestions. Please let me know if it is clear or not. Thank you. $\endgroup$ – SKM Aug 29 '16 at 16:55
  • $\begingroup$ That's somewhat better... but do you know $\mathbf{z}_n$ ? If not, do you know something about it? $\endgroup$ – Peter K. Aug 29 '16 at 19:37
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Thank you for your reply. Do I need to sort the values of H_est() in order to get a unimodal curve that has one peak? $\endgroup$ – SKM Aug 29 '16 at 21:06
  • $\begingroup$ Hi Peter, why did you start using R. I'd expect you to work on MATLAB / Python or maybe Julia. R is really a curiosity generating choice :-). $\endgroup$ – Royi Aug 29 '16 at 21:40
  • $\begingroup$ @SKM Yes, if your H_est values need to be sorted to get a unimodal plot. $\endgroup$ – Peter K. Aug 29 '16 at 21:44
  • $\begingroup$ @Drazick I needed to learn R for work, and decided I'd learn it quicker if I used it to answer things here. I was right! I found trying to write algorithms I know already in matlab / scilab in R stretches the boundaries much more than just using it in my day job. $\endgroup$ – Peter K. Aug 29 '16 at 21:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.