Here is how you can solve: 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 }
require(BB) # load this package for the nonlinear solver s <- c(-2, -4) # one row of s1 and s2 ans <- dfsane(par=c(1,1), fn=fn, s=s) ans$par # solutions for p and q You can then loop through for each row of s1 and s2 and solve it to get corresponding p and q. Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: Shaoqiong Zhao <zh...@uwm.edu> Date: Wednesday, March 10, 2010 8:19 pm Subject: [R] help about solving two equations To: r-help@r-project.org > I have two matrix s1 and s2, each of them is 1000*1. > and I have two equations: > digamma(p)-digamma(p+q)=s1, > digamma(q)-digamma(p+q)=s2, > and I want to sovle these two equations to get the value of x and y, > which are also two 1000*1 matrices. > > I write a program like this: > > f <- function(x) { > p<- x[1]; q <- x[2]; > ((digamma(p)-digamma(p+q)-s1[2,]) )^2 > +((digamma(q)-digamma(p+q)-s2[2,]) )^2 > } > s <- 1:10/10 > g <- expand.grid(p = s, q = s) > idx <- which.min(apply(g, 1, f)) > idx > g[idx,] > > 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? > > Thanks. > > Annie > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ 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.