[R] finding interpolated values along an empirical parametric curve

Michael Friendly friendly at yorku.ca
Mon Dec 5 17:02:07 CET 2011


Thanks, John
That's very clever.  I didn't realize that the splines package supports 
multivariate regression
splines and that works well enough for my purposes; plus this solution 
is very transparent.

best,
-Michael


On 12/5/2011 9:50 AM, John Fox wrote:
> Hi Michael,
>
> I can get what appears to be a good interpolation with a regression spline
> in a multivariate LM, playing around with the tuning parameter to leave 1
> residual df. Try this:
>
> library(splines)
> mod<- lm(cbind(log.det, norm.beta) ~ bs(lambda, df=4), data=pd)
> summary(mod)
>
> x<- data.frame(lambda=seq(min(pd$lambda), max(pd$lambda), length=100))
> fit<- predict(mod, newdata=x)
> points(fit[, "norm.beta"], fit[, "log.det"], pch=16, cex=0.5)
>
> x.2<- data.frame(lambda=c(lambda.HKB, lambda.LW))
> fit.2<- predict(mod, x.2)
> points(fit.2[, "norm.beta"], fit.2[, "log.det"], pch=15, col="green")
>
> That doesn't solve the problem of calculating the normal, however.
>
> I hope this helps,
>   John
>
>> -----Original Message-----
>> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
>> project.org] On Behalf Of Michael Friendly
>> Sent: December-05-11 9:16 AM
>> To: R-help
>> Subject: [R] finding interpolated values along an empirical parametric
>> curve
>>
>> Given the following data, I am plotting log.det ~ norm.beta, where the
>> points depend on a parameter, lambda (but there is no functional form).
>> I want to find the (x,y) positions along this curve corresponding to
>> two special values of lambda
>>
>> lambda.HKB<- 0.004275357
>> lambda.LW<- 0.03229531
>>
>> and draw reference lines at ~ -45 degrees (or normal to the curve) thru
>> these points.
>> How can I do this?  A complete example is below
>>
>>   >  pd
>>         lambda   log.det norm.beta
>> 0.000  0.000 -12.92710  3.806801
>> 0.005  0.005 -14.41144  2.819460
>> 0.010  0.010 -15.41069  2.423197
>> 0.020  0.020 -16.82581  2.010924
>> 0.040  0.040 -18.69819  1.611304
>> 0.080  0.080 -21.05065  1.283928
>>   >
>>
>> pd<-
>> structure(list(lambda = c(0, 0.005, 0.01, 0.02, 0.04, 0.08),
>>       log.det = c(-12.9270978142337, -14.411442487768, -
>> 15.4106886674014,
>>       -16.8258120792945, -18.6981870228698, -21.050646106925),
>>       norm.beta = c(3.8068008759562, 2.81945995964196, 2.42319655878575,
>>       2.01092421747594, 1.6113040561427, 1.28392804825009)), .Names =
>> c("lambda", "log.det", "norm.beta"), class = "data.frame", row.names =
>> c("0.000", "0.005", "0.010", "0.020", "0.040", "0.080"))
>>
>> clr<- c("black", rainbow(5, start=.6, end=.1)) lambdaf<-
>> c(expression(~widehat(beta)^OLS), ".005", ".01", ".02", ".04", ".08")
>> op<- par(mar=c(4, 4, 1, 1) + 0.2, xpd=TRUE) with(pd, {plot(norm.beta,
>> log.det, type="b",
>>       cex.lab=1.25, pch=16, cex=1.5, col=clr,
>>     xlab='shrinkage: ||b||',
>>       ylab='variance: log |(Var(b)|)')
>>       text(norm.beta, log.det, lambdaf, cex=1.25, pos=2)
>>       text(min(norm.beta), max(log.det), "Variance vs. Shrinkage",
>> cex=1.5, pos=4)
>>       })
>>
>>
>> # How to find the (x,y) positions for these values of lambda along the
>> curve of log.det ~ norm.beta ?
>> lambda.HKB<- 0.004275357
>> lambda.LW<- 0.03229531
>>
>> --
>> Michael Friendly     Email: friendly AT yorku DOT ca
>> Professor, Psychology Dept.
>> York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
>> 4700 Keele Street    Web:   http://www.datavis.ca
>> Toronto, ONT  M3J 1P3 CANADA
>>
>> ______________________________________________
>> 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.


-- 
Michael Friendly     Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    Web:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA



More information about the R-help mailing list