% function [dbar,tval,P,MLE,likelihood]=

Discussion in 'C Programming' started by Martin, Mar 30, 2011.

  1. Martin

    Martin Guest

    % phyl_paired_ttest(C, x1, SEx1, x2, SEx2,MinMaxeps,MinMaxlam)
    % Takes C (tree matrix), x1 (n x 1 vector), SEx1 (n x 1 vector)
    % x2, SEx2, [MIN MAX] values for grid search on epsilon (eps) and
    lambda
    % (lam).
    % Returns dbar, tval, P (P-value), MLEs for epsilon, sigma^2, and
    lambda,
    % and the log-likelihood.
    % Optimization is performed on a 100 x 100 grid.
    function [dbar,tval,P,MLE,likelihood]=phyl_paired_ttest(C, x1, SEx1,
    x2,
    SEx2,MinMaxeps,MinMaxlam)
    % set grid search criterion
    STEP=0.01; % *100%
    % number of parameters
    np=3;
    % number of taxa
    n=max(size(C));
    % compute differences
    d=x1-x2;
    % compute a matrix of estimated sampling variances
    E=diag(SEx1.^2+SEx2.^2);
    % create column vector of 1.0s
    one=ones(n,1);
    % set grid search for epsilon & lambda
    lambda=MinMaxlam(1):STEP*(MinMaxlam(2)-MinMaxlam(1)):MinMaxlam(2);
    epsilon=MinMaxeps(1):STEP*(MinMaxeps(2)-MinMaxeps(1)):MinMaxeps(2);
    % determine if the user has fixed epsilon or lambda
    if (MinMaxeps(2)-MinMaxeps(1))==0
    epsilon=MinMaxeps;
    np=np-1;
    end
    if (MinMaxlam(2)-MinMaxlam(1))==0
    lamda=MinMaxlam;
    np=np-1;
    end
    % OR uncomment to fix lambda at one or zero or epsilon at zero
    % lambda=[0 0]; np=np-1;
    % lambda=[1 1]; np=np-1;
    % epsilon=[0 0]; np=np-1;
    % ok, now estimate lambda & epsilon by maximizing the likelihood
    maxL=-Inf; minL=Inf;
    logL=ones(max(size(lambda)),max(size(epsilon)))*(-Inf);
    for i=1:max(size(lambda))
    for j=1:max(size(epsilon))
    % transform by lambda
    Cl=lambda_transform(lambda(i),C);
    % add sampling error
    V=Cl+epsilon(j)*E;
    % check positive definiteness
    [R p]=chol(V);
    if(p==0)
    % compute parameter estimates and likelihood
    dbar=(one'*V^-1*one)^-1*(one'*V^-1*d);
    sigma2=(d-dbar*one)'*V^-1*(d-dbar*one)*n^-1;
    logL(i,j)=-(1/2)*(d-dbar*one)'*(sigma2*V)^-1*(d-dbar*one)...
    -(1/2)*log(det(sigma2*V))-(n/2)*log(2*pi);
    if logL(i,j)>maxL
    MLpos=[i j];
    maxL=logL(i,j);
    MLE.sigma2=sigma2;
    end
    if logL(i,j)<minL
    minL=logL(i,j);
    minLpos=[i j];
    end
    end
    end
    end
    % create plot of likelihood surfaces
    colormap(sort(gray,'descend'));
    eps=char(hex2dec('65')); lam=char(hex2dec('6c'));
    subplot(2,2,1); mesh(epsilon,lambda,logL);
    xlabel(eps,'fontname','symbol'); ylabel(lam,'fontname','symbol');
    zlabel('log(L)');
    title('A)
    ');
    subplot(2,2,3); plot(lambda,logL:),MLpos(2)),'k');
    xlabel(lam,'fontname','symbol'); ylabel('log(L)');
    title('C)
    ');
    subplot(2,2,4); plot(epsilon,logL(MLpos(1),:),'k');
    xlabel(eps,'fontname','symbol'); ylabel('log(L)');
    title('D)
    ');
    % replace -Inf with minL
    for i=1:max(size(lambda))
    for j=1:max(size(epsilon))
    if logL(i,j)==-Inf
    logL(i,j)=minL;
    end
    end
    end
    subplot(2,2,2); contour(epsilon,lambda,logL,2000);
    xlabel(eps,'fontname','symbol'); ylabel(lam,'fontname','symbol');
    title('B)
    ');
    suptitle('Likelihood Surfaces');
    MLE.lambda=lambda(MLpos(1));
    MLE.epsilon=epsilon(MLpos(2));
    likelihood=maxL;
    % pause
    pause(1);
    % ok, now perform calculations using MLlambda & MLepsilon
    Cl=lambda_transform(MLE.lambda,C);
    % calculate mean difference
    V=Cl+MLE.epsilon*E;
    dbar=(one'*V^-1*one)^-1*(one'*V^-1*d);
    % calculate SE of dbar
    SEdbar=sqrt(MLE.sigma2*(n/(n-np))*(one'*V^-1*one)^-1);
    MLE.sig2e=MLE.sigma2*MLE.epsilon;
    % calculate t-statistic
    tval=dbar*SEdbar^-1;
    % calculate P-value
    P=2*(1-tcdf(abs(tval),n-np));
    % print results to screen
    fprintf(1,'dbar = %f\tt(df = %d) = %f\t',dbar,n-np,tval);
    fprintf(1,'P(t) = %f\n',P);
    fprintf(1,'MLE(lambda) = %f\tlog(L) = %f\n',MLE.lambda,maxL);
    fprintf(1,'MLE(epsilon) = %f\tlog(L) = %f\n',MLE.epsilon,maxL);
    fprintf(1,'MLE(sig2e) = %f\tlog(L) = %f\n',MLE.sig2e,maxL);
    % done
    % function Ml=lambda_transform(lambda,M)
    % Multiplies off-diagonals of square matrix M by lambda
    % and returns Ml.
    function Ml=lambda_transform(lambda,M)
    Ml=M;
    n=max(size(M));
    % transform by lambda
    for i=1:n
    for j=(i+1):n
    Ml(i,j)=M(i,j)*lambda;
    Ml(j,i)=Ml(i,j);
    end
    end
    % done
    Martin, Mar 30, 2011
    #1
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. penny336
    Replies:
    1
    Views:
    655
    tom_usenet
    Nov 7, 2003
  2. James Vanns
    Replies:
    7
    Views:
    7,010
    Evan Carew
    Jan 21, 2004
  3. komal
    Replies:
    6
    Views:
    1,404
    msalters
    Jan 25, 2005
  4. Mike P

    MLE in Python for distributions

    Mike P, Apr 3, 2008, in forum: Python
    Replies:
    0
    Views:
    397
    Mike P
    Apr 3, 2008
  5. Replies:
    10
    Views:
    207
    Oscar Benjamin
    Feb 2, 2013
Loading...

Share This Page