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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.