在R中求解方程数多于未知参数的非线性方程组

问题本质 #

你碰到的Length of fn result <> length of x!错误，核心原因是数值求解器要求方程组的方程数量必须和待求未知参数数量严格相等——你定义了7个方程，但只有4个未知参数p[1]-p[4]，维度不匹配导致求解失败。

可行解决思路 #

1. 剔除冗余方程，匹配维度 #

观察你的前6个方程，它们是不同组之间的差值等式，存在明显的线性冗余：

• 比如方程4（组2-组3的差）可以由方程1（组1-组2的差）减去方程2（组1-组3的差）推导得出，属于非独立方程；同理方程5、6也是冗余项。

只需保留3个独立的组间差值方程，加上第7个概率和为1的约束，就能凑齐4个方程，刚好匹配4个未知参数。修改后的myFun代码示例：

myFun<- function(p) {
  # 预计算每个组的核心项，减少重复运算
  term1 <- sum(gamma[,1]*x2)/p[1] - sum(gamma[,1]*x1)/(1-p[1])
  term2 <- sum(gamma[,2]*x2)/p[2] - sum(gamma[,2]*x1)/(1-p[2])
  term3 <- sum(gamma[,3]*x2)/p[3] - sum(gamma[,3]*x1)/(1-p[3])
  term4 <- sum(gamma[,4]*x2)/p[4] - sum(gamma[,4]*x1)/(1-p[4])
  
  y <- numeric(4)
  y[1] = term1 - term2  # 保留组1与组2的差值方程
  y[2] = term1 - term3  # 保留组1与组3的差值方程
  y[3] = term1 - term4  # 保留组1与组4的差值方程
  y[4] = 1 - (p[1] + p[2] + p[3] + p[4])  # 概率和为1的约束
  return(y)
}

p=c(0.2,0.4,0.1,0.3)
pi=nleqslv(p,myFun)$x

2. 超定方程组用最小二乘拟合 #

如果必须保留所有7个方程（说明这是超定方程组，方程数多于未知数），则不能用nleqslv这类求精确解的工具，改用最小二乘方法找到使所有方程残差平方和最小的参数。可以用R基础包的optim函数实现：

# 定义残差平方和计算函数
residual_sum <- function(p) {
  term1 <- sum(gamma[,1]*x2)/p[1] - sum(gamma[,1]*x1)/(1-p[1])
  term2 <- sum(gamma[,2]*x2)/p[2] - sum(gamma[,2]*x1)/(1-p[2])
  term3 <- sum(gamma[,3]*x2)/p[3] - sum(gamma[,3]*x1)/(1-p[3])
  term4 <- sum(gamma[,4]*x2)/p[4] - sum(gamma[,4]*x1)/(1-p[4])
  
  # 计算所有方程的残差
  y1 = term1 - term2
  y2 = term1 - term3
  y3 = term1 - term4
  y4 = term2 - term3
  y5 = term2 - term4
  y6 = term3 - term4
  y7 = 1 - (p[1]+p[2]+p[3]+p[4])
  
  # 返回残差平方和
  sum(c(y1,y2,y3,y4,y5,y6,y7)^2)
}

# 初始值设置，同时约束参数在(0,1)区间（避免分母为0）
p_init <- c(0.2,0.4,0.1,0.3)
result <- optim(p_init, residual_sum, method = "L-BFGS-B", 
                lower = rep(1e-6,4), upper = rep(1-1e-6,4))
pi <- result$par  # 得到最小二乘最优解

这里用L-BFGS-B方法是为了限制p的取值范围在(0,1)之间，符合概率的定义，同时避免计算中出现分母为0的错误。

3. 模型逻辑校验 #

从你的方程组来看，核心需求是让所有组的sum(gamma[,k]*x2)/p[k] - sum(gamma[,k]*x1)/(1-p[k])值相等，再加上概率和为1的约束。这种场景下，保留3个独立差值方程+1个和约束是最合理的选择，既满足数值求解的维度要求，也避免了冗余计算带来的效率损耗。

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