# Can someone show the details of how to apply AIC for sinusoidal models to specific data?

NOTE: This is a question that another user has been trying (unsuccessfully) to ask. Because the multiple questions asking, essentially, the same thing have either been deleted by me (because they were duplicates of each other) or by the user, I have tried to ask the question in such a way that will help the user and is a good fit for the site.

Suppose I have a vector of data measurements, $x(n)$, and I want to fit a model to that data: $$x(n) = \sum_{p=1}^{P} a_p \cos(2\pi f_p + \phi_p) + \epsilon(n)$$ where each of the $P$ sinusoids is parameterized by $(a_p, f_p, \phi_p)$ and $\epsilon$ is an i.i.d. Gaussian, white noise process with zero mean and variance $\sigma^2$.

My question is: How do we use AIC (or BIC or MDL) to estimate the model order, $P$ ?

The Wikipedia entry just says: $$AIC = 2P - 2\ln(L)$$ where $L$ is the likelihood function.

How do I calculate $L$ given $x(n)$?

• So the question is not how to choose $P$, but rather how to choose it using AIC (or BIC or MDL) (?) – Matt L. Mar 30 '14 at 19:20
• yes exactly,how to choose P using criteria – dato datuashvili Mar 30 '14 at 19:23
• In the previous question it was about matlab code. The problem is aic or others does not provide order information. It is aic test command in matlab but using that gives 0 order so something was wrong with the data or the method. – s.s.o Mar 30 '14 at 19:23
• this data was given from deterministic components,$cos(2*pi*f*t)$ linear combination of such models +white noise,can we apply AR model? – dato datuashvili Mar 30 '14 at 19:25
• in matlab site,there is declared matrix with size $4X4$,but in reality it could be large right? – dato datuashvili Mar 30 '14 at 19:30

The PDF (likelihood function) of the multi-tone estimation problem is: $${\cal L}({\bf x}; {\bf a}, {\bf f},{\bf \phi}, {\sigma}^2) = \frac{1}{(2\pi \sigma^2)^{\frac{N}{2}}} \exp\left( - \frac{1}{2\sigma^2} \sum_{n=0}^{N-1} ( x(n) - \sum_{p=1}^{P} {a}_p \cos(2\pi {f}_p + {\phi}_p) )^2 \right)$$

The paper by Djuric quotes the AIC and MDL values as:

The code at the end attempts to:

• Generate a 5-tone noisy signal.
• Do the ML estimates of the amplitudes, frequencies and phases for a given value of $p$.
• Find the (log) likelihood values for each value of $p$.
• Calculate the AIC and MDL values from Djuric.

The plot for $p = 1 \ldots 60$ is:

Usually, the AIC over-estimates the model order. In this particular case, AIC gets 7 and MDL gets 5 (the correct value). The minima of each curve are plotted (as a o on the AIC curve and as a + on the MDL curve).

Djuric, P.M., "A model selection rule for sinusoids in white Gaussian noise," Signal Processing, IEEE Transactions on (Volume:44 , Issue: 7 ), pp. 1744 - 1751.

scilab CODE ONLY BELOW

P = 5;
f = [ 0.1 0.12 0.231 0.333 0.478634];
phi = [0.329847 5.90324 3.09834 4.907983 1.984];
a = [ 1 1 2 2 3];
sigma = 1;

N = 100;

x = cos(2*%pi*[0:N-1]'*f + ones(N,1)* phi )*a';

xn = x + sigma*rand(N,1,'normal');

function [L,LL] = likelihood(data,ahat,fhat,phihat,sigmahat)
N = length(data);
datahat = cos(2*%pi*[0:N-1]'*fhat + ones(N,1)* phihat )*ahat';
L = 1/(2*%pi*sigmahat*sigmahat)^(N/2) * exp(-1/(2*sigmahat*sigmahat)*sum((data - datahat).^2));
LL = - log(2*%pi*sigmahat*sigmahat)*N/2 - -1/(2*sigmahat*sigmahat)*sum((data - datahat).^2);
endfunction

function [a,m] = aic(p,L,N)
a = 3*p - log(L);
m = 3*p/2*log(N) - log(L);
endfunction

function [ahat,fhat,phihat,regr] = sin_estimate(data)
N = length(data);
ahat = mx/(N/2);
phihat = atan(imag(DATA(ix)),real(DATA(ix)));

regr = data - cos(2*%pi*[0:N-1]'*fhat + ones(N,1)* phihat )*ahat
endfunction

pValues = [1:60];
likeValues = [];
loglikeValues = [];
aicValues = [];
mdlValues = [];
for p=pValues,
signal = xn;
a_est = zeros(1,p);
f_est = zeros(1,p);
phi_est = zeros(1,p);
for pp=1:p
[a_est(pp),f_est(pp),phi_est(pp),signal] =  sin_estimate(signal);
end

[l,ll] = likelihood(xn,a_est,f_est,phi_est,sigma);

likeValues = [likeValues l];
loglikeValues = [loglikeValues ll];

[aicv,mdlv] = aic(p,likeValues(p));
aicValues = [aicValues; aicv]
mdlValues = [mdlValues; mdlv]
end

clf;
plot(aicValues);
plot(mdlValues,'r');
title('AIC and MDL plots');

• can i apply this code in matlab? – dato datuashvili Mar 31 '14 at 4:35
• but how we have introduced initial values?i can't understand code – dato datuashvili Mar 31 '14 at 4:52
• i have tried following code stackoverflow.com/questions/22754460/… and what is reason of error? – dato datuashvili Mar 31 '14 at 6:20
• 1) Yes, the code should work (with some minor changes: e.g. %pi $\rightarrow$ pi). 2) The initial values are estimated using a zero-padded FFT. 3) See comment there. – Peter K. Mar 31 '14 at 11:42
• what about my code for AIC/BIC ? – dato datuashvili Mar 31 '14 at 11:45