[R] Idea/package to "linearize a curve" along the diagonal?

Jeff Newmiller jdnewmil at dcn.davis.ca.us
Mon Mar 12 22:56:24 CET 2012


You said:
> Does this make sense?

Not to me. I don't understand either of your points. It may be time for you to provide a concrete example of your own (reproducible R code, or at least specific input and output data).
---------------------------------------------------------------------------
Jeff Newmiller                        The     .....       .....  Go Live...
DCN:<jdnewmil at dcn.davis.ca.us>        Basics: ##.#.       ##.#.  Live Go...
                                      Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
/Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k
--------------------------------------------------------------------------- 
Sent from my phone. Please excuse my brevity.

Emmanuel Levy <emmanuel.levy at gmail.com> wrote:

>Hi Jeff,
>
>Thanks for your reply and the example.
>
>I'm not sure if it could be applied to the problem I'm facing though,
>for two reasons:
>
>(i) my understanding is that the inverse will associate a new Y
>coordinate given an absolute X coordinate. However, in the case I'm
>working on, the transformation that has to be applied depends on X
>*and* on its position relative to the *normal* of the fitted curve.
>This means, for instance, that both X and Y will change after
>transformation.
>
>(ii) the fitted curve can be described by a spline, but I'm not sure
>if inverse of such models can be inferred automatically (I don't know
>anything about that).
>
>The procedure I envision is the following: treat the curve "segment by
>segment", apply rotation+translation to each segment to bring it on
>the
>diagonal, and apply the same transformation to all points
>corresponding to the same segment (i.e., these are the points that are
>close and within the "normal" area covered by the segment).
>
>Does this make sense?
>All the best,
>
>Emmanuel
>
>
>On 12 March 2012 02:15, Jeff Newmiller <jdnewmil at dcn.davis.ca.us>
>wrote:
>> It is possible that I do not see what you mean, but it seems like the
>following code does what you want. The general version of this is
>likely to be rather more difficult to do, but in theory the inverse
>function seems like what you are trying to accomplish.
>>
>> x <- 1:20
>> y <- x^2 + rnorm(20)
>>
>> y.lm <- lm( y ~ I(x^2) + x )
>> plot( x, y )
>> lines( x, predict( y.lm ) )
>>
>> qa <- coef(y.lm)["I(x^2)"]
>> qb <- coef(y.lm)["x"]
>> qc <- coef(y.lm)["(Intercept)"]
>>
>> # define inverse of known model
>> f1 <- function( y ) { ( sqrt( 4*qa*( y -qc ) + qb^2 ) - qb ) / ( 2*qa
>) }
>> f2 <- function( y ) { -( sqrt( 4*qa*( y -qc ) + qb^2 ) + qb ) / (
>2*qa ) }
>>
>> plot( x, f1( y ) )
>>
>>
>>
>---------------------------------------------------------------------------
>> Jeff Newmiller                        The     .....       .....  Go
>Live...
>> DCN:<jdnewmil at dcn.davis.ca.us>        Basics: ##.#.       ##.#.  Live
>Go...
>>                                      Live:   OO#.. Dead: OO#..
> Playing
>> Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
>> /Software/Embedded Controllers)               .OO#.       .OO#.
> rocks...1k
>>
>---------------------------------------------------------------------------
>> Sent from my phone. Please excuse my brevity.
>>
>>
>>
>> Emmanuel Levy <emmanuel.levy at gmail.com> wrote:
>>
>>>Dear Jeff,
>>>
>>>I'm sorry but I do not see what you mean. It'd be great if you could
>>>let me know in more details what you mean whenever you can.
>>>
>>>Thanks,
>>>
>>>Emmanuel
>>>
>>>
>>>On 12 March 2012 00:07, Jeff Newmiller <jdnewmil at dcn.davis.ca.us>
>>>wrote:
>>>> Aren't you just reinventing the inverse of a function?
>>>>
>>>---------------------------------------------------------------------------
>>>> Jeff Newmiller                        The     .....       .....  Go
>>>Live...
>>>> DCN:<jdnewmil at dcn.davis.ca.us>        Basics: ##.#.       ##.#.
> Live
>>>Go...
>>>>                                      Live:   OO#.. Dead: OO#..
>>> Playing
>>>> Research Engineer (Solar/Batteries            O.O#.       #.O#.
> with
>>>> /Software/Embedded Controllers)               .OO#.       .OO#.
>>> rocks...1k
>>>>
>>>---------------------------------------------------------------------------
>>>> Sent from my phone. Please excuse my brevity.
>>>>
>>>> Emmanuel Levy <emmanuel.levy at gmail.com> wrote:
>>>>
>>>>>Hi,
>>>>>
>>>>>I am trying to normalize some data. First I fitted a principal
>curve
>>>>>(using the LCPM package), but now I would like to apply a
>>>>>transformation so that the curve becomes a "straight diagonal line"
>>>on
>>>>>the plot.  The data used to fit the curve would then be normalized
>by
>>>>>applying the same transformation to it.
>>>>>
>>>>>A simple solution could be to apply translations only (e.g., as
>done
>>>>>after a fit using loess), but here rotations would have to be
>applied
>>>>>as well. One could visualize this as the "stretching of a curve",
>>>>>i.e., during such an "unfolding" process, both translations and
>>>>>rotations would need to be applied.
>>>>>
>>>>>Before I embark on this (which would require me remembering long
>>>>>forgotten geometry principles) I was wondering if you can think of
>>>>>packages or tricks that could make my life simpler?
>>>>>
>>>>>Thanks for any input,
>>>>>
>>>>>Emmanuel
>>>>>
>>>>>
>>>>>Below I provide an example - the black curve is to be "brought"
>along
>>>>>the diagonal, and all data points normal to a "small segment" (of
>the
>>>>>black curve) would undergo the same transformation as it - what
>>>>>"small" is remains to be defined though.
>>>>>
>>>>>    tmp=rnorm(2000)
>>>>>    X.1 = 5+tmp
>>>>>    Y.1 = 5+ (5*tmp+rnorm(2000))
>>>>>    tmp=rnorm(1000)
>>>>>    X.2 = 9+tmp
>>>>>    Y.2 = 40+ (1.5*tmp+rnorm(1000))
>>>>>    X.3 = 7+ 0.5*runif(500)
>>>>>    Y.3 = 15+20*runif(500)
>>>>>    X = c(X.1,X.2,X.3)
>>>>>    Y = c(Y.1,Y.2,Y.3)
>>>>>
>>>>>    lpc1 = lpc(cbind(X,Y), scaled=FALSE, h=c(1,1) )
>>>>>    plot(lpc1)
>>>>>
>>>>>______________________________________________
>>>>>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.
>>>>
>>



More information about the R-help mailing list