Date: Jan 18, 2013 5:23 AM
Author: Torsten
Subject: Re: Maximum Likelihood Estimation and Confidence Intervals
"Topi Kaaresoja" wrote in message <kdb1a9$75m$1@newscl01ah.mathworks.com>...

> Hi,

>

> I have conducted an experiment where users needed to judge where two events were simultaneous or not. The time (milliseconds) between the two events were varied like this:

>

> t=[0, 10, 20, 30, 50, 70, 100, 150, 300];

>

> The question was repeated n=96 times for each value of t. The proportion of "YES" answers is in a response vector

>

> y=[0.875 0.89583 0.89583 0.86458 0.69792 0.63542 0.34375 0.1875 0];

> x=n*y;

>

> I wanted to apply a Gaussian model to that data using maximum likelihood estimation (MLE) similarly than in an MLE Tutorial (http://people.physics.anu.edu.au/~tas110/Teaching/Lectures/L3/Material/Myung03.pdf).

>

> So the PDF of one x is binomial distribution, and my model is

> a*exp(-.5*(((t-mu)/sigma).^2);

>

> So, after some calculations, the negative log-likelihood function becomes:

>

> function loglik = gauss_mle_bin(w,t,y,n)

> % gauss_mle The log-likelihood function of the gaussian model

> % based on binary distribution on separate responses (Myung MLE tutorial)

> % t = vector or m latencies

> % y = vector of m proportion of ?simultaneous? responses

> % Estimates: w(1) = mu, w(2) = sigma, w(3) = a

> x=y*n;

> mu = w(1);

> sigma = w(2);

> a = w(3);

> h = -.5*(((t-mu)/sigma).^2);

> loglik = (-1) * (sum(log(a*exp(h)).*x) + sum(log(1-a*exp(h)).*(n-x)));

>

> I get nice estimates for mu, sigma and a with fminsearch or fminunc.

> 2.4895 77.6230 0.8956

>

> My problem is: How to get confidence intervals for these parameters? Function mle does them automatically, but I have no idea how to pass my data to mle, because my data is not equally spaced.

>

> I know that I can get Hessian matrix from fminunc and by computing an inverse (variance-covariance matrix) and taking square root I will get approximations of standard errors for the parameters. From them I can calculate the 95% confidence intervals as follows:

>

> hessian = hessian returned by fminunc

> varcov = inv(hessian);

> stderrs = sqrt(diag(varcov));

> c= stderrs*sqrt(chi2inv(0.95,1));

>

> I would like to double check if these are correct with mle function (or by other means) if possible.

>

> Best,

>

> Topi Kaaresoja

> Nokia Research Center

> Finland

help nlparci

help nlpredci

Best wishes

Torsten.