[R] fit a nonlinear model using nlm()
I am trying to fit a nonlinear model using nlm().
My application of nlm() is a bit complicated.
Here is the story behind the model being fit:
The observer is trying to detect a signal corrupted by noise.
On each trial, the observer gets stim=signal+rnorm().
In the simulation below I have 500 trials. Each row of stim is a new trial.
On each trial, if the cross-correlation between the stim and
the signal is above some criterion level (crit=.5 here), the
observer says "signal" (resp=1), else he says "no signal"
(resp=0).
So the situation is this: I know resp and stim for each trial.
I would like to estimate crit and signal from these.
(You might say that I already know the signal and crit. But the
idea here is that the observer cross-correlates the stim with an
internal template that will not be identical to the actual signal.
I want to estimate this template. Also the observer's crit may
differ from the "correct" one.)
In the code below, please help me get the f() and nlm() bits right.
I want to estimate signal and crit given stim and resp.
Thanks very much for any help!
Bill
x<-1:100
con<-.1
signal<-con*cos(2*pi*3*x/length(x))
crit<-.5
noisesd<-.1
# each row is a new stim (trial). 500 trials
resp<-array(dim=500)
stim<-matrix(nrow=500,ncol=100)
for (i in 1:500)
{
stim[i,]<-signal+rnorm(n=length(signal), mean=0, sd=noisesd)
if (sum(stim[i,]*signal)>crit) (resp[i]<-1) else (resp[i]<-0)
}
f<-function(signalest)
{
r<-array(dim=500)
for (i in 1:500)
{
r[i]<-sum(stim[i,]*signalest)>critest
}
sum((r-resp)^2)
}
fit<-nlm(f, stim[1,])
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] fit a nonlinear model using nlm()
William Simpson gmail.com> writes:
>
> I am trying to fit a nonlinear model using nlm().
> The observer is trying to detect a signal corrupted by noise.
> On each trial, the observer gets stim=signal+rnorm().
>
> In the simulation below I have 500 trials. Each row of stim is a new trial.
> On each trial, if the cross-correlation between the stim and
> the signal is above some criterion level (crit=.5 here), the
> observer says "signal" (resp=1), else he says "no signal"
> (resp=0).
>
> Thanks very much for any help!
> Bill
>
>
It sounds like you are doing a classification image experiment.
You can use tapply() to get means for each x as a function of the
observer's classifications and then combine them as a function of
hits, false alarms, misses, correct rejections using the weights
1 -1, -1, 1, as in Ahumada's original approach.
You can do this with lm() if you set it up so that the noise is
the response and the classifications are a 4 level factor that prediccts
them and the contrasts are set up as above.
I think that it would be better to set it up as a glm, however,
where the responses of the observer are binary responses, the
noise and the presence/absence of the signal are predictor variables,
with a binomial family.
I have example code that does each of these if you are
interested and if it is only for simulation, see
@Article{pmid16277294,
Author="Thomas, James P and Knoblauch, Kenneth",
Title="{{F}requency and phase contributions to the detection of
temporal luminance modulation}",
Journal="J Opt Soc Am A Opt Image Sci Vis",
Year="2005",
Volume="22",
Number="10",
Pages="2257--2261",
Month="Oct"
}
for which I can send you the code also.
Ken
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
