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: Tue, 11 Nov 2008 11:26:47 -0500
Organization: The MathWorks, Inc.
Lines: 25
Message-ID: <gfcbo7$sfi$1@fred.mathworks.com>
References: <geapba$62l$1@fred.mathworks.com> <gecis6$hje$1@fred.mathworks.com> <gf6c79$87j$1@fred.mathworks.com> <gf9k8o$e3s$1@fred.mathworks.com> <gf9oam$jvi$1@fred.mathworks.com> <gface8$ka3$1@fred.mathworks.com> <gfc0ct$cnl$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 1226420807 29170 172.31.57.88 (11 Nov 2008 16:26:47 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 11 Nov 2008 16:26:47 +0000 (UTC)
User-Agent: Thunderbird 2.0.0.17 (Windows/20080914)
In-Reply-To: <gfc0ct$cnl$1@fred.mathworks.com>
Xref: news.mathworks.com comp.soft-sys.matlab:500249


Daniel wrote:
> Hi Peter , thanks for all of your help.  it seems that I am very close but I seem to missing something critical.  Following your suggestion, I defined the following functions:
> 
> myilink = @(eta) L + (1-2*L).*normcdf(min(max(eta,lowerBnd),upperBnd));
> mylink = @(mu) norminv((mu-L)./(1-2*L));
> mydlink = @(mu) 1 ./ ((1-2*L).*normpdf(norminv((mu-L)./(1-L))));

I think you want 1-2*L in the denom of the last one.


> However, when I fit using the following call (as I believe you had suggested):
> b = glmfit(motion,y,'binomial','link',{mylink mydlink myilink});
> I get the following error message:
> 
> ??? Input argument "L" is undefined.
> 
> Error in ==> direction_detect>@(mu,L)norminv((mu-L)./(1-2*L)) at 64
> mylink = @(mu,L) norminv((mu-L)./(1-2*L));

Where did that link function come from?  It isn't the one you cited above, and it has L as an argument.  L should not be an input argument.  Let the anonomous function "captures" L by creating a new mylink at each iteration inside this loop:

> lambda=0:.01:.06;
> for j=1:length(lambda)
>     L = lambda(j);