Finding difference between Mathematica code and R code for solving set of equations

I have this Mathematica code trying to solve a set of equations:

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]

And gives me the following output:

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

I am trying to do the same thing using R Software instead, by using the following code:

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)

These should be equivalent, however the R output gives me the following output:

$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

Are the in fact the same? How come Mathematica manages to find a solution while my R code does not?

Note that the values of a and p (since I presume that you are assuming a mixture distribution) given by Mathematica do not make sense theoretically.

Hence, R does not find this solutions since you are using dbinom and it does not accept negative probabilities. Thus, there seems to be no solution to your problem.