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@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.
