Results 1 to 9 of 9

Math Help - Heavy duty MLE

  1. #1
    Super Member Anonymous1's Avatar
    Joined
    Nov 2009
    From
    Big Red, NY
    Posts
    517
    Thanks
    1

    Talking 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.
    Last edited by Anonymous1; March 17th 2010 at 08:44 PM.
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Super Member Anonymous1's Avatar
    Joined
    Nov 2009
    From
    Big Red, NY
    Posts
    517
    Thanks
    1
    I'm working out the answer now. If any one is curious I will post my solution.
    Follow Math Help Forum on Facebook and Google+

  3. #3
    MHF Contributor matheagle's Avatar
    Joined
    Feb 2009
    Posts
    2,763
    Thanks
    5
    I assume you want the MLE of a.
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Super Member Anonymous1's Avatar
    Joined
    Nov 2009
    From
    Big Red, NY
    Posts
    517
    Thanks
    1
    Yes I do. It's tedious to calculate, but can you see why this problem is interesting?
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Moo
    Moo is offline
    A Cute Angle Moo's Avatar
    Joined
    Mar 2008
    From
    P(I'm here)=1/3, P(I'm there)=t+1/3
    Posts
    5,618
    Thanks
    6
    Well, there are formulas to estimate a covariance matrix... so it should be possible to narrow it to a huh ?
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Super Member Anonymous1's Avatar
    Joined
    Nov 2009
    From
    Big Red, NY
    Posts
    517
    Thanks
    1
    Did I hear someone thinking maximization of a cubic?
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Flow Master
    mr fantastic's Avatar
    Joined
    Dec 2007
    From
    Zeitgeist
    Posts
    16,948
    Thanks
    5
    Quote Originally Posted by Anonymous1 View Post
    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/
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Super Member Anonymous1's Avatar
    Joined
    Nov 2009
    From
    Big Red, NY
    Posts
    517
    Thanks
    1
    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');
    u2= load('STATES.txt');
    u3= load('HONG.txt');
    u4= load('kONG.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;
    Last edited by Anonymous1; April 2nd 2010 at 11:45 AM.
    Follow Math Help Forum on Facebook and Google+

  9. #9
    Super Member Anonymous1's Avatar
    Joined
    Nov 2009
    From
    Big Red, NY
    Posts
    517
    Thanks
    1
    Note that the final theoretical variance was:

    Var(\hat a)= \frac{(1-a^2)^2}{(1+a^2)n}
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Calculus heavy physics problem
    Posted in the Calculus Forum
    Replies: 6
    Last Post: October 11th 2011, 02:05 PM
  2. Heavy uniform sphere - angle of friction?
    Posted in the Advanced Applied Math Forum
    Replies: 0
    Last Post: June 15th 2009, 01:34 PM

Search Tags


/mathhelpforum @mathhelpforum