[R] help about solving two equations

[Berend Hasselman]

Ravi Varadhan wrote:
> 
> require(BB)
> 
> fn <- function(x, s){
> f <- rep(NA, length(x))
> f[1] <- digamma(x[1]) - digamma(x[1]+x[2]) - s[1]  
> f[2] <- digamma(x[2]) - digamma(x[1]+x[2]) - s[2] 
> f  
> } 
> 
>  nr <- 10
> smat <- matrix(runif(2*nr, -3, -1), nr, 2)  
> soln <- matrix(NA, nr, 2)
> 
> for (i in 1:nr) {
> ans <- dfsane(par=c(1,1), fn=fn, s=smat[i, ], control=list(trace=FALSE))  
> soln[i, ] <- ans$par
> }
> 
> soln # your solution
> 

As I mentioned in a previous post to the original question, it is necessary
to record if an algorithm actually found a solution. I would also advise you
to try an alternative to BB; it's a bit quicker in your case and seems to
fail less frequently.
Like this:

library(nleqslv)
library(BB)

fn <- function(x, s){
    f <- rep(NA, length(x))
    f[1] <- digamma(x[1]) - digamma(x[1]+x[2]) - s[1]
    f[2] <- digamma(x[2]) - digamma(x[1]+x[2]) - s[2]
    f
}

s <- c(-2,-4)

xstart <- c(1,1)

Krep <- 100
smat <- matrix(runif(2*Krep, -3, -1), Krep, 2)   
sobb <- matrix(NA, Krep, 2) 
sonl <- matrix(NA, Krep, 2) 

system.time( {
    noncvg <- 0
    for(k in seq(Krep)) {
       ans <- nleqslv(xstart,fn, s=smat[k,],method="Newton")
       if(ans$termcd>1){noncvg<-noncvg+1;sonl[k,] <- c(NA,NA)} else sonl[k,]
<- ans$x
    }
    msg <- sprintf("Non convergence %d ==> %.2f%%\n", noncvg,
noncvg/Krep*100)
    print(msg)
})

system.time( {
    noncvg <- 0
    for(k in seq(Krep)) {
       ans <- dfsane(par=xstart, fn=fn,
s=smat[k,],control=list(trace=FALSE))
       if(ans$convergence>0){noncvg<-noncvg+1;sobb[k,] <- c(NA,NA)} else
sobb[k,] <- ans$par
    }
    msg <- sprintf("Non convergence %d ==> %.2f%%\n", noncvg,
noncvg/Krep*100)
    print(msg)
})

sonl
sobb     

Berend

-- 
View this message in context: http://n4.nabble.com/help-about-solving-two-equations-tp1588313p1591524.html
Sent from the R help mailing list archive at Nabble.com.