Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are just
totally ridiculous so feel free to have a laugh but would appreciate if
someone can give me a hint. Otherwise I guess I'll have to give RCpp a
try.....
## Bench mark the recursion functions
loopRec <- function(x, alpha){
n <- length(x)
y <- double(n)
for(i in 1:n){
y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- matrix(eval(parse(text = paste("c(",
paste(paste(paste("rep(0, ", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2 <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1 <- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
up.mat2 <- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat <- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000, 0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
12 messages · Enrico Schumann, Michael Kao, michael.weylandt at gmail.com (R. Michael Weylandt +4 more
Hi Michael
here are two variations of your loop, and both seem faster than the
original loop (on my computer).
require("rbenchmark")
## your function
loopRec <- function(x, alpha){
n <- length(x)
y <- double(n)
for(i in 1:n){
y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2 <- function(x, alpha){
n <- length(x)
y <- numeric(n)
for(i in seq_len(n)){
y[i] <- sum(cumprod(rep.int(alpha, i)) * x[i:1])
}
y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3 <- function(x, alpha){
n <- length(x)
y <- numeric(n)
aa <- cumprod(rep.int(alpha, n))
for(i in seq_len(n)){
y[i] <- sum(aa[seq_len(i)] * x[i:1])
}
y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) 50 0.77 1.000000 0.76 0.00
3 loopRec3(1:1000, 0.5) 50 0.86 1.116883 0.85 0.00
1 loopRec(1:1000, 0.5) 50 1.84 2.389610 1.79 0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are just
totally ridiculous so feel free to have a laugh but would appreciate if
someone can give me a hint. Otherwise I guess I'll have to give RCpp a
try.....
## Bench mark the recursion functions
loopRec <- function(x, alpha){
n <- length(x)
y <- double(n)
for(i in 1:n){
y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- matrix(eval(parse(text = paste("c(", paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2 <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1 <- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
up.mat2 <- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat <- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000, 0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
Enrico Schumann Lucerne, Switzerland http://nmof.net/
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael
On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than the
original loop (on my computer).
require("rbenchmark")
## your function
loopRec <- function(x, alpha){
n <- length(x)
y <- double(n)
for(i in 1:n){
y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2 <- function(x, alpha){
n <- length(x)
y <- numeric(n)
for(i in seq_len(n)){
y[i] <- sum(cumprod(rep.int(alpha, i)) * x[i:1])
}
y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3 <- function(x, alpha){
n <- length(x)
y <- numeric(n)
aa <- cumprod(rep.int(alpha, n))
for(i in seq_len(n)){
y[i] <- sum(aa[seq_len(i)] * x[i:1])
}
y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) 50 0.77 1.000000 0.76 0.00
3 loopRec3(1:1000, 0.5) 50 0.86 1.116883 0.85 0.00
1 loopRec(1:1000, 0.5) 50 1.84 2.389610 1.79 0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are just
totally ridiculous so feel free to have a laugh but would appreciate if
someone can give me a hint. Otherwise I guess I'll have to give RCpp a
try.....
## Bench mark the recursion functions
loopRec <- function(x, alpha){
n <- length(x)
y <- double(n)
for(i in 1:n){
y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- matrix(eval(parse(text = paste("c(", paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2 <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1 <- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
up.mat2 <- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat <- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000,
0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec <- function(x, alpha){
n <- length(x)
y <- numeric(n)
if (n == 0) return(y)
y[1] <- alpha * x[1]
for(i in seq_len(n)[-1]){
y[i] <- alpha * (y[i-1] + x[i])
}
return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao <mkao006rmail at gmail.com> wrote:
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than the original loop (on my computer).
require("rbenchmark")
## your function
loopRec <- function(x, alpha){
? ?n <- length(x)
? ?y <- double(n)
? ?for(i in 1:n){
? ? ? ?y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
? ?}
? ?y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2 <- function(x, alpha){
? ?n <- length(x)
? ?y <- numeric(n)
? ?for(i in seq_len(n)){
? ? ? ?y[i] <- sum(cumprod(rep.int(alpha, i)) * x[i:1])
? ?}
? ?y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3 <- function(x, alpha){
? ?n <- length(x)
? ?y <- numeric(n)
? ?aa <- cumprod(rep.int(alpha, n))
? ?for(i in seq_len(n)){
? ? ? ?y[i] <- sum(aa[seq_len(i)] * x[i:1])
? ?}
? ?y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
?loopRec3(1:1000, 0.5),
?replications = 50, order = "relative")
... I get
? ? ? ? ? ? ? ? ? test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) ? ? ? ? ? 50 ? ?0.77 1.000000 ? ? ?0.76 ? ? 0.00
3 loopRec3(1:1000, 0.5) ? ? ? ? ? 50 ? ?0.86 1.116883 ? ? ?0.85 ? ? 0.00
1 ?loopRec(1:1000, 0.5) ? ? ? ? ? 50 ? ?1.84 2.389610 ? ? ?1.79 ? ? 0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are just
totally ridiculous so feel free to have a laugh but would appreciate if
someone can give me a hint. Otherwise I guess I'll have to give RCpp a
try.....
## Bench mark the recursion functions
loopRec <- function(x, alpha){
n <- length(x)
y <- double(n)
for(i in 1:n){
y[i] <- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- matrix(eval(parse(text = paste("c(", paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat <- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2 <- function(x, alpha){
n <- length(x)
exp.mat <- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1 <- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
up.mat2 <- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat <- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)] <- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000, 0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
______________________________________________ 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.
Hi Florent, That is great, I was working on solving the recurrence equation and this was part of that equation. Now I know how to put everything together, thanks for all the help e everybody! Cheers,
On 27/11/2011 2:05 p.m., Florent D. wrote:
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec<- function(x, alpha){
n<- length(x)
y<- numeric(n)
if (n == 0) return(y)
y[1]<- alpha * x[1]
for(i in seq_len(n)[-1]){
y[i]<- alpha * (y[i-1] + x[i])
}
return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao<mkao006rmail at gmail.com> wrote:
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than the original loop (on my computer).
require("rbenchmark")
## your function
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2<- function(x, alpha){
n<- length(x)
y<- numeric(n)
for(i in seq_len(n)){
y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
}
y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3<- function(x, alpha){
n<- length(x)
y<- numeric(n)
aa<- cumprod(rep.int(alpha, n))
for(i in seq_len(n)){
y[i]<- sum(aa[seq_len(i)] * x[i:1])
}
y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) 50 0.77 1.000000 0.76 0.00
3 loopRec3(1:1000, 0.5) 50 0.86 1.116883 0.85 0.00
1 loopRec(1:1000, 0.5) 50 1.84 2.389610 1.79 0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are just
totally ridiculous so feel free to have a laugh but would appreciate if
someone can give me a hint. Otherwise I guess I'll have to give RCpp a
try.....
## Bench mark the recursion functions
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- matrix(eval(parse(text = paste("c(", paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat<- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000, 0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
______________________________________________ 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.
You also might look at EMA() in the TTR package for a C implementation. (I think it matches your problem but can't promise) Michael
On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com> wrote:
Hi Florent, That is great, I was working on solving the recurrence equation and this was part of that equation. Now I know how to put everything together, thanks for all the help e everybody! Cheers, On 27/11/2011 2:05 p.m., Florent D. wrote:
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec<- function(x, alpha){
n<- length(x)
y<- numeric(n)
if (n == 0) return(y)
y[1]<- alpha * x[1]
for(i in seq_len(n)[-1]){
y[i]<- alpha * (y[i-1] + x[i])
}
return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao<mkao006rmail at gmail.com> wrote:
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than the original loop (on my computer).
require("rbenchmark")
## your function
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2<- function(x, alpha){
n<- length(x)
y<- numeric(n)
for(i in seq_len(n)){
y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
}
y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3<- function(x, alpha){
n<- length(x)
y<- numeric(n)
aa<- cumprod(rep.int(alpha, n))
for(i in seq_len(n)){
y[i]<- sum(aa[seq_len(i)] * x[i:1])
}
y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) 50 0.77 1.000000 0.76 0.00
3 loopRec3(1:1000, 0.5) 50 0.86 1.116883 0.85 0.00
1 loopRec(1:1000, 0.5) 50 1.84 2.389610 1.79 0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are just
totally ridiculous so feel free to have a laugh but would appreciate if
someone can give me a hint. Otherwise I guess I'll have to give RCpp a
try.....
## Bench mark the recursion functions
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- matrix(eval(parse(text = paste("c(", paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow = TRUE)
up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat<- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000, 0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
You also might look at EMA() in the TTR package for a C implementation. (I think it matches your problem but can't promise)
It's pretty close for EMA(1:1000, n=1, ratio=0.5) and 7 times faster. > EMA(1:10, n=1, ratio=0.5) [1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625 7.007812 8.003906 [10] 9.001953 > loopRec3(1:10, 0.5) [1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812 7.003906 8.001953 [10] 9.000977 This series converges very quickly to n-1 when the ratio = 0.5 and to n-3 when ratio = 0.25. It's already within 1% at the fifth iteration. There may be a simple analytical approximation that makes the process even more efficient. The asymptotic result appears to be: n-(1-ratio)/ ratio > tail(EMA(1:1000, n=1, ratio=0.5)) [1] 994 995 996 997 998 999 > tail(loopRec3(1:1000, 0.5)) [1] 994 995 996 997 998 999 > EMA(1:1000, n=1, ratio=0.1)[491:500] [1] 482 483 484 485 486 487 488 489 490 491
David.
> Michael
>
> On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com>
> wrote:
>
>> Hi Florent,
>>
>> That is great, I was working on solving the recurrence equation and
>> this was part of that equation. Now I know how to put everything
>> together, thanks for all the help e everybody!
>>
>> Cheers,
>>
>> On 27/11/2011 2:05 p.m., Florent D. wrote:
>>> You can make things a lot faster by using the recurrence equation
>>>
>>> y[n] = alpha * (y[n-1]+x[n])
>>>
>>> loopRec<- function(x, alpha){
>>> n<- length(x)
>>> y<- numeric(n)
>>> if (n == 0) return(y)
>>> y[1]<- alpha * x[1]
>>> for(i in seq_len(n)[-1]){
>>> y[i]<- alpha * (y[i-1] + x[i])
>>> }
>>> return(y)
>>> }
>>>
>>>
>>>
>>> On Sun, Nov 27, 2011 at 5:17 AM, Michael
>>> Kao<mkao006rmail at gmail.com> wrote:
>>>> Dear Enrico,
>>>>
>>>> Brilliant! Thank you for the improvements, not sure what I was
>>>> thinking using rev. I will test it out to see whether it is fast
>>>> enough for our implementation, but your help has been
>>>> SIGNIFICANT!!!
>>>>
>>>> Thanks heaps!
>>>> Michael
>>>>
>>>> On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
>>>>> Hi Michael
>>>>>
>>>>> here are two variations of your loop, and both seem faster than
>>>>> the original loop (on my computer).
>>>>>
>>>>>
>>>>> require("rbenchmark")
>>>>>
>>>>> ## your function
>>>>> loopRec<- function(x, alpha){
>>>>> n<- length(x)
>>>>> y<- double(n)
>>>>> for(i in 1:n){
>>>>> y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
>>>>> }
>>>>> y
>>>>> }
>>>>> loopRec(c(1, 2, 3), 0.5)
>>>>>
>>>>> loopRec2<- function(x, alpha){
>>>>> n<- length(x)
>>>>> y<- numeric(n)
>>>>> for(i in seq_len(n)){
>>>>> y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
>>>>> }
>>>>> y
>>>>> }
>>>>> loopRec2(c(1, 2, 3), 0.5)
>>>>>
>>>>> loopRec3<- function(x, alpha){
>>>>> n<- length(x)
>>>>> y<- numeric(n)
>>>>> aa<- cumprod(rep.int(alpha, n))
>>>>> for(i in seq_len(n)){
>>>>> y[i]<- sum(aa[seq_len(i)] * x[i:1])
>>>>> }
>>>>> y
>>>>> }
>>>>> loopRec3(c(1, 2, 3), 0.5)
>>>>>
>>>>>
>>>>> ## Check whether value is correct
>>>>> all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
>>>>> all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
>>>>>
>>>>>
>>>>> ## benchmark the functions.
>>>>> benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
>>>>> loopRec3(1:1000, 0.5),
>>>>> replications = 50, order = "relative")
>>>>>
>>>>>
>>>>> ... I get
>>>>> test replications elapsed relative user.self
>>>>> sys.self
>>>>> 2 loopRec2(1:1000, 0.5) 50 0.77 1.000000
>>>>> 0.76 0.00
>>>>> 3 loopRec3(1:1000, 0.5) 50 0.86 1.116883
>>>>> 0.85 0.00
>>>>> 1 loopRec(1:1000, 0.5) 50 1.84 2.389610
>>>>> 1.79 0.01
>>>>>
>>>>>
>>>>> Regards,
>>>>> Enrico
>>>>>
>>>>>
>>>>> Am 27.11.2011 01:20, schrieb Michael Kao:
>>>>>> Dear R-help,
>>>>>>
>>>>>> I have been trying really hard to generate the following vector
>>>>>> given
>>>>>> the data (x) and parameter (alpha) efficiently.
>>>>>>
>>>>>> Let y be the output list, the aim is to produce the the following
>>>>>> vector(y) with at least half the time used by the loop example
>>>>>> below.
>>>>>>
>>>>>> y[1] = alpha * x[1]
>>>>>> y[2] = alpha^2 * x[1] + alpha * x[2]
>>>>>> y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
>>>>>> .....
>>>>>>
>>>>>> below are the methods I have tried and failed miserably, some
>>>>>> are just
>>>>>> totally ridiculous so feel free to have a laugh but would
>>>>>> appreciate if
>>>>>> someone can give me a hint. Otherwise I guess I'll have to give
>>>>>> RCpp a
>>>>>> try.....
>>>>>>
>>>>>>
>>>>>> ## Bench mark the recursion functions
>>>>>> loopRec<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> y<- double(n)
>>>>>> for(i in 1:n){
>>>>>> y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
>>>>>> }
>>>>>> y
>>>>>> }
>>>>>>
>>>>>> loopRec(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> ## This is a crazy solution, but worth giving it a try.
>>>>>> charRec<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>>>>>> up.mat<- matrix(eval(parse(text = paste("c(",
>>>>>> paste(paste(paste("rep(0,
>>>>>> ", 0:(n - 1), ")", sep = ""),
>>>>>> paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep =
>>>>>> ","),
>>>>>> collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
>>>>>> colSums(up.mat * exp.mat)
>>>>>> }
>>>>>> vecRec(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> ## Sweep is slow, shouldn't use it.
>>>>>> matRec<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>>>>>> up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
>>>>>> byrow = TRUE), 1,
>>>>>> c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
>>>>>> up.mat[lower.tri(up.mat)]<- 0
>>>>>> colSums(up.mat * exp.mat)
>>>>>> }
>>>>>> matRec(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> matRec2<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>>>>>> up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow
>>>>>> = TRUE)
>>>>>> up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr
>>>>>> = n)
>>>>>> up.mat<- up.mat1 * up.mat2
>>>>>> up.mat[lower.tri(up.mat)]<- 0
>>>>>> colSums(up.mat * exp.mat)
>>>>>> }
>>>>>>
>>>>>> matRec2(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> ## Check whether value is correct
>>>>>> all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
>>>>>> all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
>>>>>> all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
>>>>>>
>>>>>> ## benchmark the functions.
>>>>>> benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5),
>>>>>> matRec(1:1000, 0.5),
>>>>>> matRec2(1:1000, 0.5), replications = 50,
>>>>>> order = "relative")
>>>>>>
>>>>>> Thank you very much for your help.
>>>>>>
>>>>>> ______________________________________________
>>>>>> 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.
>>>>>>
>>>> ______________________________________________
>>>> 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.
>>
>> ______________________________________________
>> 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.
>
> ______________________________________________
> 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.
David Winsemius, MD
West Hartford, CT
Hi Is there a way to link a frequently-used sub-routine (in binary form or preferably in ascii form) in a program, similar to the following in Gauss? Thanks. #include c:\programs\mylib1.g; -- Steven T. Yen, Professor of Agricultural Economics The University of Tennessee http://web.utk.edu/~syen/
From: r-help-bounces at r-project.org [r-help-bounces at r-project.org] on behalf of David Winsemius [dwinsemius at comcast.net]
Sent: Sunday, November 27, 2011 11:46 AM
To: R. Michael Weylandt <michael.weylandt at gmail.com>
Cc: r-help at r-project.org
Subject: Re: [R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
Sent: Sunday, November 27, 2011 11:46 AM
To: R. Michael Weylandt <michael.weylandt at gmail.com>
Cc: r-help at r-project.org
Subject: Re: [R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
> You also might look at EMA() in the TTR package for a C
> implementation. (I think it matches your problem but can't promise)
>
It's pretty close for EMA(1:1000, n=1, ratio=0.5) and 7 times faster.
> EMA(1:10, n=1, ratio=0.5)
[1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625
7.007812 8.003906
[10] 9.001953
> loopRec3(1:10, 0.5)
[1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812
7.003906 8.001953
[10] 9.000977
This series converges very quickly to n-1 when the ratio = 0.5 and to
n-3 when ratio = 0.25. It's already within 1% at the fifth iteration.
There may be a simple analytical approximation that makes the process
even more efficient. The asymptotic result appears to be: n-(1-ratio)/
ratio
> tail(EMA(1:1000, n=1, ratio=0.5))
[1] 994 995 996 997 998 999
> tail(loopRec3(1:1000, 0.5))
[1] 994 995 996 997 998 999
> EMA(1:1000, n=1, ratio=0.1)[491:500]
[1] 482 483 484 485 486 487 488 489 490 491
--
David.
> Michael
>
> On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com>
> wrote:
>
>> Hi Florent,
>>
>> That is great, I was working on solving the recurrence equation and
>> this was part of that equation. Now I know how to put everything
>> together, thanks for all the help e everybody!
>>
>> Cheers,
>>
>> On 27/11/2011 2:05 p.m., Florent D. wrote:
>>> You can make things a lot faster by using the recurrence equation
>>>
>>> y[n] = alpha * (y[n-1]+x[n])
>>>
>>> loopRec<- function(x, alpha){
>>> n<- length(x)
>>> y<- numeric(n)
>>> if (n == 0) return(y)
>>> y[1]<- alpha * x[1]
>>> for(i in seq_len(n)[-1]){
>>> y[i]<- alpha * (y[i-1] + x[i])
>>> }
>>> return(y)
>>> }
>>>
>>>
>>>
>>> On Sun, Nov 27, 2011 at 5:17 AM, Michael
>>> Kao<mkao006rmail at gmail.com> wrote:
>>>> Dear Enrico,
>>>>
>>>> Brilliant! Thank you for the improvements, not sure what I was
>>>> thinking using rev. I will test it out to see whether it is fast
>>>> enough for our implementation, but your help has been
>>>> SIGNIFICANT!!!
>>>>
>>>> Thanks heaps!
>>>> Michael
>>>>
>>>> On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
>>>>> Hi Michael
>>>>>
>>>>> here are two variations of your loop, and both seem faster than
>>>>> the original loop (on my computer).
>>>>>
>>>>>
>>>>> require("rbenchmark")
>>>>>
>>>>> ## your function
>>>>> loopRec<- function(x, alpha){
>>>>> n<- length(x)
>>>>> y<- double(n)
>>>>> for(i in 1:n){
>>>>> y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
>>>>> }
>>>>> y
>>>>> }
>>>>> loopRec(c(1, 2, 3), 0.5)
>>>>>
>>>>> loopRec2<- function(x, alpha){
>>>>> n<- length(x)
>>>>> y<- numeric(n)
>>>>> for(i in seq_len(n)){
>>>>> y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
>>>>> }
>>>>> y
>>>>> }
>>>>> loopRec2(c(1, 2, 3), 0.5)
>>>>>
>>>>> loopRec3<- function(x, alpha){
>>>>> n<- length(x)
>>>>> y<- numeric(n)
>>>>> aa<- cumprod(rep.int(alpha, n))
>>>>> for(i in seq_len(n)){
>>>>> y[i]<- sum(aa[seq_len(i)] * x[i:1])
>>>>> }
>>>>> y
>>>>> }
>>>>> loopRec3(c(1, 2, 3), 0.5)
>>>>>
>>>>>
>>>>> ## Check whether value is correct
>>>>> all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
>>>>> all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
>>>>>
>>>>>
>>>>> ## benchmark the functions.
>>>>> benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
>>>>> loopRec3(1:1000, 0.5),
>>>>> replications = 50, order = "relative")
>>>>>
>>>>>
>>>>> ... I get
>>>>> test replications elapsed relative user.self
>>>>> sys.self
>>>>> 2 loopRec2(1:1000, 0.5) 50 0.77 1.000000
>>>>> 0.76 0.00
>>>>> 3 loopRec3(1:1000, 0.5) 50 0.86 1.116883
>>>>> 0.85 0.00
>>>>> 1 loopRec(1:1000, 0.5) 50 1.84 2.389610
>>>>> 1.79 0.01
>>>>>
>>>>>
>>>>> Regards,
>>>>> Enrico
>>>>>
>>>>>
>>>>> Am 27.11.2011 01:20, schrieb Michael Kao:
>>>>>> Dear R-help,
>>>>>>
>>>>>> I have been trying really hard to generate the following vector
>>>>>> given
>>>>>> the data (x) and parameter (alpha) efficiently.
>>>>>>
>>>>>> Let y be the output list, the aim is to produce the the following
>>>>>> vector(y) with at least half the time used by the loop example
>>>>>> below.
>>>>>>
>>>>>> y[1] = alpha * x[1]
>>>>>> y[2] = alpha^2 * x[1] + alpha * x[2]
>>>>>> y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
>>>>>> .....
>>>>>>
>>>>>> below are the methods I have tried and failed miserably, some
>>>>>> are just
>>>>>> totally ridiculous so feel free to have a laugh but would
>>>>>> appreciate if
>>>>>> someone can give me a hint. Otherwise I guess I'll have to give
>>>>>> RCpp a
>>>>>> try.....
>>>>>>
>>>>>>
>>>>>> ## Bench mark the recursion functions
>>>>>> loopRec<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> y<- double(n)
>>>>>> for(i in 1:n){
>>>>>> y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
>>>>>> }
>>>>>> y
>>>>>> }
>>>>>>
>>>>>> loopRec(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> ## This is a crazy solution, but worth giving it a try.
>>>>>> charRec<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>>>>>> up.mat<- matrix(eval(parse(text = paste("c(",
>>>>>> paste(paste(paste("rep(0,
>>>>>> ", 0:(n - 1), ")", sep = ""),
>>>>>> paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep =
>>>>>> ","),
>>>>>> collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
>>>>>> colSums(up.mat * exp.mat)
>>>>>> }
>>>>>> vecRec(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> ## Sweep is slow, shouldn't use it.
>>>>>> matRec<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>>>>>> up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
>>>>>> byrow = TRUE), 1,
>>>>>> c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
>>>>>> up.mat[lower.tri(up.mat)]<- 0
>>>>>> colSums(up.mat * exp.mat)
>>>>>> }
>>>>>> matRec(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> matRec2<- function(x, alpha){
>>>>>> n<- length(x)
>>>>>> exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
>>>>>> up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow
>>>>>> = TRUE)
>>>>>> up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr
>>>>>> = n)
>>>>>> up.mat<- up.mat1 * up.mat2
>>>>>> up.mat[lower.tri(up.mat)]<- 0
>>>>>> colSums(up.mat * exp.mat)
>>>>>> }
>>>>>>
>>>>>> matRec2(c(1, 2, 3), 0.5)
>>>>>>
>>>>>> ## Check whether value is correct
>>>>>> all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
>>>>>> all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
>>>>>> all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
>>>>>>
>>>>>> ## benchmark the functions.
>>>>>> benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5),
>>>>>> matRec(1:1000, 0.5),
>>>>>> matRec2(1:1000, 0.5), replications = 50,
>>>>>> order = "relative")
>>>>>>
>>>>>> Thank you very much for your help.
>>>>>>
>>>>>> ______________________________________________
>>>>>> 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.
>>>>>>
>>>> ______________________________________________
>>>> 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.
>>
>> ______________________________________________
>> 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.
>
> ______________________________________________
> 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.
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
... but the original request was to build a series, not approximate its limit by building a different (but "faster") series. What I suggested is linear in terms of the size of x so you won't find a faster algorithm. What will make a difference is the implementation, e.g. C loop versus R loop. On Sun, Nov 27, 2011 at 11:46 AM, David Winsemius
<dwinsemius at comcast.net> wrote:
On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
You also might look at EMA() in the TTR package for a C implementation. (I think it matches your problem but can't promise)
It's pretty close for EMA(1:1000, n=1, ?ratio=0.5) and 7 times faster.
EMA(1:10, n=1, ?ratio=0.5)
?[1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625 7.007812 8.003906 [10] 9.001953
loopRec3(1:10, 0.5)
?[1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812 7.003906 8.001953 [10] 9.000977 This series converges very quickly to n-1 when the ratio = 0.5 and to n-3 when ratio = 0.25. It's already within 1% at the fifth iteration. There may be a simple analytical approximation that makes the process even more efficient. The asymptotic result appears to be: ?n-(1-ratio)/ratio
tail(EMA(1:1000, n=1, ?ratio=0.5))
[1] 994 995 996 997 998 999
tail(loopRec3(1:1000, 0.5))
[1] 994 995 996 997 998 999
EMA(1:1000, n=1, ?ratio=0.1)[491:500]
?[1] 482 483 484 485 486 487 488 489 490 491 -- David.
Michael On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com> wrote:
Hi Florent, That is great, I was working on solving the recurrence equation and this was part of that equation. Now I know how to put everything together, thanks for all the help e everybody! Cheers, On 27/11/2011 2:05 p.m., Florent D. wrote:
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec<- function(x, alpha){
?n<- length(x)
?y<- numeric(n)
?if (n == 0) return(y)
?y[1]<- alpha * x[1]
?for(i in seq_len(n)[-1]){
? ? y[i]<- alpha * (y[i-1] + x[i])
?}
?return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao<mkao006rmail at gmail.com>
?wrote:
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than the
original loop (on my computer).
require("rbenchmark")
## your function
loopRec<- function(x, alpha){
?n<- length(x)
?y<- double(n)
?for(i in 1:n){
? ? ?y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
?}
?y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2<- function(x, alpha){
?n<- length(x)
?y<- numeric(n)
?for(i in seq_len(n)){
? ? ?y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
?}
?y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3<- function(x, alpha){
?n<- length(x)
?y<- numeric(n)
?aa<- cumprod(rep.int(alpha, n))
?for(i in seq_len(n)){
? ? ?y[i]<- sum(aa[seq_len(i)] * x[i:1])
?}
?y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
? ? ? ? ? ? ? ? test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) ? ? ? ? ? 50 ? ?0.77 1.000000 ? ? ?0.76
0.00
3 loopRec3(1:1000, 0.5) ? ? ? ? ? 50 ? ?0.86 1.116883 ? ? ?0.85
0.00
1 ?loopRec(1:1000, 0.5) ? ? ? ? ? 50 ? ?1.84 2.389610 ? ? ?1.79
0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are
just
totally ridiculous so feel free to have a laugh but would appreciate
if
someone can give me a hint. Otherwise I guess I'll have to give RCpp
a
try.....
## Bench mark the recursion functions
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- matrix(eval(parse(text = paste("c(",
paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow =
TRUE)
up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat<- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000,
0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
David Winsemius, MD West Hartford, CT
______________________________________________ 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.
The following quickly does the recursion suggested
by Florent D.:
loopRec5 <- function(x, alpha) {
as.vector(filter(x*alpha, alpha, method="recursive"))
}
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Florent D.
Sent: Sunday, November 27, 2011 10:01 AM
To: David Winsemius
Cc: r-help at r-project.org
Subject: Re: [R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
... but the original request was to build a series, not approximate
its limit by building a different (but "faster") series.
What I suggested is linear in terms of the size of x so you won't find
a faster algorithm. What will make a difference is the implementation,
e.g. C loop versus R loop.
On Sun, Nov 27, 2011 at 11:46 AM, David Winsemius
<dwinsemius at comcast.net> wrote:
On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
You also might look at EMA() in the TTR package for a C implementation. (I think it matches your problem but can't promise)
It's pretty close for EMA(1:1000, n=1, ?ratio=0.5) and 7 times faster.
EMA(1:10, n=1, ?ratio=0.5)
?[1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625 7.007812 8.003906 [10] 9.001953
loopRec3(1:10, 0.5)
?[1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812 7.003906 8.001953 [10] 9.000977 This series converges very quickly to n-1 when the ratio = 0.5 and to n-3 when ratio = 0.25. It's already within 1% at the fifth iteration. There may be a simple analytical approximation that makes the process even more efficient. The asymptotic result appears to be: ?n-(1-ratio)/ratio
tail(EMA(1:1000, n=1, ?ratio=0.5))
[1] 994 995 996 997 998 999
tail(loopRec3(1:1000, 0.5))
[1] 994 995 996 997 998 999
EMA(1:1000, n=1, ?ratio=0.1)[491:500]
?[1] 482 483 484 485 486 487 488 489 490 491 -- David.
Michael On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com> wrote:
Hi Florent, That is great, I was working on solving the recurrence equation and this was part of that equation. Now I know how to put everything together, thanks for all the help e everybody! Cheers, On 27/11/2011 2:05 p.m., Florent D. wrote:
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec<- function(x, alpha){
?n<- length(x)
?y<- numeric(n)
?if (n == 0) return(y)
?y[1]<- alpha * x[1]
?for(i in seq_len(n)[-1]){
? ? y[i]<- alpha * (y[i-1] + x[i])
?}
?return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao<mkao006rmail at gmail.com>
?wrote:
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than the
original loop (on my computer).
require("rbenchmark")
## your function
loopRec<- function(x, alpha){
?n<- length(x)
?y<- double(n)
?for(i in 1:n){
? ? ?y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
?}
?y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2<- function(x, alpha){
?n<- length(x)
?y<- numeric(n)
?for(i in seq_len(n)){
? ? ?y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
?}
?y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3<- function(x, alpha){
?n<- length(x)
?y<- numeric(n)
?aa<- cumprod(rep.int(alpha, n))
?for(i in seq_len(n)){
? ? ?y[i]<- sum(aa[seq_len(i)] * x[i:1])
?}
?y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
? ? ? ? ? ? ? ? test replications elapsed relative user.self sys.self
2 loopRec2(1:1000, 0.5) ? ? ? ? ? 50 ? ?0.77 1.000000 ? ? ?0.76
0.00
3 loopRec3(1:1000, 0.5) ? ? ? ? ? 50 ? ?0.86 1.116883 ? ? ?0.85
0.00
1 ?loopRec(1:1000, 0.5) ? ? ? ? ? 50 ? ?1.84 2.389610 ? ? ?1.79
0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some are
just
totally ridiculous so feel free to have a laugh but would appreciate
if
someone can give me a hint. Otherwise I guess I'll have to give RCpp
a
try.....
## Bench mark the recursion functions
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- matrix(eval(parse(text = paste("c(",
paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep = ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow =
TRUE)
up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr = n)
up.mat<- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5), matRec(1:1000,
0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
David Winsemius, MD West Hartford, CT
______________________________________________ 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.
______________________________________________ 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.
You can write a function to do so and then there are all sorts of options of what to do with it - mostly they degenerate into putting it in a personal package or into your .Rprofile for convenient use. If you want to use compiled code, that's slightly trickier but good facilities exist to do so from C, C++, and Fortran. Any follow-up, after reading the R extensions manual and the documentation to ?package.skeleton, ?STARTUP, and ?function, probably deserves its own thread. Michael
On Nov 27, 2011, at 12:34 PM, "Yen, Steven T" <syen at utk.edu> wrote:
Hi Is there a way to link a frequently-used sub-routine (in binary form or preferably in ascii form) in a program, similar to the following in Gauss? Thanks. #include c:\programs\mylib1.g; -- Steven T. Yen, Professor of Agricultural Economics The University of Tennessee http://web.utk.edu/~syen/
________________________________________
From: r-help-bounces at r-project.org [r-help-bounces at r-project.org] on behalf of David Winsemius [dwinsemius at comcast.net]
Sent: Sunday, November 27, 2011 11:46 AM
To: R. Michael Weylandt <michael.weylandt at gmail.com>
Cc: r-help at r-project.org
Subject: Re: [R] generating a vector of y_t = \sum_{i = 1}^t (alpha^i * x_{t - i + 1})
On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
You also might look at EMA() in the TTR package for a C
implementation. (I think it matches your problem but can't promise)
It's pretty close for EMA(1:1000, n=1, ratio=0.5) and 7 times faster.
EMA(1:10, n=1, ratio=0.5)
[1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625
7.007812 8.003906
[10] 9.001953
loopRec3(1:10, 0.5)
[1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812
7.003906 8.001953
[10] 9.000977
This series converges very quickly to n-1 when the ratio = 0.5 and to
n-3 when ratio = 0.25. It's already within 1% at the fifth iteration.
There may be a simple analytical approximation that makes the process
even more efficient. The asymptotic result appears to be: n-(1-ratio)/
ratio
tail(EMA(1:1000, n=1, ratio=0.5))
[1] 994 995 996 997 998 999
tail(loopRec3(1:1000, 0.5))
[1] 994 995 996 997 998 999
EMA(1:1000, n=1, ratio=0.1)[491:500]
[1] 482 483 484 485 486 487 488 489 490 491
--
David.
Michael
On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com>
wrote:
Hi Florent,
That is great, I was working on solving the recurrence equation and
this was part of that equation. Now I know how to put everything
together, thanks for all the help e everybody!
Cheers,
On 27/11/2011 2:05 p.m., Florent D. wrote:
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec<- function(x, alpha){
n<- length(x)
y<- numeric(n)
if (n == 0) return(y)
y[1]<- alpha * x[1]
for(i in seq_len(n)[-1]){
y[i]<- alpha * (y[i-1] + x[i])
}
return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael
Kao<mkao006rmail at gmail.com> wrote:
Dear Enrico,
Brilliant! Thank you for the improvements, not sure what I was
thinking using rev. I will test it out to see whether it is fast
enough for our implementation, but your help has been
SIGNIFICANT!!!
Thanks heaps!
Michael
On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster than
the original loop (on my computer).
require("rbenchmark")
## your function
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2<- function(x, alpha){
n<- length(x)
y<- numeric(n)
for(i in seq_len(n)){
y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
}
y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3<- function(x, alpha){
n<- length(x)
y<- numeric(n)
aa<- cumprod(rep.int(alpha, n))
for(i in seq_len(n)){
y[i]<- sum(aa[seq_len(i)] * x[i:1])
}
y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
test replications elapsed relative user.self
sys.self
2 loopRec2(1:1000, 0.5) 50 0.77 1.000000
0.76 0.00
3 loopRec3(1:1000, 0.5) 50 0.86 1.116883
0.85 0.00
1 loopRec(1:1000, 0.5) 50 1.84 2.389610
1.79 0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following vector
given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the following
vector(y) with at least half the time used by the loop example
below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some
are just
totally ridiculous so feel free to have a laugh but would
appreciate if
someone can give me a hint. Otherwise I guess I'll have to give
RCpp a
try.....
## Bench mark the recursion functions
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- matrix(eval(parse(text = paste("c(",
paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep =
","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n, byrow
= TRUE)
up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n, nr
= n)
up.mat<- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5),
matRec(1:1000, 0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________
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.
______________________________________________
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.
______________________________________________
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.
______________________________________________
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.
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
On Nov 27, 2011, at 1:00 PM, Florent D. wrote:
... but the original request was to build a series, not approximate its limit by building a different (but "faster") series.
Right. Your function was faster that the earlier ones. But if speed were the issue, it might make more sense to use mathematics.
What I suggested is linear in terms of the size of x so you won't find a faster algorithm.
Actually, Dunlap's was 13-fold faster than yours on my machine.
> loopRec5 <- function(x, alpha) {
+ as.vector(filter(x*alpha, alpha, method="recursive"))
+ }
>
> loopRec4 <- function(x, alpha){
+ n <- length(x)
+ y <- numeric(n)
+ if (n == 0) return(y)
+ y[1] <- alpha * x[1]
+ for(i in seq_len(n)[-1]){
+ y[i] <- alpha * (y[i-1] + x[i])
+ }
+ return(y)
+ }
>
>
> ## benchmark the functions.
> benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
+ loopRec3(1:1000, 0.5),loopRec4(1:1000, 0.5),loopRec5(1:1000, 0.5),
+ replications = 50, order = "relative")
test replications elapsed relative user.self
sys.self user.child
5 loopRec5(1:1000, 0.5) 50 0.017 1.00000 0.016
0.000 0
4 loopRec4(1:1000, 0.5) 50 0.228 13.41176 0.227
0.000 0
3 loopRec3(1:1000, 0.5) 50 1.458 85.76471 1.461
0.005 0
2 loopRec2(1:1000, 0.5) 50 1.579 92.88235 1.579
0.006 0
1 loopRec(1:1000, 0.5) 50 3.006 176.82353 3.006
0.012 0
What will make a difference is the implementation, e.g. C loop versus R loop. On Sun, Nov 27, 2011 at 11:46 AM, David Winsemius <dwinsemius at comcast.net> wrote:
On Nov 27, 2011, at 9:25 AM, R. Michael Weylandt wrote:
You also might look at EMA() in the TTR package for a C implementation. (I think it matches your problem but can't promise)
It's pretty close for EMA(1:1000, n=1, ratio=0.5) and 7 times faster.
EMA(1:10, n=1, ratio=0.5)
[1] 1.000000 1.500000 2.250000 3.125000 4.062500 5.031250 6.015625 7.007812 8.003906 [10] 9.001953
loopRec3(1:10, 0.5)
[1] 0.500000 1.250000 2.125000 3.062500 4.031250 5.015625 6.007812 7.003906 8.001953 [10] 9.000977 This series converges very quickly to n-1 when the ratio = 0.5 and to n-3 when ratio = 0.25. It's already within 1% at the fifth iteration. There may be a simple analytical approximation that makes the process even more efficient. The asymptotic result appears to be: n-(1-ratio)/ratio
tail(EMA(1:1000, n=1, ratio=0.5))
[1] 994 995 996 997 998 999
tail(loopRec3(1:1000, 0.5))
[1] 994 995 996 997 998 999
EMA(1:1000, n=1, ratio=0.1)[491:500]
[1] 482 483 484 485 486 487 488 489 490 491 -- David.
Michael On Nov 27, 2011, at 9:05 AM, Michael Kao <mkao006rmail at gmail.com> wrote:
Hi Florent, That is great, I was working on solving the recurrence equation and this was part of that equation. Now I know how to put everything together, thanks for all the help e everybody! Cheers, On 27/11/2011 2:05 p.m., Florent D. wrote:
You can make things a lot faster by using the recurrence equation
y[n] = alpha * (y[n-1]+x[n])
loopRec<- function(x, alpha){
n<- length(x)
y<- numeric(n)
if (n == 0) return(y)
y[1]<- alpha * x[1]
for(i in seq_len(n)[-1]){
y[i]<- alpha * (y[i-1] + x[i])
}
return(y)
}
On Sun, Nov 27, 2011 at 5:17 AM, Michael Kao<mkao006rmail at gmail.com
wrote:
Dear Enrico, Brilliant! Thank you for the improvements, not sure what I was thinking using rev. I will test it out to see whether it is fast enough for our implementation, but your help has been SIGNIFICANT!!! Thanks heaps! Michael On 27/11/2011 10:43 a.m., Enrico Schumann wrote:
Hi Michael
here are two variations of your loop, and both seem faster
than the
original loop (on my computer).
require("rbenchmark")
## your function
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
loopRec2<- function(x, alpha){
n<- length(x)
y<- numeric(n)
for(i in seq_len(n)){
y[i]<- sum(cumprod(rep.int(alpha, i)) * x[i:1])
}
y
}
loopRec2(c(1, 2, 3), 0.5)
loopRec3<- function(x, alpha){
n<- length(x)
y<- numeric(n)
aa<- cumprod(rep.int(alpha, n))
for(i in seq_len(n)){
y[i]<- sum(aa[seq_len(i)] * x[i:1])
}
y
}
loopRec3(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), loopRec3(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), loopRec2(1:1000, 0.5),
loopRec3(1:1000, 0.5),
replications = 50, order = "relative")
... I get
test replications elapsed relative user.self
sys.self
2 loopRec2(1:1000, 0.5) 50 0.77 1.000000 0.76
0.00
3 loopRec3(1:1000, 0.5) 50 0.86 1.116883 0.85
0.00
1 loopRec(1:1000, 0.5) 50 1.84 2.389610 1.79
0.01
Regards,
Enrico
Am 27.11.2011 01:20, schrieb Michael Kao:
Dear R-help,
I have been trying really hard to generate the following
vector given
the data (x) and parameter (alpha) efficiently.
Let y be the output list, the aim is to produce the the
following
vector(y) with at least half the time used by the loop
example below.
y[1] = alpha * x[1]
y[2] = alpha^2 * x[1] + alpha * x[2]
y[3] = alpha^3 * x[1] + alpha^2 * x[2] + alpha * x[3]
.....
below are the methods I have tried and failed miserably, some
are
just
totally ridiculous so feel free to have a laugh but would
appreciate
if
someone can give me a hint. Otherwise I guess I'll have to
give RCpp
a
try.....
## Bench mark the recursion functions
loopRec<- function(x, alpha){
n<- length(x)
y<- double(n)
for(i in 1:n){
y[i]<- sum(cumprod(rep(alpha, i)) * rev(x[1:i]))
}
y
}
loopRec(c(1, 2, 3), 0.5)
## This is a crazy solution, but worth giving it a try.
charRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- matrix(eval(parse(text = paste("c(",
paste(paste(paste("rep(0,
", 0:(n - 1), ")", sep = ""),
paste("cumprod(rep(", alpha, ",", n:1, "))", sep = "") , sep
= ","),
collapse = ","), ")", sep = ""))), nc = n, byrow = TRUE)
colSums(up.mat * exp.mat)
}
vecRec(c(1, 2, 3), 0.5)
## Sweep is slow, shouldn't use it.
matRec<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat<- sweep(matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow = TRUE), 1,
c(1, cumprod(rep(1/alpha, n - 1))), FUN = "*")
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec(c(1, 2, 3), 0.5)
matRec2<- function(x, alpha){
n<- length(x)
exp.mat<- matrix(rep(x, each = n), nc = n, byrow = TRUE)
up.mat1<- matrix(cumprod(rep(alpha, n)), nc = n, nr = n,
byrow =
TRUE)
up.mat2<- matrix(c(1, cumprod(rep(1/alpha, n - 1))), nc = n,
nr = n)
up.mat<- up.mat1 * up.mat2
up.mat[lower.tri(up.mat)]<- 0
colSums(up.mat * exp.mat)
}
matRec2(c(1, 2, 3), 0.5)
## Check whether value is correct
all.equal(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec(1:1000, 0.5))
all.equal(loopRec(1:1000, 0.5), matRec2(1:1000, 0.5))
## benchmark the functions.
benchmark(loopRec(1:1000, 0.5), vecRec(1:1000, 0.5),
matRec(1:1000,
0.5),
matRec2(1:1000, 0.5), replications = 50,
order = "relative")
Thank you very much for your help.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
______________________________________________ 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.
David Winsemius, MD West Hartford, CT
______________________________________________ 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.
David Winsemius, MD West Hartford, CT