Date: Jan 18, 2013 5:23 AM
Subject: Re: Maximum Likelihood Estimation and Confidence Intervals
"Topi Kaaresoja" wrote in message <firstname.lastname@example.org>...
> 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];
> 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
> 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
> 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.
> Topi Kaaresoja
> Nokia Research Center