[R] nonlinear fitting on many voxels
Charles C. Berry
cberry at tajo.ucsd.edu
Tue Jun 10 04:43:29 CEST 2008
On Mon, 9 Jun 2008, Ayman Oweida wrote:
>
> Below is a sample code using random generated numbers that represents
> what I'm trying to do. I ran it a few times until I got the same
> error. I hope this can help in defining where I'm going
> wrong.
>
The first item under 'Technical Details' in the posting guide says
'No HTML...'
This includes stuff like "<-" (ampersand - ell - tee - ; - minus).
After I editted out your HTML, I triggered your error.
You get the error because you have NaN showing up in the parameter vector.
My suggestion to set options( error = recover ) and try to trace this
through yourself still holds.
Looking at your code I wonder if 'p' needs to obey some constraints. Does
p > 0 need to hold?? p[1] < p[2]??
Perhaps a reparameterization or use of boundary constraints would
straighten things out.
HTH,
Chuck
> library(Bolstad)
> library(nlme)
>
> Ref<-runif(10) # data for reference region
> R1Tn<-runif(10) # data for region of interest
> R10T<-rep(0.1,10) # initial value for each point
> t<- c(0:9) #time
> KR<-0.1 #constant
> VR<-0.1 #constant
>
> f.int<-array(0,length(t)) #array to contain integral values of function f
>
> f<-function(p){
> for (i in 1:(length(Ref)-1)){
>
> f.int[i+1] <- sintegral(t[1:(i+1)], (Ref[1:(i+1)] * (exp((-p[1]/p[2])*(t[i+1]-t)))[1:(i+1)]))
> }
>
> F <- (p[1]/KR)*(Ref) + (p[1]/KR)*( (KR/VR)-(p[1]/p[2]) )* f.int
> F
> }
>
> cf<-array(0,length(t)) #array to contain the first estimated parameter for each voxel
> for (j in 1:length(t)){
> res<- function(p) sum(((R1Tn[j]-R10T[j])-f(p))^2)
>
> out <- ( nlminb(c(0.1,0.25),res) )
>
> cf[j]<-out$par[1] # fill the 1st estimated parameter for each voxel into cf
> }
>
> Thank you,
> Ayman
> --- On Mon, 6/9/08, Charles C. Berry <cberry at tajo.ucsd.edu> wrote:
> From: Charles C. Berry <cberry at tajo.ucsd.edu>
> Subject: Re: [R] nonlinear fitting on many voxels
> To: "Ayman Oweida" <ajoweida at yahoo.com>
> Cc: r-help at r-project.org
> Date: Monday, June 9, 2008, 12:59 PM
>
> On Mon, 9 Jun 2008, Ayman Oweida wrote:
>
> > After many months, I am now banging my head against the wall because I
> can't find a solution to this seemingly trivial problem.&nbsp; Any help
> would be appreciated:
> >
> > I am trying to apply a nonlinear fitting routine to a 3D MR image on a
> voxel-by-voxel basis.&nbsp; I've tested the routine using simulated
> data and things went well.&nbsp; As for the real data, the fitting routine
> works variably.&nbsp; By variably, I mean the following: &nbsp; with a
> specific set of starting parameters the routine would work for the first 10
> voxels, for another set of starting values the routine would work for the first
> 1000 voxels and so on.&nbsp; NEVER was I able to get starting values that
> would allow the routine to run entirely!&nbsp; The error I would get after
> fitting the limited number of voxels is:&nbsp;
> >
> > "Error in approx(x, fx, n = 2 * n.pts + 1) :
> > &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
> need at least two non-NA values to interpolate"
> >
> > I think the error is from the library Bolstad since I'm using the
> > sintegral function as part of the fitting equation.&nbsp; I've
> tried
> > numerous nonlinear functions including: nlm, optim and nlminb, but all
> > stop after performing the fit on a limited number of voxels.&nbsp; I
> > took the values for which the fitting routine stops and I applied
> > different starting values to them and it works!&nbsp; I'm still
> using
> > the same voxel/parameter values and I am sure they are not NA values as
> > per the error!
>
> Reread the error message. It does NOT say that there are NA values. It
> only says that you 'need at least two non-NA values'.
>
> >
> > Ok, I think the problem is clear.&nbsp; any solutions?
>
> The problem is not clear.
>
> Here is some advice:
>
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented,
> minimal, self-contained, reproducible code.
>
> Often, developing the 'minimal' example helps you to perceive the
> underlying difficulty. Without it, you will only get very general advice
> and off-hand guesses about where your problem lies, which may not move you
> closer to a solution.
>
> ---
>
> So here is some of that general advice: find out why approx() sends you
> that error message.
>
> One way to do this is to set
>
> options(error=recover)
>
> before running your function. And then inspect objects in the frame in
> which the error occurred and in the frames leading up to the one in which
> the error was triggered. See
>
> ?recover
>
> HTH,
>
> Chuck
>
>
> &nbsp; I am new to this complex world of statistical analysis.
> > Your help is very very much appreciated.
> >
> > Thanks,
> > Ayman
> >
> > Montreal Neurological Institute
> > Montreal, Canada
> >
> >
> >
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help at r-project.org 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.
> >
>
> Charles C. Berry (858) 534-2098
> Dept of Family/Preventive Medicine
> E mailto:cberry at tajo.ucsd.edu UC San Diego
> http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org 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.
>
Charles C. Berry (858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
More information about the R-help
mailing list