thr3ads.net - R help - [R] help about solving two equations [Mar 2010]

If this information is useful, please help other people find it:
Share via:

Shaoqiong Zhao

2010-Mar-11 01:18 UTC

[R] help about solving two equations

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

Ravi Varadhan

2010-Mar-11 03:03 UTC

head link

[R] help about solving two equations

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: rvaradhan at jhmi.edu


----- Original Message -----
From: Shaoqiong Zhao <zhaos at uwm.edu>
Date: Wednesday, March 10, 2010 8:19 pm
Subject: [R] help about solving two equations
To: r-help at 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 at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.

Berend Hasselman

2010-Mar-11 06:59 UTC

head link

[R] help about solving two equations

Shaoqiong Zhao wrote:> 
> 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.
> 
Besides BB there is another package for solving systems of equations:
nleqslv.

Use Ravi's definition of your function.
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)
library(nleqslv)

z1 <- nleqslv(xstart,fn,s=s)
z1

or for more detail 

z2 <- nleqslv(xstart,fn,s=s, control=list(trace=1))
z2

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

Berend Hasselman

2010-Mar-11 10:20 UTC

head link

[R] help about solving two equations

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.

Ravi Varadhan

2010-Mar-12 22:41 UTC

head link

[R] help about solving two equations

Here you go:


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

Hope this helps,
Ravi.

-----Original Message-----
From: Shaoqiong Zhao [mailto:zhaos at uwm.edu] 
Sent: Friday, March 12, 2010 5:20 PM
To: Ravi Varadhan
Cc: r-help at r-project.org
Subject: Re: [R] help about solving two equations

Hello Professor Ravi,

I installed the BB package. It is so useful. Thanks a lot for developing this
package for solving nonloinear equations.

I tried couple of different ways to do the loop, but still ccant work it out.

Now I have s, which is a 1000*2 matrix,
I want to solve p,q for each row of s,  

I use the following loop but it still give me one value for p,q

for (i in 1:1000)
{
fn <- function(x, s){
 f <- rep(NA, length(x))
f[1] <- digamma(x[1]) - digamma(x[1]+x[2]) - s[i,1]
 f[2] <- digamma(x[2]) - digamma(x[1]+x[2]) - s[i,2]
f
 } 
s <- matrix(
 ans <- dfsane(par=c(1,1), fn=fn, s=s)
}
 ans$par  

I still only get one value for p,q. I tried couple of different ways to do the
loop, still didnt get any results.

Sorry to bother you again.

Best,
Annie


----- Original Message -----
From: "Ravi Varadhan" <rvaradhan at jhmi.edu>
To: "Shaoqiong Zhao" <zhaos at uwm.edu>
Cc: r-help at r-project.org
Sent: Thursday, March 11, 2010 8:02:43 PM GMT -06:00 US/Canada Central
Subject: Re: [R] help about solving two equations

You have to install it from CRAN before you can load it into your session.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Shaoqiong Zhao <zhaos at uwm.edu>
Date: Thursday, March 11, 2010 9:00 pm
Subject: Re: [R] help about solving two equations
To: Ravi Varadhan <rvaradhan at jhmi.edu>
Cc: r-help at r-project.org

> Hello Professor Ravi,
>  
>  I tried to load BB into R, but I got the following message:
>   library("BB")
>  Error in library("BB") : there is no package called 'BB'
>  > library(BB)
>  Error in library(BB) : there is no package called 'BB'
>  
>  Can you tell me why?
>  
>  Thanks a lot.
>  
>  Annie
>  
>  
>  ----- Original Message -----
>  From: "Ravi Varadhan" <rvaradhan at jhmi.edu>
>  To: "Shaoqiong Zhao" <zhaos at uwm.edu>
>  Cc: r-help at r-project.org
>  Sent: Wednesday, March 10, 2010 9:03:25 PM GMT -06:00 US/Canada Central
>  Subject: Re: [R] help about solving two equations
>  
>  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: rvaradhan at jhmi.edu
>  
>  
>  ----- Original Message -----
>  From: Shaoqiong Zhao <zhaos at uwm.edu>
>  Date: Wednesday, March 10, 2010 8:19 pm
>  Subject: [R] help about solving two equations
>  To: r-help at 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 at r-project.org mailing list
>  >  
>  >  PLEASE do read the posting guide 
>  >  and provide commented, minimal, self-contained, reproducible code. 
>

Maybe Matching Threads

Search for more apparently analagous threads

R help - Mar 2010 - help about solving two equations

[R] help about solving two equations

[R] help about solving two equations

[R] help about solving two equations

[R] help about solving two equations

[R] help about solving two equations

Maybe Matching Threads