Quantcast

newton method

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

newton method

Roslina Zakaria

Hi R-users,

Does R has a topic on newton's method?

Thank you for the info.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: newton method

Michael Kubovy
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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: newton method

Rau, Roland
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: [hidden email]
> [mailto:[hidden email]] On Behalf Of Michael Kubovy
> Sent: Monday, March 23, 2009 4:35 AM
> To: Roslina Zakaria
> Cc: [hidden email]
> 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]]
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: newton method

Yihui Xie
In reply to this post by Roslina Zakaria
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 <[hidden email]>
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 <[hidden email]> wrote:
>
> Hi R-users,
>
> Does R has a topic on newton's method?
>
> Thank you for the info.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Newton method again

Roslina Zakaria
In reply to this post by Roslina Zakaria

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.




______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newton method again

Ravi Varadhan
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: [hidden email]


----- Original Message -----
From: Roslina Zakaria <[hidden email]>
Date: Tuesday, April 7, 2009 10:35 pm
Subject: [R] Newton method again
To: Rolf Turner <[hidden email]>, R help forum <[hidden email]>


>  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.
>  
>  
>  
>  
>  ______________________________________________
>  [hidden email] mailing list
>  
>  PLEASE do read the posting guide
>  and provide commented, minimal, self-contained, reproducible code.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newton method again

Rolf Turner
In reply to this post by Roslina Zakaria

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}}

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newton method again

Roslina Zakaria
In reply to this post by Ravi Varadhan
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]]


______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newton method again

Berend Hasselman
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


Loading...