[R] newton method
Is there a non animation version of newton.method (http://bm2.genes.nig.ac.jp/RGM2/R_current/library/animation/man/newton.method.html) for finding roots of a function? It should find the roots without showing it on the GUI. thanks, Sam -- View this message in context: http://r.789695.n4.nabble.com/newton-method-tp2298100p2298100.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
saminny wrote: Is there a non animation version of newton.method (http://bm2.genes.nig.ac.jp/RGM2/R_current/library/animation/man/newton.method.html) for finding roots of a function? It should find the roots without showing it on the GUI. Have a look at package nleqslv /Berend -- View this message in context: http://r.789695.n4.nabble.com/newton-method-tp2298100p2299047.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 for single nonlinear equation
Roslina Zakaria wrote:
newton.inputsingle - function(pars,n)
{ runi - runif(974, min=0, max=1)
lendt - length(runi)
## Parameter to estimate
z - vector(length=lendt, mode= numeric)
z - pars[1]
## Constant value
alp - 2.0165 ; rho - 0.868;
c - sqrt(pi)/(gamma(alp)*(1-rho)^alp)
for (i in 1:n)
{ t1 - exp(-pars[1]/(1-rho))
t2 - (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)
bes1 - besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)
bes2 - besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)
bes3 - besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)
## Equation
f - c*t1*t2*bes1 - runi
## derivative
fprime - c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1] +
sqrt(rho)*(bes2-bes3)/(2*(1-rho)))
z[i+1] - z[i] - f/fprime
}
z
}
pars - 0.5
newton.inputsingle(pars,5)
The output :
pars - 0.5
newton.inputsingle(pars,5)
[1] 0.500 -0.4826946 -1.4653892 -2.4480838 -3.4307784 -4.4134730
Warning messages:
1: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
The warning message is the result of assigning a vector (f/fprime) to a
scalar.
You are also redefining an R built-in function: c
Furthermore it is unclear what you are trying to do.
I would advise: Keep It Simple and Stupid (KISS).
To help you along I have lifted the relevant function out of your
newton.inputsingle
(or at least what I think your function is)
zztest - function(x) {
## Constant value
alp - 2.0165 ; rho - 0.868;
czz - sqrt(pi)/(gamma(alp)*(1-rho)^alp)
t1 - exp(-x/(1-rho))
t2 - (x*(1-rho)/(2*sqrt(rho)))^(alp-0.5)
bes1 - besselI(x*sqrt(rho)/(1-rho),alp-0.5)
bes2 - besselI(x*sqrt(rho)/(1-rho),alp-1.5)
bes3 - besselI(x*sqrt(rho)/(1-rho),alp+0.5)
## Equation
f - czz*t1*t2*bes1 - .1
}
pars - 0.5
zztest(pars)
plot(zztest,0,10)
abline(0,0,lty=2)
uniroot(zztest, c(0,2))
uniroot(zztest, c(4,8))
The plot indicates how your function behaves.
The two uniroot calls solve for the two roots that appear in the plot.
You should be able to proceed on your own from here.
BTW: I do hope that this not homework.
good luck
Berend
--
View this message in context:
http://n4.nabble.com/newton-method-for-single-nonlinear-equation-tp1289991p1310861.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.
[R] Newton method
Hi r-users,
I hope somebody can help me with this code.
I would like to solve for z values using newton iteration method. I 'm not
sure which part of the code is wrong since I'm not very good at programming but
would like to learn. There seem to be some output but what I expected is a
vector of z values. Thank you so much for any help given.
newton.inputsingle - function(pars,n)
{ runi - runif(974, min=0, max=1)
lendt - length(runi)
## Parameter to estimate
z - vector(length=lendt, mode= numeric)
z - pars[1]
## Constant value
alp - 2.0165 ; rho - 0.868;
c - sqrt(pi)/(gamma(alp)*(1-rho)^alp)
for (i in 1:n)
{ t1 - exp(-pars[1]/(1-rho))
t2 - (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)
bes1 - besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)
bes2 - besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)
bes3 - besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)
## Equation
f - c*t1*t2*bes1 - runi
## derivative
fprime - c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1] +
sqrt(rho)*(bes2-bes3)/(2*(1-rho)))
z[i+1] - z[i] - f/fprime
}
z
}
pars - 0.5
newton.inputsingle(pars,5)
The output :
pars - 0.5
newton.inputsingle(pars,5)
[1] 0.500 -0.4826946 -1.4653892 -2.4480838 -3.4307784 -4.4134730
Warning messages:
1: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
2: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
3: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
4: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
5: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
[[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
your problem is because f is a vector of the same length as runi as a result, the thing you're trying to assign to z[i+1] is also a vector of that length -- View this message in context: http://n4.nabble.com/Newton-method-tp1311057p1311228.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.
[R] newton method for single nonlinear equation
Hi r-users,
I would like to solve for z values using newton iteration method. I 'm not
sure which part of the code is wrong since I'm not very good at programming but
would like to learn. There seem to be some output but what I expected is a
vector of z values. Thank you so much for any help given.
newton.inputsingle - function(pars,n)
{ runi - runif(974, min=0, max=1)
lendt - length(runi)
## Parameter to estimate
z - vector(length=lendt, mode= numeric)
z - pars[1]
## Constant value
alp - 2.0165 ; rho - 0.868;
c - sqrt(pi)/(gamma(alp)*(1-rho)^alp)
for (i in 1:n)
{ t1 - exp(-pars[1]/(1-rho))
t2 - (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)
bes1 - besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)
bes2 - besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)
bes3 - besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)
## Equation
f - c*t1*t2*bes1 - runi
## derivative
fprime - c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1] +
sqrt(rho)*(bes2-bes3)/(2*(1-rho)))
z[i+1] - z[i] - f/fprime
}
z
}
pars - 0.5
newton.inputsingle(pars,5)
The output :
pars - 0.5
newton.inputsingle(pars,5)
[1] 0.500 -0.4826946 -1.4653892 -2.4480838 -3.4307784 -4.4134730
Warning messages:
1: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
2: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
3: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
4: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
5: In z[i + 1] - z[i] - f/fprime :
number of items to replace is not a multiple of replacement length
[[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.
[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.
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 zrosl...@yahoo.com Date: Tuesday, April 7, 2009 10:35 pm Subject: [R] Newton method again To: Rolf Turner r.tur...@auckland.ac.nz, R help forum R-help@r-project.org 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
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 xieyi...@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 zrosl...@yahoo.com 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.
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.
[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 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) maxIter0){
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
qdistn functions). If you hunt for this you will eventually find it
in path to R sources/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) maxIter0){
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.
Charles C. Berry(858) 534-2098
[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 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.
