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. ______________________________________________ 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.