# Math Help - Heavy duty MLE

1. ## MLE challenge

I will be impressed if someone can get this. I will also be impressed if someone can explain why this problem is of great interest.

Two correlated random variables $x_i$ and $y_i$ are Bivariate normal

$\left[ \begin{array}{ccc} x_i \\ y_i\\ \end{array} \right]$ ~ $N(\left[ \begin{array}{ccc} 0 \\ 0\\ \end{array} \right], VC = \left[ \begin{array}{ccc} 1 & a\\ a & 1\\ \end{array} \right])$ Where $VC$ is the Variance Covariance matrix.

Given $n$ $i.i.d.$ samples the likelihood function is:

$lik= (2\pi)^{-n}|VC|^{-n/2} exp(-1/2\sum^n [x_i, y_i]VC ^{-1} \left[ \begin{array}{ccc} x_i \\ y_i\\ \end{array} \right]).$ (Assume $1> a^2$)

Some useful information:

$|VC|= (1-a^2),$ $VC^{-1}= \frac{1}{1-a^2}\left[ \begin{array}{ccc} 1 & -a\\ -a & 1\\ \end{array} \right]$

Find the MLE and asymptotic variance.

2. I'm working out the answer now. If any one is curious I will post my solution.

3. I assume you want the MLE of a.

4. Yes I do. It's tedious to calculate, but can you see why this problem is interesting?

5. Well, there are formulas to estimate a covariance matrix... so it should be possible to narrow it to a huh ?

6. Did I hear someone thinking maximization of a cubic?

7. Originally Posted by Anonymous1
Did I hear someone thinking maximization of a cubic?
Look, if you know the answer then say so and/or post it. If this was meant to be a challenge question to members then it should have been posted in the subforum whose link is given below and the rules of that subforum (plainly stated in the stickies) followed:

http://www.mathhelpforum.com/math-he...enge-problems/

8. Been working on this for a while now so I thought I'd share .

Here is a data correlation project that involves finding the solution to the above $MLE$ estimator, which entailed solving a cubic.

My instructor gave us this problem to point out that $MLE$ can run into trouble for "multiple hump" functions, $e.g.,$ cubics. This is because there may be various local maximums.

To make sure I found the proper solution to: $\frac{dl(\hat a)}{d} = 0$ $i.e.,$ the $global$ maximum, I used the naive estimator, $\hat a_0,$ as my initial guess $x_0 ,$ along with 10-step Newton's.

I hope someone finds this interesting. Cause I do!

Code:
u1= load('UNITED.txt');

m1= sum(abs(u1.^2));    u1= u1./sqrt(m1);
m2= sum(abs(u2.^2));    u2= u2./sqrt(m2);
m3= sum(abs(u3.^2));    u3= u3./sqrt(m3);
m4= sum(abs(u4.^2));    u4= u4./sqrt(m4);

k= [20, 50, 100, 200, 300, 400];
a= sum(u1.*u2);

aHatUS0= zeros(100,length(k));  aHatUS= zeros(100,length(k));
aHatHK0= zeros(100,length(k));  aHatHK= zeros(100,length(k));
for t= 1:100
for i= 1:length(k)
R= randn(length(u1),k(i));
v1= R'*u1;  v2= R'*u2; v3= R'*u3;  v4= R'*u4;
aHatUS0(t,i)= (1/k(i)).*sum(v1.*v2);
aHatHK0(t,i)= (1/k(i)).*sum(v3.*v4);

xn1= aHatUS0;   xn2= aHatHK0;   n= length(v1);
for j=1:10
a1= xn1;    a2= xn2;

f1= n*a1.^3 - a1.^2 .* sum (v1.*v2) + a1.*(-n+sum(v1.^2+v2.^2)) - sum (v1.*v2);
fp1= 3*n*a1.^2 - 2.*a1.* sum (v1.*v2) -n + sum(v1.^2+v2.^2);
f2= n*a2.^3 - a2.^2 .* sum (v3.*v4) + a2.*(-n+sum(v3.^2+v4.^2)) - sum (v3.*v4);
fp2= 3*n*a2.^2 - 2.*a2.* sum (v3.*v4) -n + sum(v3.^2+v4.^2);

xn1= a1 - f1./fp1;  xn2= a2 - f2./fp2;
end
aHatUS(t,i)=xn1(1);
aHatHK(t,i)=xn2(1);
end
end

aHatUS0Var= (m1*m2+a^2)*k.^(-1);   aHatUSVar= (1-a^2)^2./((1+a^2)*k);
aHatHK0Var= (m3*m4+a^2)*k.^(-1);   aHatHKVar= (1-a^2)^2./((1+a^2)*k);

mSEaHatUS0= aHatUS0Var + (mean(aHatUS0-a)).^2;
mSEaHatHK0= aHatHK0Var + (mean(aHatHK0-a)).^2;
mSEaHatUS= aHatUSVar + (mean(aHatUS-a)).^2;
mSEaHatHK= aHatHKVar + (mean(aHatHK-a)).^2;

figure; hold on;    grid on;
plot(k,mSEaHatUS0,'g');
plot(k,aHatUS0Var,'b');
plot(k,mSEaHatHK0,'r');
plot(k,aHatHK0Var,'y');

title('Figure 1: Empirical MSEs and Theoretical Variances')
xlabel('k'); ylabel('MSE & Variance');
h=legend('Empirical MSE of aUS0_{hat}','Theoretical Variance of aUS0_{hat}','Empirical MSE of aHK0_{hat}','Theoretical Variance of aHK0_{hat}',4);
set(h,'Location','NorthEast')
hold off;

figure; hold on; grid on;
plot(k,mSEaHatUS,'g');
plot(k,aHatUSVar,'b');
plot(k,mSEaHatHK,'r');
plot(k,aHatHKVar,'y');

title('Figure 1: Empirical MSEs and Theoretical Variances')
xlabel('k'); ylabel('MSE & Variance');
h=legend('Empirical MSE of aUS_{hat}','Theoretical Variance of aUS_{hat}','Empirical MSE of aHK_{hat}','Theoretical Variance of aHK_{hat}',4);
set(h,'Location','NorthEast')
hold off;

9. Note that the final theoretical variance was:

$Var(\hat a)= \frac{(1-a^2)^2}{(1+a^2)n}$