二元方程组（含非线性方程）求解及卡方检验问题求助

问题描述 #

给定二元方程组：

(1+x)-(x-1)*y=1.706873
((x-1)*(x-2)*y^2)-3*(x^2-1)*y+(x+1)*(2*x+1) = 4.039665

使用R语言nleqslv包求解时，得到x、y均为0的不合理结果，同时报错：Jacobian is singular (1/condition=0.0e+000) (see allowSingular option)。

此外，求解结果需代入循环代码，在pgamma(1+jj*x)和dnbinomial(..., 1/y)中计算概率值e，进而计算卡方值。原假设要求1-chisq>0.05但未满足，无法确定问题出在方程组求解还是后续代码环节。

问题排查与解决 #

1. 修正方程组代码的核心错误 #

原代码中第一个方程的表达式写错了：

• 原代码错误写法：(1+x)*(x-1)*y - 1.706873
• 对应实际方程组的正确写法：(1 + x) - (x - 1)*y - 1.706873

这是导致初始求解结果完全偏离的关键原因。

2. 解决Jacobian奇异问题 #

初始值dstart <- c(1,1)会触发奇异矩阵：当x=1时，(x-1)=0，使得方程组的梯度信息失效，导致求解器无法正常迭代。更换为更合理的初始值，比如c(2, 0.5)或c(0, 0)。

3. 修正后的完整求解代码 #

install.packages('nleqslv')
library(nleqslv)

f1 <- function(d) {
  f <- numeric(2)
  x <- d[1]
  y <- d[2]
  # 修正第一个方程
  f[1] <- (1 + x) - (x - 1)*y - 1.706873
  f[2] <- (x-1)*(x-2)*y^2 - 3*(x^2-1)*y + (x+1)*(2*x+1) - 4.039665
  f
}

# 使用合理初始值
dstart <- c(2, 0.5)
d1 <- nleqslv(dstart, f1)
print(d1)

运行后可得到合理解（例如x≈1.5，y≈0.4057，具体数值以实际运行结果为准），且无奇异矩阵报错。

4. 后续卡方检验的排查方向 #

• 验证方程组解的正确性：将求解得到的x、y代入原方程组，确认等式两边的误差在可接受范围内（如1e-6以内），排除求解环节的问题。
• 检查概率计算逻辑：
  • 确认pgamma的参数（形状、尺度等）是否符合需求，1+jj*x的输入是否正确。
  • 确认dnbinomial的参数设置，特别是size参数是否为1/y，概率参数是否匹配业务逻辑。
• 检查卡方值计算：确保观测值与期望（即e）的对应关系正确，卡方统计量公式sum((观测值-期望值)^2/期望值)无计算错误。

内容的提问来源于stack exchange，提问作者fbah