求解联立方程 (>2) 而无需转换为 R 或 Python 中的矩阵

Question

我有一组 4 个联立方程：

0.059z = x
0.06w = y
z+w = 8093
x+y = 422

到目前为止，我发现的所有解决方案似乎都适用于每个方程中存在所有变量的方程，然后转换为矩阵并使用求解函数。

有没有更简单的方法在 R 或 Python 中使用原始形式的方程来做到这一点？

另外，如何确保解决方案中只返回正数？

希望这是有道理的...非常感谢您的帮助

Answer 1

您可以为此使用sympy：

from sympy import symbols, linsolve, Eq
x,y,z,w = symbols('x y z w')
linsolve([Eq(0.059*z, x), Eq(0.06*w, y), Eq(z+w, 8093), Eq(x+y, 422)], (x, y, z, w))

输出：

关于负值的 cmets - 方程组只有一个解，它对 y 和 w 具有负值。如果有多个解决方案，sympy 将返回它们，您可以将解决方案从那里过滤为仅正值。

Answer 2

在R，也许你可以像下面这样尝试：

library(rootSolve)
library(zeallot)

model <- function(v){
  c(x,y,z,w) %<-% v
  return(c(0.059*z-x, 0.06*w-y, z+w-8093, x+y-422))
}
res <- multiroot(f = model, start = c(0,0,0,0))

那么你可以得到解决方案res为

> res
[1]   3751.22  -3329.22  63580.00 -55487.00

Answer 3

这里发生了一些事情。首先正如 CDJB 所说：如果有任何积极的解决方案，那么 sympy 会找到它们。我搜索了这些数字并找到了this paper，这表明您应该使用 7088 而不是 8093。我们可以进行快速的健全性检查：

def pct(value):
    return f"{value:.1%}"

print(pct(422 / 8093))  # ~5.2%
print(pct(422 / 7088))  # ~6.0%

确认您将难以将平均 ~5.9% 和 ~6.0% 提高到 ~5.2%，并在其他答案中解释否定解决方案。此外，这些可能是计数，因此您的所有变量也需要是整数。

一旦使用了这个正确的分母，我会评论说有 许多 解决方案（我的计数是 11645）例如：

cases = [1, 421]
pop = [17, 7071]

rates = [pct(c / p) for c, p in zip(cases, pop)]

提供适当的输出，如下所示：

cases = [2, 420]
pop = [34, 7054]

这是因为数据被四舍五入到小数点后两位。您可能也不想使用上述任何一种，它们只是我得到的前两个有效解决方案。

我们可以定义一个 Python 函数来枚举所有解决方案：

from math import floor, ceil

def solutions(pop, cases, rate1, rate2, err):
    target = (pct(rate1), pct(rate2))
    for pop1 in range(1, pop):
        pop2 = pop - pop1
        c1_lo = ceil(pop1 * (rate1 - err))
        c1_hi = floor(pop1 * (rate1 + err))

        for c1 in range(c1_lo, c1_hi+1):
            c2 = cases - c1
            if (pct(c1 / pop1), pct(c2 / pop2)) == target:
                yield c1, c2, pop1, pop2

all_sols = list(solutions(7088, 422, 0.059, 0.060, 0.0005))

我从上面得到了 11645 的计数。

不确定对此有何建议，但您可以发送bootstrap 来查看您的统计数据因不同的解决方案而变化多少。另一种选择是进行贝叶斯分析，它可以让您将先验置于人口规模之上，从而大大减少这一点。