[R] BB package

[Berend Hasselman]

Roslina Zakaria wrote:
> 
> Hi I would like to solve a system of nonlinear equations below using
> dfsane function
>  
> mn <- 142.36; vr <- 9335.69 ; sk <- 0.81;  kur <- 0.25
> test_fn <- function(p) 
> {    f <- rep(NA, length(p))
>   
>   f[1] <- p[1]*(p [2]+p[3])- mn
>   f[2] <- - vr + 2*p[1]*p[2]*p[3]*(p[4]-1)+p[1]*(p[2]+p[3])^2 
>   f[3] <- - sk + (p[1]*(p
> [2]+p[3])^3*(p[1]+1)*(p[1]+2)-6*p[1]*p[2]*p[3]*(p[2]+ p[3])*(1- p[4])*( p
> [1]+1))/(p[1]*(p[2]+p[3])^2*(p[2]+p[3])-2*p[1]*p[2]*p[3]*(1-p[4]))^1.5 
>   f[4] <- - kur + (p[1]*(p[2]+p[3])^4*(p[1]+3)*(p[1]+2)*(p[1]+1)
> -12*p[1]*p[2]*p[3]*(1-p[4])*(p[2]+p[3])^2*(p[1]+2)*(p[1]+1)+12*p[1]*p[2]^2*p[3]^2*(1-p[4])^2*(p[1]+1))/(p[1]*(p[2]+p[3])^2*(p[2]+p[3])-2*p[1]*p[2]*p[3]*(1-p[4]))^2 
>   f
>   }
>  
> p0 <- c(1.3, 50,60,0.8)
> dfsane(par=p0, fn = test_fn, control=list(trace=FALSE))
> 

I tried solving your system of equations with package nleqslv.
It reports an ill conditioned/singular jacobian and cannot solve the system.
I did a rough calculation of the jacobian in the initial point you provided.

> n <- length(p0)
> J <- matrix(data=rep(0,n*n),nrow=n)
> 
> eps <- 1e-8
> f0 <- test_fn(p0)
> J[,1] <- (test_fn(p0+c(eps,0  ,0  ,0  )) - f0)/eps
> J[,2] <- (test_fn(p0+c(0  ,eps,0  ,0  )) - f0)/eps
> J[,3] <- (test_fn(p0+c(0  ,0  ,eps,0  )) - f0)/eps
> J[,4] <- (test_fn(p0+c(0  ,0  ,0  ,eps)) - f0)/eps
> 
> J
             [,1]          [,2]          [,3]         [,4]
[1,] 1.100000e+02  1.300000e+00  1.300000e+00 0.000000e+00
[2,] 1.090000e+04  2.548002e+02  2.600003e+02 7.800000e+03
[3,] 2.014600e-03 -7.268630e-05 -7.097656e-05 2.569023e-03
[4,] 4.256762e-04 -3.317069e-05 -3.225198e-05 1.378761e-03
> 

Columns 2 and 3 are almost identical.
p[2] and p[3] seem to be interchangeable.

If you do svd(J) then you will see the singular values 

$d
[1] 1.340860e+04 6.401512e+01 2.129760e-04 5.305005e-10

smallest singular value/largest singular value is very close to machine
precision.

Your system of equations is almost singular.

Berend
-- 
View this message in context: http://r.789695.n4.nabble.com/BB-package-tp2074937p2076045.html
Sent from the R help mailing list archive at Nabble.com.