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.
