R help - Mar 2009 - newton method

Hi R-users,

Does R has a topic on newton's method?

Thank you for the info.

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

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 <xieyihui at gmail.com>
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 <zroslina at yahoo.com>
wrote:
>
> Hi R-users,
>
> Does R has a topic on newton's method?
>
> Thank you for the info.
>

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.

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