如何在循环和if语句中使用递归公式
How to use recursive formula in loop and if statement
我可能把这个问题复杂化了。
问题是这样的
变量
a = lower triangular matrix/dataframe 20x20
b = a 1x20 matrix/vector
c = the previous row result of the formula (recursive bit)
我想要一个下三角矩阵 (I,j) 递归定义
伪代码
if(I<j) = 0
else if (I==j) = 1
else (I>j) = sum(a * b * c )/ b[I] - b[J]（其中 I 是当前行位置，J 是当前列位置）
我希望公式/R 的矢量化如何在较小的矩阵上工作以使生活更轻松的示例。
我=行
j = 列
b(我,j)
b(1,2) 指的是矩阵中的位置
示例矩阵的伪代码
Matrix C
1 | 2 |3 |4 |5 |6 |
1|i=j=1 |i<j=0 |...|...|...|...|
2|i>j=a(2,1)*b(1)*c(1,1)/b(2)-b(1) |i=j=1 |...|...|...|...|
3|i>j= (a(3,1)*b(1)*c(1,1) + a(3,2)*b(2)*c(2,1))/(b(3)-b(1) |i>j= a(3,2)*b(2)*c(2,2)/(b(3)-b(2) |...|...|...|...|
4|i>j = a(4,1)*b(1)*c(1,1)+a(4,2)*b(2)*c(2,1)+a(4,3)*b(3)*c(3,1))/b(4)-b(1) |i>j= a(4,2)*b(2)*c(2,2) + a(4,3)*b(3)*c(3,2)/b(4)-b(2) |...|...|...|...|
5|... |... |...|...|...|...|
6|... |... |...|...|...|...|
关于我目前的代码，
先创建变量如下图：
my_int <- 20
nr <- as.integer(my_int)
#create a n x n matrix with zeroes
a <- matrix(0, nr, nr)
# For each row and for each column, assign values based on position
# These values are the product of two indexes
for(i in 1:dim(a)[1]) {
for(j in 1:dim(a)[2]) {
a[i,j] = if(i<j) {
0
}else if(i==j) {
1
}else {
3
}
}
}
# make into dataframe
mymat <- data.frame(mymat)
# create b variable
z <- rep(1:20)
# made it into lower diagonal matrix to make it easier to work with
b <- matrix(0, length(z), length(z))
b[lower.tri(b, diag = TRUE)] <- z[sequence(length(z):1)]
b
# create matrix for formula to operate with
# create variable "c"
c <- matrix(0, nr, nr)
# For each row and for each column, assign values based on position
# These values are the product of two indexes
for(i in 1:dim(c)[1]) {
for(j in 1:dim(c)[2]) {
mymat2[i,j] = if(i<j) {
0
}else if(i==j) {
1
}else {
5 # place holder for now
}
}
}
用于计算结果的公式，我认为由于 R 的矢量化而起作用
sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,])
然后我的问题是如何将其插入 if 语句中以创建递归定义的矩阵，如下所示
# calculate recursively defined lower triangular matrix
c <- matrix(0, nr, nr)
# For each row and for each column, assign values based on position
# These values are the product of two indexes
for(i in 1:dim(c)[1]) {
for(j in 1:dim(c)[2]) {
mymat2[i,j] = if(i<j) {
0
}else if(i==j) {
1
}else {
sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,]) # formula for calculation of values for lower triangular matrix
}
}
}
这给出了错误
Error in mymat2[i, j] <- if (i < j) { : number of items to replace is not a multiple of replacement length
我可以 link 一个 excel 电子表格，如果它有帮助的话，这个公式可以工作。只有通过大量手动输入等 excel 才能实现
使用真实数据的预期结果示例
a = 5x5 lower triangle matrix
0|0 |0 |0 |0
5|0 |0 |0 |0
5|0.56|0 |0 |0
5|0.20|0.61|0 |0
5|0.06|0.16|0.61|0
b = 1x5 matrix/vector
0.27917|0.499|0.83|1.191|1.48
c = recursive matrix results
1 |0 |0 |0 |0
6.36|1 |0 |0 |0
5.77|0.84|1 |0 |0
5.50|0.77|1.43|1 |0
5.3 |0.72|1.80|2.46|1
我如何在 excel
中计算“c”的示例
1 |0 |0 |0 |0
=(INDEX(a,2,1)*INDEX(b,1)*INDEX(c,1,1))/(INDEX(b,2)-INDEX(b,1))|1 |0 |0 |0
=(INDEX(a,3,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,3,2)*INDEX(b,2)*INDEX(c,2,1))/(INDEX(b,3)-INDEX(b,1))|=(INDEX(a,3,2)*INDEX(b,2)*INDEX(c,2,2))/(INDEX(b,3)-INDEX(b,2))|1 |0 |0
=(INDEX(a,4,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,4,2)*INDEX(b,2)*INDEX(c,2,1)+INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,1))/(INDEX(b,4)-INDEX(b,1))|=(INDEX(a,4,2)*INDEX(b,2)*INDEX(c,2,2)+INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,2))/(INDEX(b,4)-INDEX(b,2))|=(INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,3))/(INDEX(b,4)-INDEX(b,3))|1 |0
=(INDEX(a,5,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,5,2)*INDEX(b,2)*INDEX(c,2,1)+INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,1)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,1))/(INDEX(b,5)-INDEX(b,1)) |=(INDEX(a,5,2)*INDEX(b,2)*INDEX(c,2,2)+INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,2)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,2))/(INDEX(b,5)-INDEX(b,2))|=(INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,3)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,3))/(INDEX(b,5)-INDEX(b,3))|=(INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,4))/(INDEX(b,5)-INDEX(b,4))|1
以上代码显示了如何通过索引和手动设置大量单元格位置在 excel 中进行计算。
这是一个可能的动态规划（？）解决方案。我只验证了几个 'cells'，因此您需要验证它是否提供了您预期的结果。
# FUNCTIONS
vpad <- function(vec) {
extend <- 20 - length(vec)
c(vec, rep(0, extend))
}
clip <- function(x) { x[x >= 0 & x <= 20] }
pad_a <- function(a, iter) {
mat <- cbind(a[, iter:20], matrix(0, 20, iter-1))
mat * lower.tri(mat)
}
pad_b <- function(b, iter) {
tmp <- vpad(b[iter:20])
mat <- matrix(rep(tmp, each=20), 20, 20)
mat * lower.tri(mat)
}
pad_c <- function(c, iter) {
L <- 21 - iter
i <- clip(seq(L) + iter - 1)
j <- seq(i)
v <- vpad(mapply(function(i, j) c[i, j], i, j))
mat <- matrix(rep(v, each=20), 20, 20)
mat <- mat * (lower.tri(mat) + diag(20))
mat
end <- 21 - iter
mat1 <- rbind(matrix(0, iter - 1, 20), mat[1:end, ])
mat1 * lower.tri(mat1)
}
# MAKE FAKE DATA
set.seed(1)
m <- matrix(sample(400)/400, 20, 20)
a <- m * lower.tri(m)
b <- sample(20)/21
# INITIALIZE VARS
add_prev <- matrix(rep(0, 200), 20, 20)
init_mat <- lower.tri(m) + diag(20)
overwrite_mat <- init_mat
for (i in 1:20) {
overwrite_mat <- (pad_a(a, i) * pad_b(b, i) * pad_c(overwrite_mat, i)) + add_prev
add_prev <- overwrite_mat
}
denom1 <- matrix(rep(b, times=20), 20, 20)
denom1 <- denom1 * lower.tri(denom1) + diag(20)
denom2 <- matrix(rep(b, each=20), 20, 20)
denom2 <- denom2 * lower.tri(denom2) + diag(20)
final <- ((overwrite_mat / denom1) - denom2) + diag(20)
final <- replace(final, is.nan(final), 0)
输出（只显示左上角的（5,5）矩阵）
[,1] [,2] [,3] [,4] [,5]
[1,] 0.0000000 0.00000000 0.00000000 0.00000000
[2,] -0.5129073 0.00000000 0.00000000 0.00000000
[3,] 0.6738060 -0.05926190 0.00000000 0.00000000
[4,] 9.2856557 5.97660268 -0.15059524 0.00000000
[5,] 0.1971807 0.01591345 -0.14775033 -0.08718254
我可能把这个问题复杂化了。
问题是这样的
变量
a = lower triangular matrix/dataframe 20x20
b = a 1x20 matrix/vector
c = the previous row result of the formula (recursive bit)
我想要一个下三角矩阵 (I,j) 递归定义
伪代码
if(I<j) = 0
else if (I==j) = 1
else (I>j) = sum(a * b * c )/ b[I] - b[J]（其中 I 是当前行位置，J 是当前列位置）
我希望公式/R 的矢量化如何在较小的矩阵上工作以使生活更轻松的示例。
我=行
j = 列
b(我,j)
b(1,2) 指的是矩阵中的位置
示例矩阵的伪代码
Matrix C
1 | 2 |3 |4 |5 |6 |
1|i=j=1 |i<j=0 |...|...|...|...|
2|i>j=a(2,1)*b(1)*c(1,1)/b(2)-b(1) |i=j=1 |...|...|...|...|
3|i>j= (a(3,1)*b(1)*c(1,1) + a(3,2)*b(2)*c(2,1))/(b(3)-b(1) |i>j= a(3,2)*b(2)*c(2,2)/(b(3)-b(2) |...|...|...|...|
4|i>j = a(4,1)*b(1)*c(1,1)+a(4,2)*b(2)*c(2,1)+a(4,3)*b(3)*c(3,1))/b(4)-b(1) |i>j= a(4,2)*b(2)*c(2,2) + a(4,3)*b(3)*c(3,2)/b(4)-b(2) |...|...|...|...|
5|... |... |...|...|...|...|
6|... |... |...|...|...|...|
关于我目前的代码，
先创建变量如下图：
my_int <- 20
nr <- as.integer(my_int)
#create a n x n matrix with zeroes
a <- matrix(0, nr, nr)
# For each row and for each column, assign values based on position
# These values are the product of two indexes
for(i in 1:dim(a)[1]) {
for(j in 1:dim(a)[2]) {
a[i,j] = if(i<j) {
0
}else if(i==j) {
1
}else {
3
}
}
}
# make into dataframe
mymat <- data.frame(mymat)
# create b variable
z <- rep(1:20)
# made it into lower diagonal matrix to make it easier to work with
b <- matrix(0, length(z), length(z))
b[lower.tri(b, diag = TRUE)] <- z[sequence(length(z):1)]
b
# create matrix for formula to operate with
# create variable "c"
c <- matrix(0, nr, nr)
# For each row and for each column, assign values based on position
# These values are the product of two indexes
for(i in 1:dim(c)[1]) {
for(j in 1:dim(c)[2]) {
mymat2[i,j] = if(i<j) {
0
}else if(i==j) {
1
}else {
5 # place holder for now
}
}
}
用于计算结果的公式，我认为由于 R 的矢量化而起作用
sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,])
然后我的问题是如何将其插入 if 语句中以创建递归定义的矩阵，如下所示
# calculate recursively defined lower triangular matrix
c <- matrix(0, nr, nr)
# For each row and for each column, assign values based on position
# These values are the product of two indexes
for(i in 1:dim(c)[1]) {
for(j in 1:dim(c)[2]) {
mymat2[i,j] = if(i<j) {
0
}else if(i==j) {
1
}else {
sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,]) # formula for calculation of values for lower triangular matrix
}
}
}
这给出了错误
Error in mymat2[i, j] <- if (i < j) { : number of items to replace is not a multiple of replacement length
我可以 link 一个 excel 电子表格，如果它有帮助的话，这个公式可以工作。只有通过大量手动输入等 excel 才能实现
使用真实数据的预期结果示例
a = 5x5 lower triangle matrix
0|0 |0 |0 |0
5|0 |0 |0 |0
5|0.56|0 |0 |0
5|0.20|0.61|0 |0
5|0.06|0.16|0.61|0
b = 1x5 matrix/vector
0.27917|0.499|0.83|1.191|1.48
c = recursive matrix results
1 |0 |0 |0 |0
6.36|1 |0 |0 |0
5.77|0.84|1 |0 |0
5.50|0.77|1.43|1 |0
5.3 |0.72|1.80|2.46|1
我如何在 excel
中计算“c”的示例1 |0 |0 |0 |0
=(INDEX(a,2,1)*INDEX(b,1)*INDEX(c,1,1))/(INDEX(b,2)-INDEX(b,1))|1 |0 |0 |0
=(INDEX(a,3,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,3,2)*INDEX(b,2)*INDEX(c,2,1))/(INDEX(b,3)-INDEX(b,1))|=(INDEX(a,3,2)*INDEX(b,2)*INDEX(c,2,2))/(INDEX(b,3)-INDEX(b,2))|1 |0 |0
=(INDEX(a,4,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,4,2)*INDEX(b,2)*INDEX(c,2,1)+INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,1))/(INDEX(b,4)-INDEX(b,1))|=(INDEX(a,4,2)*INDEX(b,2)*INDEX(c,2,2)+INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,2))/(INDEX(b,4)-INDEX(b,2))|=(INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,3))/(INDEX(b,4)-INDEX(b,3))|1 |0
=(INDEX(a,5,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,5,2)*INDEX(b,2)*INDEX(c,2,1)+INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,1)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,1))/(INDEX(b,5)-INDEX(b,1)) |=(INDEX(a,5,2)*INDEX(b,2)*INDEX(c,2,2)+INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,2)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,2))/(INDEX(b,5)-INDEX(b,2))|=(INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,3)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,3))/(INDEX(b,5)-INDEX(b,3))|=(INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,4))/(INDEX(b,5)-INDEX(b,4))|1
以上代码显示了如何通过索引和手动设置大量单元格位置在 excel 中进行计算。
这是一个可能的动态规划（？）解决方案。我只验证了几个 'cells'，因此您需要验证它是否提供了您预期的结果。
# FUNCTIONS
vpad <- function(vec) {
extend <- 20 - length(vec)
c(vec, rep(0, extend))
}
clip <- function(x) { x[x >= 0 & x <= 20] }
pad_a <- function(a, iter) {
mat <- cbind(a[, iter:20], matrix(0, 20, iter-1))
mat * lower.tri(mat)
}
pad_b <- function(b, iter) {
tmp <- vpad(b[iter:20])
mat <- matrix(rep(tmp, each=20), 20, 20)
mat * lower.tri(mat)
}
pad_c <- function(c, iter) {
L <- 21 - iter
i <- clip(seq(L) + iter - 1)
j <- seq(i)
v <- vpad(mapply(function(i, j) c[i, j], i, j))
mat <- matrix(rep(v, each=20), 20, 20)
mat <- mat * (lower.tri(mat) + diag(20))
mat
end <- 21 - iter
mat1 <- rbind(matrix(0, iter - 1, 20), mat[1:end, ])
mat1 * lower.tri(mat1)
}
# MAKE FAKE DATA
set.seed(1)
m <- matrix(sample(400)/400, 20, 20)
a <- m * lower.tri(m)
b <- sample(20)/21
# INITIALIZE VARS
add_prev <- matrix(rep(0, 200), 20, 20)
init_mat <- lower.tri(m) + diag(20)
overwrite_mat <- init_mat
for (i in 1:20) {
overwrite_mat <- (pad_a(a, i) * pad_b(b, i) * pad_c(overwrite_mat, i)) + add_prev
add_prev <- overwrite_mat
}
denom1 <- matrix(rep(b, times=20), 20, 20)
denom1 <- denom1 * lower.tri(denom1) + diag(20)
denom2 <- matrix(rep(b, each=20), 20, 20)
denom2 <- denom2 * lower.tri(denom2) + diag(20)
final <- ((overwrite_mat / denom1) - denom2) + diag(20)
final <- replace(final, is.nan(final), 0)
输出（只显示左上角的（5,5）矩阵）
[,1] [,2] [,3] [,4] [,5]
[1,] 0.0000000 0.00000000 0.00000000 0.00000000
[2,] -0.5129073 0.00000000 0.00000000 0.00000000
[3,] 0.6738060 -0.05926190 0.00000000 0.00000000
[4,] 9.2856557 5.97660268 -0.15059524 0.00000000
[5,] 0.1971807 0.01591345 -0.14775033 -0.08718254