[R] Intersection for two curves

[Muhammad Rahiz]

Thanks again David,

I've made the changes. The optimize() function works alright but it does 
not give me the intersection. I suspect it must have to do with the 
user-defined function which I am totally clueless and am inexperienced 
with. Mind showing the way.

Thanks again and thanks for your patience...

*   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

x<- read.table("test.data.txt",header=TRUE)
ds <- x[,2] ; cr <- x[,3]
plot(ds,cr)

fc <- function(x,a,b){a*exp(-b*x)}   
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)

n <- 1/2.71

int <- function(x) {coef(fm)[1] + x*coef(fm)[2]}
in1 <- optimize(f=function(x) abs(int(x)-n),c(0,120))

abline(v=in1$minimum,col="red",lty=2)
abline(h=n,col="red",lty=2)

Muhammad



David Winsemius wrote:
> On Apr 23, 2010, at 8:06 PM, Muhammad Rahiz wrote:
>
>   
>> 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?
>>     
>
> Should work for any process that can produce a function.
>
>   
>> 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'
>>     
>
> Why did you did not change the name of the function?
> And the code below would not have produced that error, so why did you  
> post it?
>
>   
>> 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())
>>     
>
> PLEASE, DON'T DO THAT!!!!. Or rather you can do it in your workspace  
> but don't post it. It's not fair to a person who may not read your  
> code line by line before pasting it into their workspace and having it  
> wiped out. Do you expect us to completely clear out our workspaces  
> just so we can answer your questions? At the very least comment it out  
> so it doesn't blow away half a days work for someone who is trying to  
> be helpful.
>
>   
>> 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))
>>     
>
> Huh? I don't really see anything in there that would be a function. I  
> think you would need to specify "fc" and coefficient values for "a"  
> and "b".
>
>   
>> abline(v=in1$minimize)   # ????
>>     
>
> I think you would use in1$minimum once you get this far.
>
>   
>> 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.
>>>
>>>
>>>       
>
> David Winsemius, MD
> West Hartford, CT
>
>