[R] efficient sine interpolation

Frede Aakmann Tøgersen frtog at vestas.com
Tue May 13 08:43:57 CEST 2014


Hi

Why not use functional data analysis. There is an R package called fda to help you.

See http://www.psych.mcgill.ca/misc/fda/ex-weather-a1.html for an example exploring temperature from different locations.

Yours sincerely / Med venlig hilsen


Frede Aakmann Tøgersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender. 


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On Behalf Of Ortiz-Bobea, Ariel
> Sent: 13. maj 2014 05:42
> To: r-help at r-project.org
> Subject: [R] efficient sine interpolation
> 
> Hello,
> 
> I'm trying to fit a sine curve over successive temperature readings (i.e.
> minimum and maximum temperature) over several days and for many
> locations. The code below shows a hypothetical example of 5000 locations
> with  7 days of temperature data. Not very efficient when you have many
> more locations and days.
> 
> The linear interpolation takes 0.7 seconds, and the sine interpolations take 2
> to 4 seconds depending on the approach.
> 
> Any ideas on how to speed this up? Thanks in advance.
> 
> Ariel
> 
> ### R Code ######
> 
> # 1- Prepare data fake data
>   days<- 7
>   n <- 5000*days
>   tmin <- matrix(rnorm(n, mean=0) , ncol=days, nrow=5000)
>   tmax <- matrix(rnorm(n, mean=10), ncol=days, nrow=5000)
>   m <- matrix(NA, ncol=days*2, nrow=5000)
>   m[,seq(1,ncol(m),2)]  <- tmin
>   m[,seq(2,ncol(m)+1,2)]<- tmax
>   # check first row
>   plot(1:ncol(m), m[1,], type="l")
> 
> # 2 -linear interpolation: 0.66 seconds
>   xout <- seq(0,ncol(m),0.25/24*2)[-1] # time step = 0.25 hours or 15 minutes
>   system.time( m1 <- t(apply(m,1, function(y) approx(x=1:ncol(m), y=y,
> xout=xout, method="linear")$y)) )
>   # Check first row
>   plot(1:ncol(m), m[1,], type="l")
>   points(xout, m1[1,], col="red", cex=1)
> 
> 
> # 3- sine interpolation
>   sine.approx1 <- function(index, tmin, tmax) {
>     b <- (2*pi)/24  # period = 24 hours
>     c <- pi/2 # horizontal shift
>     xout <- seq(0,24,0.25)[-1]
>     yhat <- apply(cbind(tmin[index,],tmax[index,]), 1, function(z) diff(z)/2 *
> sin(b*xout-c) + mean(z))
>     #yhat <- yhat[-nrow(yhat),]
>     yhat <- c(yhat)
>     #plot(yhat, type="l")
>   }
>   sine.approx2 <- function(index, tmin, tmax) {
>     b <- (2*pi)/24  # period = 24 hours
>     c <- pi/2 # horizontal shift
>     xout1 <- seq(0 ,12,0.25)
>     xout2 <- seq(12,24,0.25)[-1]
>     xout2 <- xout2[-length(xout2)]
>     yhat1 <- apply(cbind(tmin[index,]                       ,tmax[index,]    ), 1,
> function(z) diff(z)/2 * sin(b*xout1-c) + mean(z))
>     yhat2 <- apply(cbind(tmax[index,][-length(tmax[index,])],tmin[index,][-
> 1]), 1, function(z) diff(z)/2 * sin(b*xout2+c) + mean(z))
>     yhat2 <- cbind(yhat2,NA)
>     yhat3 <- rbind(yhat1,yhat2)
>     #yhat3 <- yhat3[-nrow(yhat3),]
>     yhat3 <- c(yhat3)
>     yhat <- yhat3
>     #plot(c(yhat1))
>     #plot(c(yhat2))
>     #plot(yhat, type="l")
>   }
> 
>   # Single sine: 2.23 seconds
>   system.time( m2 <- t(sapply(1:nrow(m), function(i) sine.approx1(i,
> tmin=tmin, tmax=tmax))) )
> 
>   # Double sine: 4.03 seconds
>   system.time( m3 <- t(sapply(1:nrow(m), function(i) sine.approx2(i,
> tmin=tmin, tmax=tmax))) )
> 
>   # take a look at approach 1
>   plot(seq(-1,ncol(m)-1,1)[-1], m[1,], type="l")
>   points(xout, m2[1,], col="red", cex=1)
> 
>   # take a look at approach 2
>   plot(seq(-1,ncol(m)-1,1)[-1], m[1,], type="l")
>   points(xout, m3[1,], col="blue", cex=1)
> 
> 
> ---
> Ariel Ortiz-Bobea
> Fellow
> Resources for the Future
> 1616 P Street, N.W.
> Washington, DC 20036
> 202-328-5173
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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