使命
我正在尝试找到一种在 R 中进行迭代比例拟合的方法。该过程的逻辑是这样的：有一个表格，其中包含例如一些变量的样本分布。让我们说它是这个：
sample1 <- structure(c(6L, 14L, 46L, 16L, 6L, 21L, 62L, 169L, 327L, 174L,
44L, 72L, 43L, 100L, 186L, 72L, 23L, 42L), .Dim = c(6L, 3L), .Dimnames = list(
c("Primary", "Lowersec", "Highersec", "Highershort", "Higherlong",
"University"), c("B", "F", "W")))
另一个表来自其他来源，比如另一个示例：
sample2 <- structure(c(171796L, 168191L, 240671L, 69168L, 60079L, 168169L,
954045L, 1040981L, 1872732L, 726410L, 207366L, 425786L, 596239L,
604826L, 991640L, 323215L, 134066L, 221696L), .Dim = c(6L, 3L
), .Dimnames = list(c("Primary", "Lowerse", "Highersec", "Highershort",
"Higherlong", "University"), c("B", "F", "W")))
现在，我们想保留在 中找到的变量之间的关系sample1，但我们想将这些关系应用于我们在 中找到的边际分布sample2。迭代比例拟合按照此处所述进行此操作（我无法提供更好的解释）。我曾尝试在 LEM 中执行此操作，结果如下：
B F W
Primary 124204.64 960173.6 637701.7
Lowerse 119749.12 1081459.0 612789.9
Highersec 336934.21 1792001.6 976107.2
Highershort 90512.27 736464.1 291816.6
Higherlong 43486.91 238593.0 119431.0
University 163186.85 418628.6 233835.5
我不是 100% 肯定这个结果，但很有可能（比如 99%）它是。第一个表中的优势比保留在结果表中，而边际分布（行和列总和）与第二个输入表相同。
问题
奇怪的是，这种非常有用的算法在 R 中并不容易获得，至少不是以用户友好的形式。一个可能相关的功能是cat::ipf()。但是，我无法弄清楚如何使用该margins=参数。在这个问题上，我当然并不孤单。帮助示例使用了一个 3 维表，这使事情变得更加混乱。
此外，还有一些用户编写的功能，一个在这里，一个在这里找到。不幸的是，第一个给出了错误的结果。第二个也是非常不透明的，需要专门预先格式化的 CSV 文件作为输入，而不是 R 矩阵对象。
问题
- 谁能解释一下如何实际使用该
cat::ipf()功能？ - 是否有任何替代函数来实现 IPF 调整任务，使用矩阵作为输入？
- （已解决）可以修复此功能以提供正确的结果吗？
谢谢你。
附录：我能够从（3）中的函数获得正确的输出。经过一番考虑，结果发现该函数不接受矩阵作为边际分布的输入，而只是这些边际分布的列表。所以实际上这个问题已经解决了。然而，1 和 2 的正确答案将对更大的社区有用，因为 IPF 在对数线性模型中非常重要。
对于那些将来搜索的人来说，在这里复制一个有效的 IPF 函数似乎是一个好主意：
ipf <- function(Margins_, seedAry, maxiter=100, closure=0.001) {
#Check to see if the sum of each margin is equal
MarginSums. <- unlist(lapply(Margins_, sum))
if(any(MarginSums. != MarginSums.[1])) warning("sum of each margin
not equal")
#Replace margin values of zero with 0.001
Margins_ <- lapply(Margins_, function(x) {
if(any(x == 0)) warning("zeros in marginsMtx replaced with
0.001")
x[x == 0] <- 0.001
x
})
#Check to see if number of dimensions in seed array equals the number of
#margins specified in the marginsMtx
numMargins <- length(dim(seedAry))
if(length(Margins_) != numMargins) {
stop("number of margins in marginsMtx not equal to number of
margins in seedAry")
}
#Set initial values
resultAry <- seedAry
iter <- 0
marginChecks <- rep(1, numMargins)
margins <- seq(1, numMargins)
#Iteratively proportion margins until closure or iteration criteria are met
while((any(marginChecks > closure)) & (iter < maxiter)) {
for(margin in margins) {
marginTotal <- apply(resultAry, margin, sum)
marginCoeff <- Margins_[[margin]]/marginTotal
marginCoeff[is.infinite(marginCoeff)] <- 0
resultAry <- sweep(resultAry, margin, marginCoeff, "*")
marginChecks[margin] <- sum(abs(1 - marginCoeff))
}
iter <- iter + 1
}
#If IPF stopped due to number of iterations then output info
if(iter == maxiter) cat("IPF stopped due to number of iterations\n")
#Return balanced array
resultAry
}
请注意，此函数不接受边际分布的矩阵。以下将有所帮助：
m1 <- rowSums(sample2)
m2 <- colSums(sample2)
m <- list(m1,m2)
然后提供m作为第一个参数。