查找用于求解方程组的Mathematica代码和R代码之间的差异

Question

我有这个Mathematica代码试图求解一组方程：

f[k_, n_, p_, a_, b_] :=p*(Binomial[n, k]*a^k*(1 - a)^(n - k)) + (1 - p)*Binomial[n, k]*b^k*(1 - b)^(n - k));
    mom1 = Sum[k^1 *f[k, n, p, a, b], {k, 0, n}] - 3.3;
    mom2 = Sum[k^2*f[k, n, p, a, b], {k, 0, n}] - 13.04;
    mom3 = Sum[k^3*f[k, n, p, a, b], {k, 0, n}] - 58.08;
    mom4 = Sum[k^4*f[k, n, p, a, b], {k, 0, n}] - 281.96;
estimate = NSolve[{mom1 == 0, mom2 == 0, mom3 == 0, mom4 == 0, n > 6}, {p, a, b, n}, Reals]

并给我以下输出：

{{p -> -0.0000925709, a -> -1.15159, b -> 0.343157, n -> 9.61271}}

我试图通过使用以下代码使用R软件来做同样的事情：

init=c(0.2, 9, 0.2, 0.2);
GetMoment<-function(x, k, ord){
  (k^ord)*DensityFn(x, k)
}
DensityFn<-function(x, k){
  x[1]*dbinom(k, x[2], x[3])+(1-x[1])*dbinom(k, x[2], x[4])
}
target<-function(x){
  y <- integer(4)
  x[2]=floor(x[2]);
  y[1]=sum(GetMoment(x, 0:x[2], 1))-3.3;
  y[2]=sum(GetMoment(x, 0:x[2], 2))-13.04;
  y[3]=sum(GetMoment(x, 0:x[2], 3))-58.08;
  y[4]=sum(GetMoment(x, 0:x[2], 4))-281.96;
  y
}
out=nleqslv(init, target);
print(out)

这些应该等效，但是R输出为我提供以下输出：

$x
[1] 0.2 9.0 0.2 0.2
$fvec
[1]    -2.438   -24.314  -226.394 -2120.652
$termcd
[1] 6
$message
[1] "Jacobian is singular (1/condition=0.0e+000) (see allowSingular option)"
$scalex
[1] 1 1 1 1
$nfcnt
[1] 0
$njcnt
[1] 1
$iter
[1] 1

实际上是一样的吗？Mathematica如何设法找到解决方案，而我的R代码没有？

Answer 1

请注意， a和 p的值(因为我认为您假设您假设的混合物分布(，而Mathematica给出了理论上没有意义的。

因此，r找不到此解决方案，因为您使用的是dbinom，并且不接受负概率。因此，似乎没有解决您的问题的方法。