[R] newton method
Hi R-users, Does R has a topic on newton's method? Thank you for the info. __ 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.
Re: [R] newton method
Take a look at the functionsnlm(), optim() in the stats package and maxNR() in the maxLik package. On Mar 22, 2009, at 11:15 PM, Roslina Zakaria wrote: > Does R has a topic on newton's method? _ Professor Michael Kubovy University of Virginia Department of Psychology Postal Address: P.O.Box 400400, Charlottesville, VA 22904-4400 Express Parcels Address: Gilmer Hall, Room 102, McCormick Road, Charlottesville, VA 22903 Office:B011;Phone: +1-434-982-4729 Lab:B019; Phone: +1-434-982-4751 WWW:http://www.people.virginia.edu/~mk9y/ Skype name: polyurinsane [[alternative HTML version deleted]] __ 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.
Re: [R] newton method
Hi, you might be also interested in a general overview as given here: http://cran.r-project.org/web/views/Optimization.html Hope this helps, Roland > -Original Message- > From: r-help-boun...@r-project.org > [mailto:r-help-boun...@r-project.org] On Behalf Of Michael Kubovy > Sent: Monday, March 23, 2009 4:35 AM > To: Roslina Zakaria > Cc: r-help@r-project.org > Subject: Re: [R] newton method > > Take a look at the functionsnlm(), optim() in the stats package and > maxNR() in the maxLik package. > > On Mar 22, 2009, at 11:15 PM, Roslina Zakaria wrote: > > > Does R has a topic on newton's method? > > > _ > Professor Michael Kubovy > University of Virginia > Department of Psychology > Postal Address: > P.O.Box 400400, Charlottesville, VA 22904-4400 > Express Parcels Address: > Gilmer Hall, Room 102, McCormick Road, Charlottesville, VA 22903 > Office:B011; Phone: +1-434-982-4729 > Lab:B019; Phone: +1-434-982-4751 > WWW:http://www.people.virginia.edu/~mk9y/ > Skype name: polyurinsane > > > > > > [[alternative HTML version deleted]] > > __ > 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. > -- This mail has been sent through the MPI for Demographic Research. Should you receive a mail that is apparently from a MPI user without this text displayed, then the address has most likely been faked. If you are uncertain about the validity of this message, please check the mail header or ask your system administrator for assistance. __ 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.
Re: [R] newton method
I'm not sure what you meant by "a topic on newton's method" (algorithm? demo?), but the demonstration in the package 'animation' might help: install.packages('animation') par(pch = 20) ani.options(nmax = 50) newton.method(function(x) 5 * x^3 - 7 * x^2 - 40 * x + 100, 7.15, c(-6.2, 7.1)) Regards, Yihui -- Yihui Xie Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086 Mobile: +86-15810805877 Homepage: http://www.yihui.name School of Statistics, Room 1037, Mingde Main Building, Renmin University of China, Beijing, 100872, China On Mon, Mar 23, 2009 at 11:15 AM, Roslina Zakaria wrote: > > Hi R-users, > > Does R has a topic on newton's method? > > Thank you for the info. > __ 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.
[R] Newton method again
Hi Rolf, I would like to extend the problem that I asked you before regarding the newton method using 4 functions with 4 parameters. My functions involve the modified bessel function of the first kind which I can type them without any problem. The big problem is the Jacobian matrix. I use Maple 11.0 to get the algebraic expression for the Jacobian matrix. The problem is that the Jacobian matrix is so complicated and it is quite impossible to type the expression as I end up of more than 50 pages for the algebraic form of the Jacobian. Can you suggest any method to solve the problems? Thank you so much for your attention and help. Roslina. __ 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.
[R] Newton method iteration problem
Hi all, I am coding for finding the root of f(x)= phi(x) -alpha where phi(x) is the cumulative density function and alpha is constant . The problem right now is I can't get the "initialX" representing the root out of the while loop when ending , it seems to me it disappear when the loop ends accroding to the error message. I need help . Please suggest the cause or solution to this problem. Thanks. # code #generate target function (phi(x)-alpha) (allow input x and alpha) target<-function(x,alpha){ pnorm(x)-alpha } #generate the first derivative of the of the target function firstDerivative<-function(x){ exp(-(x^2)/2)/sqrt(2*pi) } # Finding the root by Newton method rootFinding<-function(initialX,setAlpha){ while(target(initialX,setAlpha)!=0){ initialX<-initialX-(target(initialX,setAlpha)/firstfirstDerivative(initialX) } initialX } -- View this message in context: http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13439031 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.
Re: [R] Newton method again
Hi, I do not know what your problem is, but it seems like you want to solve a system of non-linear equations. Look at the function dfsane() in the package "BB". It doesn't require you to specify jacobians. require(BB) ?dfsane 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: Roslina Zakaria Date: Tuesday, April 7, 2009 10:35 pm Subject: [R] Newton method again To: Rolf Turner , R help forum > Hi Rolf, > > I would like to extend the problem that I asked you before regarding > the newton method using 4 functions with 4 parameters. My functions > involve the modified bessel function of the first kind which I can > type them without any problem. > The big problem is the Jacobian matrix. > I use Maple 11.0 to get the algebraic expression for the Jacobian > matrix. The problem is that the Jacobian matrix is so complicated and > it is quite impossible to type the expression as I end up of more than > 50 pages for the algebraic form of the Jacobian. Can you suggest any > method to solve the problems? > > Thank you so much for your attention and help. > > > > Roslina. > > > > > __ > 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.
Re: [R] Newton method again
On 8/04/2009, at 2:31 PM, Roslina Zakaria wrote: Hi Rolf, I would like to extend the problem that I asked you before regarding the newton method using 4 functions with 4 parameters. My functions involve the modified bessel function of the first kind which I can type them without any problem. The big problem is the Jacobian matrix. I use Maple 11.0 to get the algebraic expression for the Jacobian matrix. The problem is that the Jacobian matrix is so complicated and it is quite impossible to type the expression as I end up of more than 50 pages for the algebraic form of the Jacobian. Can you suggest any method to solve the problems? I'm a frayed knot. I would suggest that if you really need to do this then you use some other numerical optimizer (optim(), nlm(), nls(), ...). These don't necessarily require an analytic expression for the Jacobian. If this is another homework exercise that *insists* that you use Newton's method, then you're in a bit of a bind. There are two possibilities that I can think of. One is to program up your own numerical Jacobian calculator, and instead of having your objective function return the analytic expression for the Jacobian at the current parameter values, have it return the numerical approximation to this Jacobian. This would be dicey; calculating such a numerical approximation requires a substantial background in numerical analysis. I wouldn't want to try this myself. The second poss. is to get Maple to return the expression for the Jacobian as Fortran code (which I *think* Maple will do --- but don't quote me on this! I'm not a Maple user.) You can then compile and dynamically load (see ?SHLIB and ?dyn.load) the Fortran code and call this inside your objective function, using .Fortran(). (See ?.Fortran .) This would be much safer than the first possibility, and is the way I would go (given that Maple will return Fortran code) but it requires you to come to terms with the use of .Fortran() which is a bit demanding of the beginner (as I gather you are). If this is indeed a homework exercise, then your instructor should be giving you (quite) a bit more guidance, IMHO. Good luck. cheers, Rolf ## Attention:\ This e-mail message is privileged and confid...{{dropped:9}} __ 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.
Re: [R] Newton method again
Hi Ravi, I did ask you some question regarding newton method sometime ago.. Now I have fixed the problem and I also wrote 2 looping code (ff1 and ff2) to evaluate the modified Bessel function of the first kind and call them in the newton code. But I dont't understand why it gives the error message but still give the result but it is incorrect(too big and too small). ff1 <- function(bb,eta,z){ r <- length(z) for (i in 1:r) { sm <- sum(besselI(z[i]*bb,eta)/besselI(z[i]*bb,eta+1))*z[i]} sm } ff1(bb,eta,z) ff2 <- function(bb,eta,z,k){ r <- length(z) for (i in 1:r) { sm1 <- sum((z[i]*bb/2)*(psigamma((0:k)+eta+1,deriv=0)/(factorial(0:k)*gamma((0:k)+eta+1 sm2 <- sum((besselI(z[i]*bb,eta)*log(z[i]*bb/2) - sm1)/besselI(z[i]*bb,eta))} sm2 } ff2(bb,eta,z,10) newton.input3 <- function(pars) { ## parameters to be approximated , note: eta <- alpha3-0.5 eta <- pars[1] bt3 <- pars[2] bt4 <- pars[3] rho <- pars[4] b1 <- (pars[2]-pars[3])^2+4*pars[2]*pars[3]*pars[4] b2 <- sqrt(b1) bb <- b2/(2*pars[2]*pars[3]*(1-pars[4])) bf2 <- (pars[3]+2*pars[2]*pars[4]-pars[2])/(2*pars[2]^2*(pars[4]-1)*b2) bf3 <- (pars[2]+2*pars[3]*pars[4]-pars[3])/(2*pars[3]^2*(pars[4]-1)*b2) bf4 <- (2*pars[2]*pars[3]*pars[4]+pars[2]^2+pars[3]^2)/(2*pars[2]*pars[3]*(pars[4]-1)^2*b2) zsm <- sum(z) psigm <- psigamma(pars[1]+0.5,deriv=0) pdz <- log(prod(z)) erh <- (1+2*pars[1])*(pars[4]-1) brh1 <- 2*pars[2]*pars[3]*pars[4]+pars[2]^2+pars[3]^2 brh2 <- 2*pars[2]*pars[3]*(pars[4]-1)^2 k <- 1000 ## function fn1a <- pdz -r*(2*psigm + log(b1))/2 fn2a <- (pars[2]*r*erh+zsm)/(2*pars[2]^2*(1-pars[4])) fn3a <- (pars[3]*r*erh+zsm)/(2*pars[3]^2*(1-pars[4])) fn4a <- (pars[2]*pars[3]*r*erh+(pars[2]+pars[3])*zsm)/(-2*pars[2]*pars[3]*(pars[4]-1)^2) ## function that involve modified Bessel function of 1st kind fn1b <- ff2(bb,eta,z,k) fn2b <- bf2*ff1(bb,eta,z) fn3b <- bf3*ff1(bb,eta,z) fn4b <- bf4*ff1(bb,eta,z) ## final function fn1 <- fn1a + fn1b fn2 <- fn2a + fn2b fn3 <- fn3a + fn3b fn4 <- fn4a + fn4b fval <- c(fn1,fn2,fn3,fn4) ## output list(fval=fval) } library(BB) start <- c(0.7,0.8,0.6,0.4) dfsane(pars=start,fn=newton.input3) newton.input3(start) > library(BB) > start <- c(0.7,0.8,0.6,0.4) > dfsane(pars=start,fn=newton.input3) Error in dfsane(pars = start, fn = newton.input3) : element 1 is empty; the part of the args list of 'length' being evaluated was: (par) > newton.input3(start) $fval [1] 103.0642 452.5835 823.6637 -1484.3209 There were 50 or more warnings (use warnings() to see the first 50) > Here is my data: > z [1] 4.2 11.2 0.8 20.4 16.6 3.8 1.2 4.0 10.8 10.2 6.6 25.6 18.2 4.6 15.0 1.2 12.0 25.4 6.4 1.6 4.8 10.0 3.0 [24] 7.0 1.8 15.0 8.6 11.2 5.4 1.8 23.2 10.8 25.4 6.0 6.0 5.0 1.4 11.0 8.4 7.4 6.4 2.6 8.6 15.8 > The answer that I should get is c(0.4960, 6.6265,2.6470,0.4378) Thank you so much for help and time. [[alternative HTML version deleted]] __ 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.
Re: [R] Newton method again
Roslina Zakaria wrote: > > Hi Ravi, > I did ask you some question regarding newton method sometime ago.. Now I > have fixed the problem and I also wrote 2 looping code (ff1 and ff2) to > evaluate the modified Bessel function of the first kind and call them in > the newton code. But I dont't understand why it gives the error message > but still give the result but it is incorrect(too big and too small). > > ff1 <- function(bb,eta,z){ > r <- length(z) > for (i in 1:r) { > sm <- sum(besselI(z[i]*bb,eta)/besselI(z[i]*bb,eta+1))*z[i]} > sm > } > ff1(bb,eta,z) > ff2 <- function(bb,eta,z,k){ > r <- length(z) > for (i in 1:r) { > sm1 <- > sum((z[i]*bb/2)*(psigamma((0:k)+eta+1,deriv=0)/(factorial(0:k)*gamma((0:k)+eta+1 > sm2 <- sum((besselI(z[i]*bb,eta)*log(z[i]*bb/2) - > sm1)/besselI(z[i]*bb,eta))} > sm2 > } > ff2(bb,eta,z,10) > > newton.input3 <- function(pars) > { ## parameters to be approximated , note: eta <- alpha3-0.5 > eta <- pars[1] > bt3 <- pars[2] > bt4 <- pars[3] > rho <- pars[4] > b1 <- (pars[2]-pars[3])^2+4*pars[2]*pars[3]*pars[4] > b2 <- sqrt(b1) > bb <- b2/(2*pars[2]*pars[3]*(1-pars[4])) > bf2 <- > (pars[3]+2*pars[2]*pars[4]-pars[2])/(2*pars[2]^2*(pars[4]-1)*b2) > bf3 <- > (pars[2]+2*pars[3]*pars[4]-pars[3])/(2*pars[3]^2*(pars[4]-1)*b2) > bf4 <- > (2*pars[2]*pars[3]*pars[4]+pars[2]^2+pars[3]^2)/(2*pars[2]*pars[3]*(pars[4]-1)^2*b2) > zsm <- sum(z) > psigm <- psigamma(pars[1]+0.5,deriv=0) > pdz <- log(prod(z)) > erh <- (1+2*pars[1])*(pars[4]-1) > brh1 <- 2*pars[2]*pars[3]*pars[4]+pars[2]^2+pars[3]^2 > brh2 <- 2*pars[2]*pars[3]*(pars[4]-1)^2 > k <- 1000 > ## function > fn1a <- pdz -r*(2*psigm + log(b1))/2 > fn2a <- (pars[2]*r*erh+zsm)/(2*pars[2]^2*(1-pars[4])) > fn3a <- (pars[3]*r*erh+zsm)/(2*pars[3]^2*(1-pars[4])) > fn4a <- > (pars[2]*pars[3]*r*erh+(pars[2]+pars[3])*zsm)/(-2*pars[2]*pars[3]*(pars[4]-1)^2) > > ## function that involve modified Bessel function of 1st kind > fn1b <- ff2(bb,eta,z,k) > fn2b <- bf2*ff1(bb,eta,z) > fn3b <- bf3*ff1(bb,eta,z) > fn4b <- bf4*ff1(bb,eta,z) > > ## final function > fn1 <- fn1a + fn1b > fn2 <- fn2a + fn2b > fn3 <- fn3a + fn3b > fn4 <- fn4a + fn4b > fval <- c(fn1,fn2,fn3,fn4) > ## output > list(fval=fval) > } > library(BB) > start <- c(0.7,0.8,0.6,0.4) > dfsane(pars=start,fn=newton.input3) > newton.input3(start) > >> library(BB) >> start <- c(0.7,0.8,0.6,0.4) >> dfsane(pars=start,fn=newton.input3) > Error in dfsane(pars = start, fn = newton.input3) : element 1 is empty; > the part of the args list of 'length' being evaluated was: > (par) >> newton.input3(start) > $fval > [1] 103.0642 452.5835 823.6637 -1484.3209 > There were 50 or more warnings (use warnings() to see the first 50) >> > > Here is my data: >> z > [1] 4.2 11.2 0.8 20.4 16.6 3.8 1.2 4.0 10.8 10.2 6.6 25.6 18.2 4.6 > 15.0 1.2 12.0 25.4 6.4 1.6 4.8 10.0 3.0 > [24] 7.0 1.8 15.0 8.6 11.2 5.4 1.8 23.2 10.8 25.4 6.0 6.0 5.0 1.4 > 11.0 8.4 7.4 6.4 2.6 8.6 15.8 >> > > You could also try package nleqslv which implements Newton and Broyden methods for solving systems of equations. I have tried to run your problem but you are not providing all the information required. Moreover your example contains errors: for example where are the arguments defined in the call of ff1 on the line "ff1(bb,eta,z)" right after the definition of ff1? Where is the variable "r" used in the lines calculating fn1a, fn2a etc. in function newton.input3? Is it the same as in ff1 and ff2? length(z)? When I insert r<-length(z) in newton.input3() I get the results shown in your post for $fval. The warnings are being given by factorial(0:k): In factorial(0:k) : value out of range in 'gammafn' Why are you assigning pars[1], pars[2] etc to scalars and then afterwards not or hardly using them? You code is inefficient since you are calling ff1 in newton.input3 three times with exactly the same input. I have tried to run your code in nleqslv but it appears to run very slowly so I can't help you any further at this point in time. What is the purpose of the loop in function ff1 for (i in 1:r) { sm <- sum(besselI(z[i]*bb,eta)/besselI(z[i]*bb,eta+1))*z[i]} (on returning from the function sm will contain the value obtained for i=r) ? Given the presentation of your problem, I cannot make head or tail of what you are trying to do so I can't help you any further. Berend Hasselman -- View this message in context: http://www.nabble.com/newton-method-tp22653758p23875467.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.
Re: [R] Newton method iteration problem
On Fri, 26 Oct 2007, kevinchang wrote: > > Hi all, > > I am coding for finding the root of f(x)= phi(x) -alpha where phi(x) is the > cumulative density function and alpha is constant . The problem right now is > I can't get the "initialX" representing the root out of the while loop when > ending , it seems to me it disappear when the loop ends accroding to the > error message. I need help . Please suggest the cause or solution to this > problem. Thanks. Learn to type without making errors? Learn to format (space and indent) your code so errors will be more obvious to you?? Your code worked for me once I corrected the typos and syntax errors. It even agrees with qnorm( setAlpha ). If all you want is root finding capability, I suggest you see ?uniroot and friends. HTH, Chuck > > # code > > #generate target function (phi(x)-alpha) (allow input x and alpha) > target<-function(x,alpha){ > pnorm(x)-alpha > } > > > #generate the first derivative of the of the target function > firstDerivative<-function(x){ > exp(-(x^2)/2)/sqrt(2*pi) > } > > # Finding the root by Newton method > rootFinding<-function(initialX,setAlpha){ > while(target(initialX,setAlpha)!=0){ > initialX<-initialX-(target(initialX,setAlpha)/firstfirstDerivative(initialX) > } > initialX > } > > > -- > View this message in context: > http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13439031 > 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. > Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ 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.
Re: [R] Newton method iteration problem
On Fri, 26 Oct 2007, Charles C. Berry wrote: > On Fri, 26 Oct 2007, kevinchang wrote: > >> >> Hi all, >> >> I am coding for finding the root of f(x)= phi(x) -alpha where phi(x) is the >> cumulative density function and alpha is constant . The problem right now is >> I can't get the "initialX" representing the root out of the while loop when >> ending , it seems to me it disappear when the loop ends accroding to the >> error message. I need help . Please suggest the cause or solution to this >> problem. Thanks. > > Learn to type without making errors? Learn to format (space and indent) > your code so errors will be more obvious to you?? > > Your code worked for me once I corrected the typos and syntax errors. > > It even agrees with qnorm( setAlpha ). Not quite. I neglected to mention that I added a 'maxIter' variable that terminates the while loop after a few dozen passes. Your loop will not terminate otherwise (unless you are uncommonly lucky). There are lots of tricks for the unwary in writing optimization software. So it is best to leave the details to the experts whenever possible. Fortunately, R has well crafted optimizers available to the user. :-) > > If all you want is root finding capability, I suggest you see > >?uniroot > > and friends. > > HTH, > > Chuck > >> >> # code >> >> #generate target function (phi(x)-alpha) (allow input x and alpha) >> target<-function(x,alpha){ >> pnorm(x)-alpha >> } >> >> >> #generate the first derivative of the of the target function >> firstDerivative<-function(x){ >> exp(-(x^2)/2)/sqrt(2*pi) >> } >> >> # Finding the root by Newton method >> rootFinding<-function(initialX,setAlpha){ >> while(target(initialX,setAlpha)!=0){ >> initialX<-initialX-(target(initialX,setAlpha)/firstfirstDerivative(initialX) >> } >> initialX >> } >> >> >> -- >> View this message in context: >> http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13439031 >> 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. >> > > Charles C. Berry(858) 534-2098 > Dept of Family/Preventive Medicine > E mailto:[EMAIL PROTECTED]UC San Diego > http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 > > __ > 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. > Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ 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.
Re: [R] Newton method iteration problem
Thanks Chales for pointing out the errors. I fixed an errorr and R accepted my "rootFinding" code. But the problem right now is that my code will work only if the intialX value is close enough to the solution. Otherwise, R says there is missing true/false value in the codition test of while loop. I can't figure out what happens. Some advice , please?? #generate target function (phi(x)-alpha) (allow input x and alpha) target<-function(x,alpha){ pnorm(x)-alpha } #generate the first derivative of the of the target function firstDerivative<-function(x){ exp(-(x^2)/2)/sqrt(2*pi) } # Finding the root by Newton method rootFinding<-function(initialX,setAlpha,maxIter){ while((target(initialX,setAlpha)!=0) && maxIter>0){ initialX<-initialX-(target(initialX,setAlpha)/firstDerivative(initialX)) maxIter<-maxIter-1 } initialX } Charles C. Berry wrote: > > On Fri, 26 Oct 2007, kevinchang wrote: > >> >> Hi all, >> >> I am coding for finding the root of f(x)= phi(x) -alpha where phi(x) is >> the >> cumulative density function and alpha is constant . The problem right now >> is >> I can't get the "initialX" representing the root out of the while loop >> when >> ending , it seems to me it disappear when the loop ends accroding to the >> error message. I need help . Please suggest the cause or solution to >> this >> problem. Thanks. > > Learn to type without making errors? Learn to format (space and indent) > your code so errors will be more obvious to you?? > > Your code worked for me once I corrected the typos and syntax errors. > > It even agrees with qnorm( setAlpha ). > > If all you want is root finding capability, I suggest you see > > ?uniroot > > and friends. > > HTH, > > Chuck > >> >> # code >> >> #generate target function (phi(x)-alpha) (allow input x and alpha) >> target<-function(x,alpha){ >> pnorm(x)-alpha >> } >> >> >> #generate the first derivative of the of the target function >> firstDerivative<-function(x){ >> exp(-(x^2)/2)/sqrt(2*pi) >> } >> >> # Finding the root by Newton method >> rootFinding<-function(initialX,setAlpha){ >> while(target(initialX,setAlpha)!=0){ >> initialX<-initialX-(target(initialX,setAlpha)/firstfirstDerivative(initialX) >> } >> initialX >> } >> >> >> -- >> View this message in context: >> http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13439031 >> 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. >> > > Charles C. Berry(858) 534-2098 > Dept of Family/Preventive > Medicine > E mailto:[EMAIL PROTECTED]UC San Diego > http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 > > __ > 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. > > -- View this message in context: http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13442728 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.
Re: [R] Newton method iteration problem
On Sat, 27 Oct 2007, kevinchang wrote: > > Thanks Chales for pointing out the errors. > I fixed an errorr and R accepted my "rootFinding" code. > But the problem right now is that my code will work only if the intialX > value is close enough to the solution. > Otherwise, R says there is missing true/false value in the codition test of > while loop. I can't figure out what happens. Some advice , please?? I think I warned you that this is a tricky business. Finding good starting values is part of the art as is dealing with the consequences of bad starting values. If you are simply trying to implement root finding, my suggestion to use uniroot() and the like still holds. If you are trying to learn about numerical optimization, then (in addition to reading a book or two on the topic), I suggest you do this: ?debug ## find out about this helpful function first debug( rootFinding ) rootFinding( 5, 0.05, 10 ) < now inspect each object just after it is created > < also check out the values of target() and firstDerivative() > When you see where this algorithm 'hit the wall', you may wish to ponder how Newton and his algorithm could have got it so wrong. ;-) I also suggest you look at the source used by qnorm() (and the other q functions). If you hunt for this you will eventually find it in /src/nmath/qnorm.c If you cannot grok what is going on from reading that source, then referring to the algorithms referenced there (AS 111 and AS 241) will likely be instructive. In the code, you will see is a comment that suggests that finding qnorm( value.near.zero.or.one ) is a bit delicate. Chuck > > #generate target function (phi(x)-alpha) (allow input x and alpha) > target<-function(x,alpha){ > pnorm(x)-alpha > } > > #generate the first derivative of the of the target function > firstDerivative<-function(x){ > exp(-(x^2)/2)/sqrt(2*pi) > } > > # Finding the root by Newton method > rootFinding<-function(initialX,setAlpha,maxIter){ > while((target(initialX,setAlpha)!=0) && maxIter>0){ > initialX<-initialX-(target(initialX,setAlpha)/firstDerivative(initialX)) > maxIter<-maxIter-1 > } > initialX > } > > > > > > > > > > > > Charles C. Berry wrote: >> >> On Fri, 26 Oct 2007, kevinchang wrote: >> >>> >>> Hi all, >>> >>> I am coding for finding the root of f(x)= phi(x) -alpha where phi(x) is >>> the >>> cumulative density function and alpha is constant . The problem right now >>> is >>> I can't get the "initialX" representing the root out of the while loop >>> when >>> ending , it seems to me it disappear when the loop ends accroding to the >>> error message. I need help . Please suggest the cause or solution to >>> this >>> problem. Thanks. >> >> Learn to type without making errors? Learn to format (space and indent) >> your code so errors will be more obvious to you?? >> >> Your code worked for me once I corrected the typos and syntax errors. >> >> It even agrees with qnorm( setAlpha ). >> >> If all you want is root finding capability, I suggest you see >> >> ?uniroot >> >> and friends. >> >> HTH, >> >> Chuck >> >>> >>> # code >>> >>> #generate target function (phi(x)-alpha) (allow input x and alpha) >>> target<-function(x,alpha){ >>> pnorm(x)-alpha >>> } >>> >>> >>> #generate the first derivative of the of the target function >>> firstDerivative<-function(x){ >>> exp(-(x^2)/2)/sqrt(2*pi) >>> } >>> >>> # Finding the root by Newton method >>> rootFinding<-function(initialX,setAlpha){ >>> while(target(initialX,setAlpha)!=0){ >>> initialX<-initialX-(target(initialX,setAlpha)/firstfirstDerivative(initialX) >>> } >>> initialX >>> } >>> >>> >>> -- >>> View this message in context: >>> http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13439031 >>> 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. >>> >> >> Charles C. Berry(858) 534-2098 >> Dept of Family/Preventive >> Medicine >> E mailto:[EMAIL PROTECTED] UC San Diego >> http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 >> >> __ >> 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. >> >> > > -- > View this message in context: > http://www.nabble.com/Newton-method-iteration-problem-tf4701085.html#a13442728 > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://sta