Path: news.mathworks.com!not-for-mail
From: Peter Perkins <Peter.PerkinsRemoveThis@mathworks.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: glm with probit model
Date: Mon, 10 Nov 2008 10:33:43 -0500
Organization: The MathWorks, Inc.
Lines: 38
Message-ID: <gf9k8o$e3s$1@fred.mathworks.com>
References: <geapba$62l$1@fred.mathworks.com> <gecis6$hje$1@fred.mathworks.com> <gf6c79$87j$1@fred.mathworks.com>
NNTP-Posting-Host: perkinsp.dhcp.mathworks.com
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
X-Trace: fred.mathworks.com 1226331224 14460 172.31.57.88 (10 Nov 2008 15:33:44 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Mon, 10 Nov 2008 15:33:44 +0000 (UTC)
User-Agent: Thunderbird 2.0.0.17 (Windows/20080914)
In-Reply-To: <gf6c79$87j$1@fred.mathworks.com>
Xref: news.mathworks.com comp.soft-sys.matlab:499979


Daniel wrote:

> If I understand the m-file you provided, the following represent the three functions that are required. 
> deta=(1-2*lambda)*normpdf(x,mu,sigma);    %This is the derivative of the link function
> eta=norminv((p_lapse-lambda)/(1-2*lambda),mu,sigma);     %This is a form of the link function
> mu=lambda + (1-2*lambda)*normcdf(x,mu,sigma);  %This is a form of the inverse link function
> 
> However, when I try to use the m-file that you provided as a template, I start to run into problems.
> 
> I can almost figure out how to implement the derivative function, since I can write this equation in a simple closed form solution:
> deta=(1-2*L)*exp(-((mu)).^2)/sqrt(2*pi);  
> (But this is missing the normalization by sigma which I can't see how to write in the desired format for the GLM, but I probably am just missing something.)
> 
> However, the inverse link function cannot be written in a closed form solution.  In fact it includes the normal cdf, which is a from of the error function and I just can't see how to write this in a closed form that is appropriate for the glmfit.  Seems like the same is true for the link function.

I'm not sure where those 2's are coming from.

These are the probit link, link derivative, and inverse link functions that GLMFIT uses (edit private/stattestlink, called by glmfit):

    link = @(mu) norminv(mu);
    dlink = @(mu) 1 ./ normpdf(norminv(mu));
    % keep mu = ilink(eta) in [eps, (1-eps)];
    lowerBnd = norminv(eps(dataClass)); upperBnd = -lowerBnd;
    ilink = @(eta) normcdf(min(max(eta,lowerBnd),upperBnd));

Given lambda, I believe you want 

    ilink = @(eta) lambda + ...
            (1-lambda).*normcdf(min(max(eta,lowerBnd),upperBnd));

and therefore

    link = @(mu) norminv((mu-lambda)./(1-lambda));
    dlink = @(mu) 1 ./ ((1-lambda).*normpdf(norminv((mu-lambda)./(1-lambda))));

but you'll want to check my work.  If your lambda is really .01, then you may not run into the problem with initial values that I alluded to, unless youhave more than 100 observations.

Hope this helps.