如何在循环和 if 语句中使用递归公式答案

【问题标题】：How to use recursive formula in loop and if statement如何在循环和 if 语句中使用递归公式
【发布时间】：2020-07-22 01:17:59
【问题描述】：

我可能把这个问题复杂化了。

问题是这样的

变量

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 的矢量化如何在较小的矩阵上显示以使生活更轻松的示例。

I = 行
j = 列
b(I,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

如果有帮助，我可以链接一个 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中通过索引和手动设置大量的单元格位置来计算它。

【问题讨论】：

在您的矩阵中，您能否修复单元格（4,2）并显示单元格（4,1）应该是什么样子？
如果您包含一个具有预期输出的 5x5 矩阵，将会很有帮助。我不确定您的sum(...) 线路是否正在做您希望它做的事情。此外，sum(...) 行中的分母将使结果成为您示例中长度为 20 的向量。错误是 mymat2[i, j] 应该只分配一个值，但这会分配 20。
您好，感谢您的评论，我已经进行了更改并添加了数据示例。 @CPak 我已更正 (4,2) 感谢您发现并添加到 (4,1)。
@Cole 我添加了一个数据应该做什么的例子。我认为您对总和的理解是正确的，但解释错误。
对不起，我看不出它是如何概括的。不确定这是否是您正在做的线性代数，但使用与此相关的术语。最后，一些用于演示制作矩阵的初始代码可以简化为n = 5L; mat = diag(n); mat[lower.tri(mat)] = 3。祝你帖子的其余部分好运。

标签： r recursion matrix formula recursive-datastructures

【解决方案1】：

这是一个可能的动态规划（？）解决方案。我只验证了几个“单元格”，因此您需要验证它是否提供了您预期的结果。

# 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

【讨论】：

非常感谢！真正令人印象深刻的解决方案。