[R] Intersection for two curves

[David Winsemius]

On Apr 24, 2010, at 8:10 AM, Muhammad Rahiz wrote:

> 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)}
# Needed to split that line
>   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))

# There is still no function in that call to optimize. Need fc(x,  
a=a,b=b)
# This was pointed out in my last reply:

in1 <- optimize(f=function(x) {abs( fc(x, a=co[1], b=co[2]) -n)}, c(0,  
120) )

>
> abline(v=in1$minimum,col="red",lty=2)
> abline(h=n,col="red",lty=2)
>
> Muhammad
-- 
David.
>
> 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
>>
>>

David Winsemius, MD
West Hartford, CT