[R] solve nonlinear equation using BBsolve

[Berend Hasselman]

Roslina Zakaria wrote:
> 
> Hi r-users,
> 
> I would like to solve system of nonlinear equation using BBsolve function
> and 
> below is my code.  I have 4 parameters and I have 4 eqns.
> 
> mgf_gammasum <- function(p)
> {
> t  <- rep(NA, length(p))
> mn <- 142.36
> vr <- 9335.69
> sk <- 0.8139635
> kur <- 3.252591
> rh  <- 0.896
> # cumulants
> k1 <- p[1]*(p[2]+p[3])
> k2 <- p[1]*(2*p[2]*p[3]*p[4] +p[2]^2+p[3]^2)
> k3 <- 2*p[1]*(p[2]+p[3])*(p[2]^2 + p[3]^2 - p[2]*p[3] + 3*p[2]*p[3]*p[4])
> k4 <- 6*p[1]*((p[2]+p[3])^2*(p[2]^2 + p[3]^2 - 2*p[2]*p[3] +
> 4*p[2]*p[3]*p[4])+ 
> 2*p[2]^2*p[3]^2*(1-p[4])^2)
> t[1] <- k1 - mn
> t[2] <- k2 - vr
> t[3] <- k3/(k2^1.5) - sk
> t[4] <- k4/(k2^2)   - kur
> t
> }
> 
> I tried this 
> p0 <- rep(0, 4)
> BBsolve(par = p0, fn = mgf_gammasum)
> dfsane(par = p0, fn = mgf_gammasum, control = list(trace = FALSE))
> sane(par = p0, fn = mgf_gammasum, control = list(trace = FALSE))
> 
> and got the error message:
> 
>> BBsolve(par = p0, fn = mgf_gammasum)
> Error in optim(par = par, fn = U, method = "Nelder-Mead", control =
> list(maxit = 
> 100),  : 
> 

I did this:
 
> p0 <- rep(0, 4) 
> mgf_gammasum(p0)
[1]  -142.36 -9335.69      NaN      NaN

Your starting values are invalid.

We've had a discussion about this function before (somewhere in May of this
year).

I tried this

p0 <- c(1,1,1,1)
mgf_gammasum(p0)

BBsolve(par = p0, fn = mgf_gammasum) 


==> BBsolve can't find a solution
So let's try something else

library(nleqslv)
nleqslv(p0,mgf_gammasum)

==> Oops, singular jacobian

# Try to calculate the jacobian in the starting point

n <- length(p0)
J <- matrix(data=rep(0,n*n),nrow=n)

eps <- 1e-8
f0 <- mgf_gammasum(p0)
J[,1] <- (mgf_gammasum(p0+c(eps+eps*p0[1],0  ,0  ,0  )) -
f0)/(eps+eps*p0[1])
J[,2] <- (mgf_gammasum(p0+c(0  ,eps+eps*p0[2],0  ,0  )) -
f0)/(eps+eps*p0[2])
J[,3] <- (mgf_gammasum(p0+c(0  ,0  ,eps+eps*p0[3],0  )) -
f0)/(eps+eps*p0[3])
J[,4] <- (mgf_gammasum(p0+c(0  ,0  ,0  ,eps+eps*p0[4])) -
f0)/(eps+eps*p0[4])
                                               
J

svd(J)


Your system of equations is almost singular (at least in this starting
point).
Solution is to find a better starting point assuming there is one.

Berend
-- 
View this message in context: http://r.789695.n4.nabble.com/solve-nonlinear-equation-using-BBsolve-tp3052167p3052178.html
Sent from the R help mailing list archive at Nabble.com.