[R] help about solving two equations

[Berend Hasselman]

Shaoqiong Zhao wrote:
> 
> ....
> I am not sure if this is correct and also this can only solve one row. How
> to get the whole 1000 rows of p and q?
> 

You can try something 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
}

xstart <- c(1,1)

Krep <- 100
s1 <- rnorm(Krep,mean=-2,sd=.5)
s2 <- rnorm(Krep,mean=-4,sd=.5)
sz <- cbind(s1,s2)
xans <- matrix(ncol=2,nrow=Krep)

system.time( {
    noncvg <- 0
    for(k in seq(Krep)) {
       ans <- nleqslv(xstart,fn, s=sz[k,],method="Newton")
       xans[k,] <- ans$x
       if(ans$termcd>1){noncvg<-noncvg+1}
    }
    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=sz[k,],control=list(trace=FALSE))
       xans[k,] <- ans$par
       if(ans$convergence>0){noncvg<-noncvg+1}
    }
    msg <- sprintf("Non convergence %d ==> %.2f%%\n", noncvg,
noncvg/Krep*100)
    print(msg)
})

It seems that nleqslv finds more solutions than dfsane and is also quite a
bit quicker.

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