Hi all,
I'm plotting to get the intersection value of three curves. Defining
the x-axis as dsm, the following code works;
dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0)
s3 = seq(0.05,1.05,0.1)
plot(dsm,s3,col="blue",las=1,ylab="fraction",xlab="distance
(km)")
fc <- function(x,a,b){a*exp(-b*x)}
fm <- nls(s3~fc(dsm,a,b),start=c(a=1,b=0))
co <- coef(fm)
curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1)
r <- range(dsm,na.rm=TRUE)
n <- 1/2.71
val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n),c(r[1],r[2]))
abline(v=val$minimum,lty=2,col="blue",lwd=1)
abline(h=n,lty=2,col="red",lwd=1)
text(100,0.1,paste(round(val$minimum,2),"km",sep=" "))
When I flip the axes, i.e. x-axis = s3, and change the nls start up value
to a=800, the three curves does not intersect.
dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0)
s3 = seq(0.05,1.05,0.1)
plot(s3,dsm,col="blue",las=1,xlab="fraction",ylab="distance
(km)")
fc <- function(x,a,b){a*exp(-b*x)}
fm <- nls(dsm~fc(s3,a,b),start=c(a=800,b=0))
co <- coef(fm)
curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1)
r <- range(dsm,na.rm=TRUE)
n <- 1/2.71
val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n),c(r[1],r[2]))
abline(h=val$minimum,lty=2,col="blue",lwd=1)
abline(v=n,lty=2,col="red",lwd=1)
text(100,0.1,paste(round(val$minimum,2),"km",sep=" "))
It's merely reversing the axes - the code should work, shouldn't it?
Any suggestions why this happens and how I can correct it?
Thanks.
Muhammad
You can do this.
dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0)
s3 = seq(0.05,1.05,0.1)
plot(s3,dsm,col="blue",las=1,xlab="fraction",ylab="distance
(km)")
fc <- function(x,a,b){a*exp(-b*x)}
fm <- nls(dsm~fc(s3,a,b),start=c(a=800,b=0))
co <- coef(fm)
curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1)
r <- range(s3,na.rm=TRUE) # range should be changed
n.dsm <- fc(1/2.71, co[1], co[2]) # calculate a different intersection
val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n.dsm),c(r[1],r[2]))
abline(h=n.dsm,lty=2,col="blue",lwd=1)
abline(v=n,lty=2,col="red",lwd=1)
text(0.1, 500, paste(round(n.dsm,2),"km",sep=" ")) # text
should be placed
at different location
However, note that you do not need `optimize'. The answer that you are
seeking is `n.dsm', which is obtained directly from the nls fit.
Ravi.
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Muhammad Rahiz
Sent: Tuesday, October 19, 2010 5:45 AM
To: x.r-help
Subject: [R] nls & optimize
Hi all,
I'm plotting to get the intersection value of three curves. Defining
the x-axis as dsm, the following code works;
dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0)
s3 = seq(0.05,1.05,0.1)
plot(dsm,s3,col="blue",las=1,ylab="fraction",xlab="distance
(km)")
fc <- function(x,a,b){a*exp(-b*x)}
fm <- nls(s3~fc(dsm,a,b),start=c(a=1,b=0))
co <- coef(fm)
curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1)
r <- range(dsm,na.rm=TRUE)
n <- 1/2.71
val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n),c(r[1],r[2]))
abline(v=val$minimum,lty=2,col="blue",lwd=1)
abline(h=n,lty=2,col="red",lwd=1)
text(100,0.1,paste(round(val$minimum,2),"km",sep=" "))
When I flip the axes, i.e. x-axis = s3, and change the nls start up value
to a=800, the three curves does not intersect.
dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0)
s3 = seq(0.05,1.05,0.1)
plot(s3,dsm,col="blue",las=1,xlab="fraction",ylab="distance
(km)")
fc <- function(x,a,b){a*exp(-b*x)}
fm <- nls(dsm~fc(s3,a,b),start=c(a=800,b=0))
co <- coef(fm)
curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1)
r <- range(dsm,na.rm=TRUE)
n <- 1/2.71
val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n),c(r[1],r[2]))
abline(h=val$minimum,lty=2,col="blue",lwd=1)
abline(v=n,lty=2,col="red",lwd=1)
text(100,0.1,paste(round(val$minimum,2),"km",sep=" "))
It's merely reversing the axes - the code should work, shouldn't it?
Any suggestions why this happens and how I can correct it?
Thanks.
Muhammad
______________________________________________
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.
Let f be your estimated function. Suppose we have a root function, say root(). You are looking for b = root(f-a) where a is some constant. Now suppose we consider the inverse of f, call it f.inv. Then the following holds: a = root(f.inv-b). In your code, you find b = root(f-a) and c = root(f.inv-a). Then you ask, why is b not equal to c? Answer: Why should they? -tgs On Tue, Oct 19, 2010 at 5:44 AM, Muhammad Rahiz < muhammad.rahiz@ouce.ox.ac.uk> wrote:> Hi all, > > I'm plotting to get the intersection value of three curves. Defining the > x-axis as dsm, the following code works; > > dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0) > s3 = seq(0.05,1.05,0.1) > > plot(dsm,s3,col="blue",las=1,ylab="fraction",xlab="distance (km)") > > fc <- function(x,a,b){a*exp(-b*x)} > fm <- nls(s3~fc(dsm,a,b),start=c(a=1,b=0)) > co <- coef(fm) > curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1) > > r <- range(dsm,na.rm=TRUE) > n <- 1/2.71 > val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n),c(r[1],r[2])) > > abline(v=val$minimum,lty=2,col="blue",lwd=1) > abline(h=n,lty=2,col="red",lwd=1) > text(100,0.1,paste(round(val$minimum,2),"km",sep=" ")) > > When I flip the axes, i.e. x-axis = s3, and change the nls start up value > to a=800, the three curves does not intersect. > > dsm = c(800,600,NA,525,NA,450,400,NA,NA,NA,0) > s3 = seq(0.05,1.05,0.1) > > plot(s3,dsm,col="blue",las=1,xlab="fraction",ylab="distance (km)") > > fc <- function(x,a,b){a*exp(-b*x)} > fm <- nls(dsm~fc(s3,a,b),start=c(a=800,b=0)) > co <- coef(fm) > curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="black",lwd=1) > > r <- range(dsm,na.rm=TRUE) > n <- 1/2.71 > val <- optimize(f=function(x) abs(fc(x,a=co[1],b=co[2])-n),c(r[1],r[2])) > > abline(h=val$minimum,lty=2,col="blue",lwd=1) > abline(v=n,lty=2,col="red",lwd=1) > text(100,0.1,paste(round(val$minimum,2),"km",sep=" ")) > > It's merely reversing the axes - the code should work, shouldn't it? > Any suggestions why this happens and how I can correct it? > > Thanks. > > Muhammad > > ______________________________________________ > R-help@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. >[[alternative HTML version deleted]]