Re: [R] Non-linear system of equations

From: Ravi Varadhan <rvaradhan_at_jhmi.edu>
Date: Fri, 25 Apr 2008 09:35:54 -0400

Hi Radka,

The problem lies in your function my.mm(). Especially, the following lines:

p <- .Machine$double.eps
q <- .Machine$double.eps

Even if you take out these 2 lines, still the definition is not correct, I think, since there is a trivial answer x* = c(0, q).

I am not sure how you intended to use these, but they are the problem. Therefore, you need to define your function correctly.

Ravi.



Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan_at_jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html  



-----Original Message-----
From: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-project.org] On Behalf Of Radka Pancheva
Sent: Friday, April 25, 2008 8:13 AM
To: R-help_at_r-project.org
Subject: Re: [R] Non-linear system of equations

Hello Paul,

Thank you for your quick answer. I have tried to use your advice and to estimate the parameters of beta distribution with moments matching. This is my code:

ex <- 0.3914877
ex2 <- 0.2671597

my.mm <- function(x){

	p <- x[1]
	q <- x[2]
	p <-  .Machine$double.eps
	q <-  .Machine$double.eps

	F <- rep(NA,2)

	F[1] <- p/(p + q)
	F[2]<- (p*q + (p + q + 1)*p^2)/((p + q + 1)*(p + q)^2)

	return(F)

}

p0 <- c(ex,ex2)

dfsane(par=p0, fn=my.mm,control=list(maxit=50000))

and I became the following output:

.
iteration: 3640 ||F(xn)|| = 0.7071068 iteration: 3641 ||F(xn)|| = 0. 7071068 .
iteration: 49990 ||F(xn)|| = 0. 7071068 iteration: 50000 ||F(xn)|| = 0. 7071068
$par

[1] -446.2791 -446.4034

$residual

[1] 0.5

$fn.reduction

[1] 0

$feval

[1] 828495

$iter

[1] 50001

$convergence

[1] 1

$message

[1] "Maximum limit for iterations exceeded"

I have tried maxiter=100000 but the output is the same. I know that ex and ex2 are bringing the problems, but I am stuck with them. How can I make it convergent?

Thank you,

Evgeniq

 >2008/4/25 Radka Pancheva <radica_at_abv.bg>:  >> I am trying to estimate the parameters of a bimodal normal distribution using moments matching, so I have to solve a non-linear system of equations. How can I solve the following simple example?

 >>
 >>  x^2 - y^2 = 6
 >>  x - y = 3
 >>
 >>  I heard about nlsystemfit, but I don't know how to run it exactly. I
have tried the following code, but it doesn't really work:
 >>
 >>
 >>  f1 <-y~ x[1]^2-x[2]^2-6
 >>  f2 <-z~ x[1]-x[2]-3
 >>  f  <- list(f1=0,f2=0)
 >>  nlsystemfit("OLS",f,startvals=c(0,0))
 >
 >You could try the recent package BB by Ravi Varadhan. The code could
 >be the following:
 >
 >library(BB)
 >
 >f <- function(x) {
 >  x1 <- x[1]
 >  x2 <- x[2]
 >
 >  F <- rep(NA, 2)
 >
 >  F[1] <- x1^2 - x2^2 - 6
 >  F[2] <- x1 - x2 - 3
 >
 >  return(F)
 >}
 >
 >p0 <- c(1,2)
 >dfsane(par=p0, fn=f,control=list(maxit=3000))
 >
 >I got the solution:
 >
 >x1 = 2.5
 >x2 = -0.5
 >
 >Paul
 >
 >______________________________________________
 >R-help_at_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.  >

R-help_at_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.

R-help_at_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. Received on Fri 25 Apr 2008 - 13:51:56 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 25 Apr 2008 - 14:30:32 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive