Stack Overflow Asked by Biggerthanpenny on November 17, 2021
I may be overcomplicating this problem.
The problem is as such
Variables
a = lower triangular matrix/dataframe 20x20
b = a 1x20 matrix/vector
c = the previous row result of the formula (recursive bit)
I want to have a lower triangular matrix (I,j) where it is recursively defined that
pseudocode
if(I<j) = 0
else if (I==j) = 1
else (I>j) = sum(a * b * c )/ b[I] - b[J]
(where I is current row position and J is current column position)
Example of how I hope the formula/ R’s vectorisation will work shown on a smaller matrix to make life easier.
I = row
j = column
b(I,j)
b(1,2) refers to position in matrix
Pseudo code of example matrix
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|... |... |...|...|...|...|
Onto the code I have so far,
Create variable first as shown below:
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
}
}
}
Formula for calculating results, I think works thanks to R’s vectorisation
sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,])
My issue is then how to insert this in the if statements for the creation of the recursively defined matrix as shown below
# 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
}
}
}
This gives error
Error in mymat2[i, j] <- if (i < j) { : number of items to replace is not a multiple of replacement length
I can link an excel spreadsheet where this formula works if it would help. It is only achievable in excel by a lot of manual inputs etc.
Example of expected results using real data
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
Example of how I calculated "c" in excel
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
Above code shows how it can be calculated in excel via index and manually setting a large amount of cell positions.
Here's a possibly dynamic programming(?) solution. I've only verified a few 'cells' so you'll need to verify that it's giving results you expected.
# 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)
Output (just showing upper left (5,5) matrix)
[,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
Answered by CPak on November 17, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP