[R] Intersection for two curves

[Muhammad Rahiz]

Thanks David & Peter,

The locator() works but not practical as I have to repeat the process 
many times.

Does the code works on linear regression only? When i tried to find the 
intersection at a non-linear curve, i get the following error
Error in optimize(f = function(x) abs(xyf(ds) - n), c(0, 13)) :
   invalid function value in 'optimize'

I'd like to know what the error means and how I can correct it.

I have my sample data and code as follows;


*   s1         cm
1    0 0.57223196
2   10 0.33110049
3   20 0.11163181
4   30 0.10242237
5   40 0.09254315
6   50 0.02739370
7   60 0.02567137
8   70 0.02492397
9   80 0.03206637
10  90 0.02487381
11 100 0.01747584
12 110 0.15977533
13 120 0.13317708

rm(list=ls())
x<- read.table("test.data.txt",header=TRUE)

ds <- x[,2] # distance
cr <- x[,3] # correlation values
plot(ds,cr)
n <- 1/2.71
abline(h=n)

fc <- function(x,a,b){a*exp(-b*x)}    # where a & b are constants
fm <- nls(cr ~ fc(ds,a,b),start=c(a=1,b=0))
co <- coef(fm)
curve(fc(x,a=co[1],b=co[2]),add=TRUE,col="red",lwd=1.5)

int <- function(x) coef(fm)[1] + x*coef(fm)[2]
in1 <- optimize(f=function(x) c(0,120), abs(int(ds)-n))
abline(v=in1$minimize)

thanks,
Muhammad




On 04/23/2010 07:32 PM, Peter Ehlers wrote:
> On 2010-04-23 11:46, David Winsemius wrote:
>    
>> On Apr 23, 2010, at 1:06 PM, Muhammad Rahiz wrote:
>>
>>      
>>> Does anyone know of a method that I can get the intersection where the
>>> red and blue curves meet i.e. the value on the x-axis?
>>>
>>> x<- 1:10
>>> y<- 10:1
>>> plot(x,y)
>>> abline(lm(y~x),col="blue")
>>> abline(h=2.5,col="red")
>>>        
>> Two ways :
>>
>>   >  xy<- lm(y~x)
>>   >  xyf<- function(x) coef(xy)[1] +x*coef(xy)[2]
>>
>> # absolute difference
>>   >  optimise(f=function(x) abs(xyf(x)-2.5), c(1,10) )
>> $minimum
>> [1] 8.49998
>>
>> $objective
>> (Intercept)
>> 1.932015e-05
>>
>> #N minimize squared difference
>>   >  optimise(f=function(x) (xyf(x)-2.5)^2, c(1,10) )
>> $minimum
>> [1] 8.5
>>
>> $objective
>> (Intercept)
>> 3.155444e-30
>>
>>      
> Another (crude) way is to use locator(). I usually maximize
> the plot window for this.
>
>