Hello R users, 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)) Thank You in advance for your help. Evgeniq Petrova
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
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 >stat.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide R-project.org/posting-guide.html >and provide commented, minimal, self-contained, reproducible code. >
> x^2 - y^2 = 6 > x ? y = 3 > >You can also try this # function f <- function(x) { y <- numeric(2) y[1] <- x[1]^2-x[2]^2-6 y[2] <- x[1]-x[2]-3 y } # function values transformed to scalar # minimising fnorm this way is not the best method of finding a solution for f(x)=0 # there may be values for x which minimise fnorm but do not set f(x) = 0 # but you can always try fnorm <- function(z) {p <- f(z);return(crossprod(p))} #starting values xstart <- c(0,0) # You can use nlm or nlminb nlm(fnorm,xstart) nlminb(xstart,fnorm) Sometimes minpack.lm can be used to find roots. Do library(minpack.lm) nls.lm(xstart,f) and in this case you'll see that it doesn't work. In this case stick to nlm and/or nlminb. Berend -- View this message in context: nabble.com/Non-linear-system-of-equations-tp16893056p16895541.html Sent from the R help mailing list archive at Nabble.com.
Try this:> library(Ryacas) > x <- Sym("x") > y <- Sym("y") > Solve(List(x^2+y^2==6, x-y==3), List(x,y))expression(list(list(x == root(6 - y^2, 2), y == y))) 2008/4/25 Radka Pancheva <radica at abv.bg>:> Hello R users, > > 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)) > > Thank You in advance for your help. > > Evgeniq Petrova > > ______________________________________________ > R-help at r-project.org mailing list > stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
Hello, Thank you all for your answers. I'll try them all and I'm sure it will work. > >Try this: > >> library(Ryacas) >> x <- Sym("x") >> y <- Sym("y") >> Solve(List(x^2+y^2==6, x-y==3), List(x,y)) >expression(list(list(x == root(6 - y^2, 2), y == y))) > > >2008/4/25 Radka Pancheva <radica at abv.bg>: >> Hello R users, >> >> 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)) >> >> Thank You in advance for your help. >> >> >> ______________________________________________ >> R-help at r-project.org mailing list >> stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> >