[R] Beginners Question: Make nlm work
Johannes Graumann
graumann at its.caltech.edu
Thu Aug 26 03:57:55 CEST 2004
Aaaahhhh - there was my conceptual problem ... Thanks!
... but 'Oh, God!'! R sucks at fitting (when directed to do so by the
simple minded)! Or am I really that off?
I attach a plot that contains my data as well as two modeled curves:
- the nice one fitted by gnuplot (nonlinear least-squares (NLLS)
Marquardt-Levenberg)
- the nasty one fitted by the code below.
Both algorithms were started with the same values for p[1] and p[2].
Comments?
Is there any way to access the 2 elements of 'out$estimate' from the
program?
Joh
setwd("~/Biology/R_versuch")
mydata<-read.table("YJG45-7_Growth.dat")
#plot(mydata$V1,mydata$V2,xlab="Time
(h)",ylab=expression(OD[600][~nm]),las=1) x<-mydata$V1
y<-mydata$V2
VH <- function(p) sum(
(y - p[1]/(1+((p[1]-0.008)/0.008)*exp(-(p[2]*x))))^2)
plot(x, y,
axes=FALSE,
type="p",
xlab="Time (min)",
ylab=expression(OD[600][~nm]),
las=1,
pch=20,
tck=0.015,
mgp=c(1.5,0.25,0),
)
axis(1,
at=NULL,
tick=TRUE,
tck=0.015,
mgp=c(1.5,0.25,0)
)
axis(2,
at=NULL,
tick=TRUE,
tck=0.015,
mgp=c(1.5,0.25,0)
)
axis(3,
at=c(300,600,900,1200,1500,1800,2100),
labels=c(5,10,15,20,25,30,35),
tick=TRUE,
tck=0.015,
mgp=c(1.5,0.25,0))
axis(4,
labels=FALSE,
at=NULL,
tick=TRUE,
tck=0.015,
mgp=c(1.5,0.25,0)
)
box()
out <- nlm(VH, p = c(3, 4e-3), hessian = TRUE)
#Plot Results from Gnuplot fit
yfit1 <- 4.21632/(1+((4.21632-0.008)/0.008)*exp(-(0.00414403*x)))
lines(spline(x, yfit1))
#Plot Results from 'out'
yfit2 <-
3.000002050/(1+((3.000002050-0.008)/0.008)*exp(-(0.004587924*x)))
lines(spline(x,yfit2))
On Wed, 25 Aug 2004 21:06:24 -0400
"Jim Brennan" <jfbrennan at rogers.com> wrote:
>
> ----- Original Message -----
> From: "Johannes Graumann" <johannes_graumann at web.de>
> To: <r-help at stat.math.ethz.ch>
> Sent: Wednesday, August 25, 2004 7:07 PM
> Subject: [R] Beginners Question: Make nlm work
>
>
> > Hello,
> >
> > I'm new to this and am trying to teach myself some R by plotting
> > biological data. The growth curve in question is supposed to be
> > fitted to the Verhulst equation, which may be transcribed as
> > follows: f(x)=a/(1+((a-0.008)/0.008)*exp(-(b*x)))
> > - for a known population density (0.008) at t(0).
> >
> > I am trying to rework the example from "An Introduction to R" (p.
> > 72) for my case and am failing miserably. Could somebody glance over
> > the code below and nudge me into the right direction - I must have
> > some major conceptual problem which can't be solved by me staring at
> > it ... Since I'm repeating something I have done with gnuplot I know
> > that 3 and 4e-3 as starting values for the fit are appropriate ...
> >
> > Thanks for any hint,
> >
> > Joh
> >
> > setwd("~/Biology/R_versuch")
> > mydata<-read.table("YJG45-7_Growth.dat")
> > x<-mydata$V1
> > y<-mydata$V2
> > VH <- function(p) y ~ p[1]/(1+((p[1]-0.008)/0.008)*exp(-(p[2]*x)))
>
> Shouldn't this be
> VH <- function(p) sum((y -
> p[1]/(1+((p[1]-0.008)/0.008)*exp(-(p[2]*x))))^2) Where you are
> minimizing the squared deviations from what your given equation is
> when you pass it to nlm?
>
> Maybe if you send some data that would make it more clear.
>
> Jim
>
>
> > plot(x, y, xlab="Time (h)",ylab=expression(OD[600][~nm]),las=1)
> > out <- nlm(VH, p = c(3, 4e-3), hessian = TRUE)
> >
> > ______________________________________________
> > R-help at 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
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plot1.eps
Type: application/postscript
Size: 6160 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20040825/acb8a8b9/plot1.eps
More information about the R-help
mailing list