Re: [R] uniroot
Hello,
Yes, it's FAQ 7.31 but it's not uniroot's fault.
The proper way of checking the result would be to call the function fun,
not to take the digits output by the print method and compute the
function's expression with them.
rui@rui:~$ R -q -f rhelp.R
fun <- function(x) {x^x -23}
# Clearly the root lies somewhere between 2.75 and 3.00
x0 <- uniroot(fun, lower = 2.75, upper = 3.00, tol = 0.001)
# uniroot result
x0$f.root
#[1] 0.0001136763
# check the root, right
fun(x0$root)
#[1] 0.0001136763
# OP result, wrong
2.923125^2.923125 - 23
#[1] 0.000125
sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=pt_PT.UTF-8 LC_NUMERIC=C
[3] LC_TIME=pt_PT.UTF-8LC_COLLATE=pt_PT.UTF-8
[5] LC_MONETARY=pt_PT.UTF-8LC_MESSAGES=pt_PT.UTF-8
[7] LC_PAPER=pt_PT.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=pt_PT.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_4.1.1
Also, why change the default tol to a lesser one?
Hope this helps,
Rui Barradas
Às 18:42 de 27/08/21, Jeff Newmiller escreveu:
Yes. This kind of issue is covered in any decent undergraduate course in
numerical methods... it is not specific to R. It is also related to FAQ 7.31.
https://en.m.wikipedia.org/wiki/Root-finding_algorithms
https://en.m.wikipedia.org/wiki/Floating-point_arithmetic#Representable_numbers,_conversion_and_rounding
On August 27, 2021 10:30:38 AM PDT, Thomas Subia via R-help
wrote:
Colleagues,
I've been using uniroot to identify a root of an equation.
As a check, I always verify that calculated root.
This is where I need some help.
Consider the following script
fun <- function(x) {x^x -23}
# Clearly the root lies somewhere between 2.75 and 3.00
uniroot(fun, lower = 2.75, upper = 3.00, tol = 0.001)
# output
$root
[1] 2.923125
$f.root
[1] 0.0001136763
# Let's verify this root.
2.923125^2.923125 - 23
0.000125
This result is different than what was calculated with uniroot
0.000125# verified check using x = 2.923125
0.0001136763# using $f.root
Does this imply that the root output of 2.923125 may need more significant
digits displayed?
I suspect that whatever root is calculated, that root may well be dependent
on what interval one defines where the root may occur
and what tolerance one has input.
I am not sure that is the case, nevertheless, it's worth asking the
question.
Some guidance would be appreciated.
Thanks!
Thomas Subia
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] uniroot
Yes. This kind of issue is covered in any decent undergraduate course in
numerical methods... it is not specific to R. It is also related to FAQ 7.31.
https://en.m.wikipedia.org/wiki/Root-finding_algorithms
https://en.m.wikipedia.org/wiki/Floating-point_arithmetic#Representable_numbers,_conversion_and_rounding
On August 27, 2021 10:30:38 AM PDT, Thomas Subia via R-help
wrote:
>Colleagues,
>
>I've been using uniroot to identify a root of an equation.
>As a check, I always verify that calculated root.
>This is where I need some help.
>
>Consider the following script
>
>fun <- function(x) {x^x -23}
>
># Clearly the root lies somewhere between 2.75 and 3.00
>
>uniroot(fun, lower = 2.75, upper = 3.00, tol = 0.001)
>
># output
>$root
>[1] 2.923125
>
>$f.root
>[1] 0.0001136763
>
># Let's verify this root.
>
>2.923125^2.923125 - 23
>
>0.000125
>
>This result is different than what was calculated with uniroot
>0.000125 # verified check using x = 2.923125
>0.0001136763 # using $f.root
>
>Does this imply that the root output of 2.923125 may need more significant
>digits displayed?
>
>I suspect that whatever root is calculated, that root may well be dependent
>on what interval one defines where the root may occur
>and what tolerance one has input.
>I am not sure that is the case, nevertheless, it's worth asking the
>question.
>
>Some guidance would be appreciated.
>
>Thanks!
>
>Thomas Subia
>
>__
>R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>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.
--
Sent from my phone. Please excuse my brevity.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] uniroot
Colleagues,
I've been using uniroot to identify a root of an equation.
As a check, I always verify that calculated root.
This is where I need some help.
Consider the following script
fun <- function(x) {x^x -23}
# Clearly the root lies somewhere between 2.75 and 3.00
uniroot(fun, lower = 2.75, upper = 3.00, tol = 0.001)
# output
$root
[1] 2.923125
$f.root
[1] 0.0001136763
# Let's verify this root.
2.923125^2.923125 - 23
0.000125
This result is different than what was calculated with uniroot
0.000125# verified check using x = 2.923125
0.0001136763# using $f.root
Does this imply that the root output of 2.923125 may need more significant
digits displayed?
I suspect that whatever root is calculated, that root may well be dependent
on what interval one defines where the root may occur
and what tolerance one has input.
I am not sure that is the case, nevertheless, it's worth asking the
question.
Some guidance would be appreciated.
Thanks!
Thomas Subia
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] uniroot problem
Thanks a lot - f was renamed FF and things are OK BW Troels -Oprindelig meddelelse- Fra: Berwin A Turlach Sendt: 19. december 2018 10:27 Til: Troels Ring Emne: Re: [R] uniroot problem G'day Troels, On Wed, 19 Dec 2018 10:03:09 +0100 "Troels Ring" wrote: > Dear friends and helpers - in the script below, uniroot is called with > a function CHB that calls a function, Charge. On its own, CHB > apparently does what is expected, but from within uniroot, problems > appear. An error is thrown > > Error in f(lower, ...) : could not find function "f" > > So CHB is not seen or understood from within uniroot? Read the help page of uniroot. The first argument is called "f", it is the function for which the root is searched. In your call: > uniroot(CHB,interval=c(1e-19,5),tol=.Machine$double.eps,maxiter=10, > Na=Na,Cl=Cl,K=K,TOT=TOT,f=f)$root You pass the object "f" (a vector "f <- c(1,1,1)" created earlier in your code) to the argument "f". Presumably there is no function called "f" in your search path, and so R correctly complains that it cannot find the function whose root you are looking for. In R, arguments are passed first by exact matching of actual and formal arguments, then by partial matching and then by position. The easiest fix is probably to rename the object "f" to something else. Cheers, Berwin __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] uniroot problem
Dear friends and helpers - in the script below, uniroot is called with a
function CHB that calls a function, Charge. On its own, CHB apparently does
what is expected, but from within uniroot, problems appear. An error is
thrown
Error in f(lower, ...) : could not find function "f"
So CHB is not seen or understood from within uniroot?
I'm on windows 10, 64 bit R version 3.5.1 (2018-07-02)
All best wishes
Troels
kw <- 1e-14
TOT <- 1
Pk1 <- 10^-2.16
Pk2 <- 10^-7.21
Pk3 <- 10^-12.32
K <- c(Pk1,Pk2,Pk3)
f <- c(1,1,1)
H <- 10^-7.4
Charge <- function(TOT,f,K,H)
{
num <- c()
num[1] <- K[1]/(f[1]^2*H)
for (i in 2:length(K)) num[i] <- i*prod(K[1:i])/(f[1]^i*f[i]*H^i)
num <- sum(num)
denum <- c()
denum[1] <- 1+ K[1]/(f[1]^2*H)
for (i in 2:length(K)) denum[i] <- prod(K[1:i])/(f[1]^i*f[i]*H^i)
denum <- sum(denum)
num/denum
}
Na <- 0.140
Cl <- 0.1
CHB <- function(Na,Cl,H,K,f,TOT) {Na-Cl+H-kw/(f[1]^2*H)-Charge(TOT,f,K,H)}
H <- uniroot(CHB,interval=c(1e-19,5),tol=.Machine$double.eps,maxiter=10,
Na=Na,Cl=Cl,K=K,TOT=TOT,f=f)$root
#Error in f(lower, ...) : could not find function "f"
CHB(Na,Cl,10^-7.4,K,f,TOT) # -1.567668 OK right!
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Uniroot error message with in intergration
Many thanks Ellison
I have modified it as you suggested but I have this error message Error
in f(lower, ...) : unused argument(s) (N = 54)
I am not sure which arguments I have missed?
*y - function(t,n){
diff - 0.5
df1 - 2*n-2
ncp1 - sqrt((diff^2*n)/2)
p - 1- pt(t,df=df1)
test - qt((1-p),df=df1,ncp=ncp1)*(1/sqrt(2))
return(test)
}
integ - function(n){
1-integrate(y,lower=0,upper=2.7,n)$value -0.8
}
uniroot(integ,lower=0,upper=1000,N=n)
traceback()*
--
View this message in context:
http://r.789695.n4.nabble.com/Uniroot-error-message-with-in-intergration-tp4634247p4634278.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] Uniroot error message with in intergration
Dear all
I am trying to calculate the value of n using uniroot. Here is the message
I am having:
Error in uniroot(integ, lower = 0, upper = 1000, n) :
'interval' must be a vector of length 2
Please would you be able to give me an indication on why I am having this
error message.
Many thanks
EXAMPLE BELOW:
## t = statistics test from t -distribution [-Inf,Inf]
## df == degree of freedom
###p - p value
diff - 0.5
y - function(t,n){
df - 2*n-2
ncp1 - sqrt((diff^2*n)/2)
p - 1- pt(t,df=df1)
test - qt((1-p),df=df1,ncp=ncp1)*(1/sqrt(2))
return(test)
}
integ - function(t,n){
1-integrate(y,lower=0,upper=3.6,n)$value -0.5
}
uniroot(integ,lower=0,upper=1000)
--
View this message in context:
http://r.789695.n4.nabble.com/Uniroot-error-message-with-in-intergration-tp4634247.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] Uniroot error message with in intergration
On Fri, Jun 22, 2012 at 4:00 PM, jeka12386 jeka12...@hotmail.co.uk wrote:
Dear all
I am trying to calculate the value of n using uniroot. Here is the message
I am having:
Error in uniroot(integ, lower = 0, upper = 1000, n) :
'interval' must be a vector of length 2
Please would you be able to give me an indication on why I am having this
error message.
Many thanks
EXAMPLE BELOW:
## t = statistics test from t -distribution [-Inf,Inf]
## df == degree of freedom
###p - p value
diff - 0.5
y - function(t,n){
df - 2*n-2
ncp1 - sqrt((diff^2*n)/2)
p - 1- pt(t,df=df1)
test - qt((1-p),df=df1,ncp=ncp1)*(1/sqrt(2))
return(test)
}
integ - function(t,n){
1-integrate(y,lower=0,upper=3.6,n)$value -0.5
}
uniroot(integ,lower=0,upper=1000)
Running this I get a different error: Error in f(x, ...) : argument
n is missing, with no default
Perhaps you need to define n somewhere? It's also not clear to me what
happens to the t argument of integ? Is it supposed to be the y of
integrate?
Best,
Michael
--
View this message in context:
http://r.789695.n4.nabble.com/Uniroot-error-message-with-in-intergration-tp4634247.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-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] Uniroot error message with in intergration
I have defined t at the beginning of my query. I have added n on uniroot below and still getting the same error message uniroot(integ,lower=0,upper=1000,n) -- View this message in context: http://r.789695.n4.nabble.com/Uniroot-error-message-with-in-intergration-tp4634247p4634264.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] Uniroot error message with in intergration
Error in uniroot(integ, lower = 0, upper = 1000, n) :
'interval' must be a vector of length 2
Please would you be able to give me an indication on why I am having this
error message.
Because uniroot has a second parameter called 'interval' which overrides lower
and upper and you have given uniroot a scalar called n as the first non-named
argument. R has used its usual positional matching rules, which say that Any
unmatched formal arguments are bound to unnamed supplied arguments, in order.
your first unmatched formal argument (interval) has been matched to your first
unnamed argument (n) so uniroot has tried to use n as your interval.
You need to add your own parameters as named parameters; eg if your function
took something called N, you would specify N=n
S Ellison
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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] Uniroot error
Dear All
I am trying to find a uniroot of a function within another function (see
example) but I am getting an error message (f()values at end points not of
opposite sign). I was wondering if you would be able to advise how redefine
my function so that I can find the solution. In short my first function
calculates the intergrale which is function of t , I need to find the
uniroot of n defined in the second function.
y - function(t){
(dnorm(t,mean=(diff*sqrt(n/2)),sd=sqrt(rho)))*(pnorm((qnorm((1-alpha),mean=0,sd=1)-t)/(sqrt(1-rho^2
}
inter - function(n){
integrate(y,lower=-Inf,upper=Inf)$value-0.8
}
rho - 0.5
alpha - 0.0125
diff - 0.5
n1 - uniroot(inter,lower=1,upper=10)$root
--
View this message in context:
http://r.789695.n4.nabble.com/Uniroot-error-tp4539816p4539816.html
Sent from the R help mailing list archive at Nabble.com.
[[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] Uniroot error
On 07-04-2012, at 19:48, jeka12386 wrote:
Dear All
I am trying to find a uniroot of a function within another function (see
example) but I am getting an error message (f()values at end points not of
opposite sign). I was wondering if you would be able to advise how redefine
my function so that I can find the solution. In short my first function
calculates the intergrale which is function of t , I need to find the
uniroot of n defined in the second function.
y - function(t){
(dnorm(t,mean=(diff*sqrt(n/2)),sd=sqrt(rho)))*(pnorm((qnorm((1-alpha),mean=0,sd=1)-t)/(sqrt(1-rho^2
}
inter - function(n){
integrate(y,lower=-Inf,upper=Inf)$value-0.8
}
rho - 0.5
alpha - 0.0125
diff - 0.5
n1 - uniroot(inter,lower=1,upper=10)$root
When I run your example as presented by you, R issues an error message:
Error in dnorm(t, mean = (diff * sqrt(n/2)), sd = sqrt(rho))
object 'n' not found
Calls: uniroot - f - integrate - Anonymous - f - dnorm
Execution halted
Indeed, the argument n of function inter is not passed to function y, which is
using n (and is not global).
Your example I changed to this:
y - function(t, n) {
(dnorm(t, mean=(diff*sqrt(n/2)),
sd=sqrt(rho)))*(pnorm((qnorm((1-alpha),mean=0,sd=1)-t)/(sqrt(1-rho^2
}
inter - function(n) {
integrate(y,lower=-Inf,upper=Inf, n=n)$value-0.8
}
rho - 0.5
alpha - 0.0125
diff - 0.5
n1 - uniroot(inter,lower=1,upper=10)
n1
Changes: n added to arguments of y and argument n of function inter passed to y.
Output:
$root
[1] 9.210123
$f.root
[1] -6.86966e-11
$iter
[1] 9
$estim.prec
[1] 6.103516e-05
Finally: do realize that argument t of function y will be a vector (that's what
integrate does).
Berend
__
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] uniroot function question
I have one equation, two unknowns, so I am trying to build the solution set
by running through possible values for one unknown, and then using uniroot
to solve for the accompanying second solution, then graphing the two
vectors.
p0 = .36
f = function(x) 0.29 * exp(5.66*(x - p0))
f.integral = integrate(f, p0, 1)
p1 = p0 + .01
i = 1
n = (1 - p0)/.01
p1.vector = rep(0,n)
p2.vector = rep(0,n)
for (i in 1:n) {
p1.vector[i] = p1
fcn = function(p2) p1*f(p1) + (.20/5.66)*(exp(5.66*(p2 - p0)) -
exp(5.66*(p1 - p0))) + (1 - p2)*f(p2) - as.numeric(f.integral$value)
sol = uniroot(try, lower = p1, upper = 1)
p2.vector[i] = p2
i = i+1
p1 = p1 + .01
}
plot(p1.vector,p2.vector)
p1, p2 both have to be between p0 and 1, p1 p2. Is there a better way to
do this? I keep getting the error that my lower and upper bounds are not of
opposite sign, but I don't know how to find the correct interval values in
that case. This may not even be a uniroot question (although I don't know
how to find those values in a general sense), if there is a better way to do
the same thing.
Any ideas?
--
View this message in context:
http://r.789695.n4.nabble.com/uniroot-function-question-tp4195737p4195737.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] uniroot function question
kchkchkch wrote
I have one equation, two unknowns, so I am trying to build the solution
set by running through possible values for one unknown, and then using
uniroot to solve for the accompanying second solution, then graphing the
two vectors.
p0 = .36
f = function(x) 0.29 * exp(5.66*(x - p0))
f.integral = integrate(f, p0, 1)
p1 = p0 + .01
i = 1
n = (1 - p0)/.01
p1.vector = rep(0,n)
p2.vector = rep(0,n)
for (i in 1:n) {
p1.vector[i] = p1
fcn = function(p2) p1*f(p1) + (.20/5.66)*(exp(5.66*(p2 - p0)) -
exp(5.66*(p1 - p0))) + (1 - p2)*f(p2) - as.numeric(f.integral$value)
sol = uniroot(try, lower = p1, upper = 1)
p2.vector[i] = p2
i = i+1
p1 = p1 + .01
}
plot(p1.vector,p2.vector)
p1, p2 both have to be between p0 and 1, p1 p2. Is there a better way
to do this? I keep getting the error that my lower and upper bounds are
not of opposite sign, but I don't know how to find the correct interval
values in that case. This may not even be a uniroot question (although I
don't know how to find those values in a general sense), if there is a
better way to do the same thing.
Shouldn't the line with uniroot be (why try as function?)
sol = uniroot(fcn, lower = p1, upper = 1)
Before the line with for(i in 1:n) insert the following
fcn = function(p2) p1*f(p1) + (.20/5.66)*(exp(5.66*(p2 - p0)) - exp(5.66*(p1
- p0))) + (1 - p2)*f(p2) - as.numeric(f.integral$value)
curve(fcn, from=0, to=1)
And it should be obvious what the problem is.
Berend
--
View this message in context:
http://r.789695.n4.nabble.com/uniroot-function-question-tp4195737p4196220.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] uniroot function question
Heh. Yes, Berend Hasselman wrote Shouldn't the line with uniroot be (why try as function?) sol = uniroot(fcn, lower = p1, upper = 1) Before the line with for(i in 1:n) insert the following fcn = function(p2) p1*f(p1) + (.20/5.66)*(exp(5.66*(p2 - p0)) - exp(5.66*(p1 - p0))) + (1 - p2)*f(p2) - as.numeric(f.integral$value) curve(fcn, from=0, to=1) And it should be obvious what the problem is. Berend I had been trying various solutions (so I'd been trying stating the functions a bunch of different ways)--hence the try in there. but then I fixed it as you suggested--name the function correctly (insert eyeroll at myself here), and move the fcn = line to outside the for loop--and I still get the same error. -- View this message in context: http://r.789695.n4.nabble.com/uniroot-function-question-tp4195737p4196520.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] uniroot function question
kchkchkch wrote Heh. Yes, Berend Hasselman wrote Shouldn't the line with uniroot be (why try as function?) sol = uniroot(fcn, lower = p1, upper = 1) Before the line with for(i in 1:n) insert the following fcn = function(p2) p1*f(p1) + (.20/5.66)*(exp(5.66*(p2 - p0)) - exp(5.66*(p1 - p0))) + (1 - p2)*f(p2) - as.numeric(f.integral$value) curve(fcn, from=0, to=1) And it should be obvious what the problem is. Berend I had been trying various solutions (so I'd been trying stating the functions a bunch of different ways)--hence the try in there. but then I fixed it as you suggested--name the function correctly (insert eyeroll at myself here), and move the fcn = line to outside the for loop--and I still get the same error. Yes of course. But you are misinterpreting what I wrote. I said to insert two lines (NOT move) before the start of the for loop. The second line to insert was curve() which will draw a plot of your function for the first case. And you will see that the function fcn does not cross the horizontal axis which it should do if it is to have a root for fcn(x)=0. So you need to take a closer look at your function. Berend -- View this message in context: http://r.789695.n4.nabble.com/uniroot-function-question-tp4195737p4198877.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] Uniroot - error
Hi,
I have tried to use uniroot to solve a value (value a in my function) that
gives f=0, and I repeat this process for 1 times(stimulations). However
error occures from the 4625th stimulation - Error in uniroot(f, c(0, 2),
maxiter = 1000, tol = 0.001) :
f() values at end points not of opposite sign
I have also tried interval of (lower=min(U), upper=max(U)) and it won't work
as well.
Can anyone help me as I have struggled for few days already and I have to
finish it soon. Thanks.
numsim=1
set.seed(12345)
P = c()
for (m in 1:numsim) {
Y = rnorm(140,0.0125,(0.005^(1/2)))
U = exp(X1)
.
.(sorry i have to skip the code in between otherwise
..my assignment will get penalty for plagarism according to
those screening sotware)
S = sum(.)
f = function(a){sum(F*(answer^(910:1))) - S}
answer = uniroot(f, c(0,2), maxiter=1000,tol=0.001)$root
P[m] = answer^26 - 1
}
all the vectors are correct; it works without stimulation; it also works for
loop(1:4624) but after 4625 there's error.
--
View this message in context:
http://r.789695.n4.nabble.com/Uniroot-error-tp3502628p3502628.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] Uniroot - error
You should ask your instructor or teaching assistant for help. R-help
is not for doing homework.
Duncan Murdoch
On 06/05/2011 9:00 AM, CarJabo wrote:
Hi,
I have tried to use uniroot to solve a value (value a in my function) that
gives f=0, and I repeat this process for 1 times(stimulations). However
error occures from the 4625th stimulation - Error in uniroot(f, c(0, 2),
maxiter = 1000, tol = 0.001) :
f() values at end points not of opposite sign
I have also tried interval of (lower=min(U), upper=max(U)) and it won't work
as well.
Can anyone help me as I have struggled for few days already and I have to
finish it soon. Thanks.
numsim=1
set.seed(12345)
P = c()
for (m in 1:numsim) {
Y = rnorm(140,0.0125,(0.005^(1/2)))
U = exp(X1)
.
.(sorry i have to skip the code in between otherwise
..my assignment will get penalty for plagarism according to
those screening sotware)
S = sum(.)
f = function(a){sum(F*(answer^(910:1))) - S}
answer = uniroot(f, c(0,2), maxiter=1000,tol=0.001)$root
P[m] = answer^26 - 1
}
all the vectors are correct; it works without stimulation; it also works for
loop(1:4624) but after 4625 there's error.
--
View this message in context:
http://r.789695.n4.nabble.com/Uniroot-error-tp3502628p3502628.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-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] Uniroot - error
sorry I am not asking someone to do my homework, as I have finished all the procedure. I am just wondering why this technical error occurs, so I can fix it myself. By the way i don't have any instructor or teaching assistant for help, so any suggestion for the error will be appreciated. Thanks very much. -- View this message in context: http://r.789695.n4.nabble.com/Uniroot-error-tp3502628p3502773.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] Uniroot - error
On Fri, May 6, 2011 at 6:33 AM, CarJabo carly_...@hellokitty.com wrote: sorry I am not asking someone to do my homework, as I have finished all the procedure. I am just wondering why this technical error occurs, so I can fix it myself. My guess would be it has something to do with the random data generated at the 4625th simulation, but you have not posted a reproducible example (as requested in the posting guide) so it is not really possible to say. By the way i don't have any instructor or teaching assistant for help, so any suggestion for the error will be appreciated. If whatever institution you are taking this at simply gives assignments, grades them and penalizes for plagiarism without having any instructors or teachers, I recommend moving to an institution where classes are taught by someone who can answer questions etc. As Dr. Murdoc (and the posting guide) said and say, R-help is not for homework. Thanks very much. Good luck, Josh -- View this message in context: http://r.789695.n4.nabble.com/Uniroot-error-tp3502628p3502773.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. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] uniroot speed and vectorization?
curiosity---given that vector operations are so much faster than scalar operations, would it make sense to make uniroot vectorized? if I read the uniroot docs correctly, uniroot() calls an external C routine which seems to be a scalar function. that must be slow. I am thinking a vectorized version would be useful for an example such as of - function(x,a) ( log(x)+x+a ) uniroot( of, c( 1e-7, 100 ), a=rnorm(100) ) I would have timed this, but I would have used a 'for' loop, which is probably not the R way of doing this. has someone already written a package that does this? /iaw Ivo Welch (ivo.we...@brown.edu, ivo.we...@gmail.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] uniroot speed and vectorization?
On Sat, Apr 02, 2011 at 08:24:07AM -0400, ivo welch wrote:
curiosity---given that vector operations are so much faster than
scalar operations, would it make sense to make uniroot vectorized? if
I read the uniroot docs correctly, uniroot() calls an external C
routine which seems to be a scalar function. that must be slow. I am
thinking a vectorized version would be useful for an example such as
of - function(x,a) ( log(x)+x+a )
uniroot( of, c( 1e-7, 100 ), a=rnorm(100) )
Hi.
The slowest part of the solution using uniroot() is the repeated
evaluation of the R level input function. If this can be vectorized,
then a faster algorithm could be possible. The following is an
example of a vectorized bisection, which is simpler and less
efficient than zeroin used in uniroot().
of - function(x,a) { log(x)+x+a }
a - rnorm(1000)
x1 - rep(1e-7, times=length(a))
x2 - rep(100, times=length(a))
stopifnot(of(x1, a) 0)
stopifnot(of(x2, a) 0)
for (i in 1:60) {
x3 - (x1 + x2)/2
pos - of(x3, a) 0
y1 - ifelse(pos, x1, x3)
y2 - ifelse(pos, x3, x2)
x1 - y1
x2 - y2
}
print(range(of(x1, a)))
print(range(of(x2, a)))
It can be more efficient to find approximations of the roots
using a few iterations of the above approach and then switch
to the Newton method, which can be vectorized easily.
Hope this helps for a start.
Petr Savicky.
__
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] uniroot
Thanks for your advice. There was an error in the equation that is was copying. Doug -- View this message in context: http://r.789695.n4.nabble.com/uniroot-tp3260090p3264288.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] uniroot
Hi,
I am using the uniroot function in order to carry out a bivariate Monte
Carlo simulation using the logistics model.
I have defined the function as:
BV.FV - function(x,y,a,A)
(((x^(-a^-1)+y^(-a^-1))^(a-1))*(y^(a-1/a))*(exp(-((1^(-a^-1)+y^(-a^-1))^a)+y^-1)))-A
and the procedure is as follows:
Randomly generate values of A~(0,1), y0 = -(lnA)^-1
Where: A=Pr{Xxi|Y=yi-1} and a is the dependency between X and y (0.703)
Use y0 to determine x where x = x(i) and y = y0(i-1) in the above equation.
Unfortunately when the randomly defined A gets to approximately 0.46 the
intervals I provide no longer have the opposite signs (both -ve).
I have tried various different upper limits but I still can't seem to find a
root.
Does anyone have any suggestions?
Thanks,
Doug
--
View this message in context:
http://r.789695.n4.nabble.com/uniroot-tp3260090p3260090.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] uniroot
dpender wrote:
Hi,
I am using the uniroot function in order to carry out a bivariate Monte
Carlo simulation using the logistics model.
I have defined the function as:
BV.FV - function(x,y,a,A)
(((x^(-a^-1)+y^(-a^-1))^(a-1))*(y^(a-1/a))*(exp(-((1^(-a^-1)+y^(-a^-1))^a)+y^-1)))-A
and the procedure is as follows:
Randomly generate values of A~(0,1), y0 = -(lnA)^-1
Where: A=Pr{Xxi|Y=yi-1} and a is the dependency between X and y (0.703)
Use y0 to determine x where x = x(i) and y = y0(i-1) in the above
equation.
Unfortunately when the randomly defined A gets to approximately 0.46 the
intervals I provide no longer have the opposite signs (both -ve).
Try
curve(BV.FV(x,y=y0,a=.703,A=.7),from=1,to=10)
for different values of a, A and y0.
I have a feeling that your function or something else is not quite correct.
Berend
--
View this message in context:
http://r.789695.n4.nabble.com/uniroot-tp3260090p3260438.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] uniroot
On Fri, Feb 04, 2011 at 04:35:00AM -0800, dpender wrote:
Hi,
I am using the uniroot function in order to carry out a bivariate Monte
Carlo simulation using the logistics model.
I have defined the function as:
BV.FV - function(x,y,a,A)
(((x^(-a^-1)+y^(-a^-1))^(a-1))*(y^(a-1/a))*(exp(-((1^(-a^-1)+y^(-a^-1))^a)+y^-1)))-A
and the procedure is as follows:
Randomly generate values of A~(0,1), y0 = -(lnA)^-1
Where: A=Pr{Xxi|Y=yi-1} and a is the dependency between X and y (0.703)
Use y0 to determine x where x = x(i) and y = y0(i-1) in the above equation.
Unfortunately when the randomly defined A gets to approximately 0.46 the
intervals I provide no longer have the opposite signs (both -ve).
I have tried various different upper limits but I still can't seem to find a
root.
Does anyone have any suggestions?
Hi.
The expression for BV.FV() contains only one occurrence of x and it may be
separated from the equation. The solution for x may be expressed as
get.x - function(y,a,A)
((A/(y^(a-1/a))/(exp(-((1^(-a^-1)+y^(-a^-1))^a)+y^-1)))^(1/(a-1))-y^(-a^-1))^(-a)
a - 0.703
A - 0.45
y - - 1/log(A)
uniroot(BV.FV, c(1, 20), y=y, a=a, A=A)$root
[1] 3.477799
get.x(y, a, A)
[1] 3.477794
As pointed out in another reply, the equation does not always have
a solution. It seems to tend to infinity in the following
for (A in seq(0.45, 0.4677, length=20)) {
y - - 1/log(A)
cat(get.x(y, a, A), uniroot(BV.FV, c(1, 1000), y=y, a=a, A=A)$root,
\n)
}
3.477794 3.477791
3.637566 3.637565
3.812851 3.812851
4.006213 4.006214
4.220838 4.220834
4.460738 4.460739
4.731048 4.731052
5.038453 5.038473
5.391845 5.391843
5.803329 5.803329
6.289873 6.28988
6.876084 6.876084
7.5992 7.599201
8.518665 8.518683
9.736212 9.736212
11.44329 11.44328
14.05345 14.05345
18.68155 18.68155
29.97659 29.9766
219.3156 219.3155
Hope this helps.
Petr Savicky.
__
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] uniroot vs.optimize
I looked at the descriptions for uniroot and optimize and they are somewhat different but the book reference is the same and I am wondering if there are reasons to pick one over the other? Thank you. Kevin __ 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] uniroot vs.optimize
Well, there is a difference between finding a root, and finding a minimum or a maximum. So you would use one or the other depending on which you need to do. -Don At 7:00 PM -0800 11/19/09, rkevinbur...@charter.net wrote: I looked at the descriptions for uniroot and optimize and they are somewhat different but the book reference is the same and I am wondering if there are reasons to pick one over the other? Thank you. Kevin __ 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. -- - Don MacQueen Lawrence Livermore National Laboratory Livermore, CA, USA 925-423-1062 m...@llnl.gov __ 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] Uniroot and Newton-Raphson Anomaly
I have the following function for which I need to find the root of a:
f - function(R,a,c,q) sum((1 - (1-R)^a)^(1/a)) - c * q
To give context for the problem, this is a psychometric issue where R is
a vector denoting the percentage of students scoring correct on test
item i in class j, c is the proportion correct on the test by student k,
and q is the number of items on the test in total.
I have programmed this using Newton-Raphson as follows:
numer - function(R,a,c,q){
result - sum((1-(1-R)^a)^(1/a))-c*q
result
}
denom - function(R,a,c,q){
result - sum(-((1 - (1 - R)^a)^(1/a) * (log((1 - (1 - R)^a)) *
(1/a^2)) +
(1 - (1 - R)^a)^((1/a) - 1) * ((1/a) * ((1 - R)^a * log((1 -
R))
result
}
aConst - function(R, c, q, con){
a - .5 # starting value for a
change - 1
while(abs(change) con) {
r1 - numer(R,a,c,q)
r2 - denom(R,a,c,q)
change - r1/r2
a - a - change
}
a
}
The function now works as follows. Assume there are two test items on
the test. Assume the proportion correct on the items in class j is .2
and .4, and assume student k scored correctly on one item only.
aConst(R = c(.2,.4), c=.5, q=2, con=.001)
Now, one might notice that the first derivative of the function (in the
function denom) has in it log(1-R). In cases where all students in a
class answer the item correct, R = 1, and this creates an anomaly in
that log(0) is undefined. One cheap trick, I think, is to correct the
vector R such that any values of 1 become .999 and any values of 0
become .001 (0 is also an anomaly).
I am accustomed to Newton-Raphson and don't use bisection or the uniroot
function. So, given the anomaly above, I am thinking a change of mindset
might be necessary, although I am not certain if the same issue that
affects NR will propagate to other root finding algorithms.
Now, to use uniroot, my understanding is that I need to start with two
guesses for a lower and upper limit of a such that:
f(x_l)*f(x_u) 0
Where f(x_l) and f(x_u) are the lower and upper limits, respectively.
With this, I can try:
f - function(R,a,c,q) sum((1 - (1-R)^a)^(1/a)) - c * q
uniroot(f, c(.5,2), R = c(.2,.4), c = .5, q=2)
Which gives the same root as my NR function. Now, the issue I am running
into is that this function is applied to each student in the data, which
can be in the hundreds of thousands, and the only value that is fixed is
q. The values of c vary by student and the values in the vector R vary
by class. So, as I have applied this to a larger dataset, I often can't
find the root because the values I use for the search interval are not
universal to maintain the necessary condition that f(x_l)*f(x_u) 0 for
student k' in class j'. As a result, the error that the endpoints are
not of opposite sign appears.
Now, it is not the error that confuses me, that I believe is rather
transparent in meaning. It is how to apply a search interval that can be
universally applied when implementing this function to many students
when the values of R and c vary. Obviously, with hundreds of thousands
of students I cannot toy around with the function for each student to
find a search interval to maintain f(x_l)*f(x_u) 0. So, after
pondering this over the weekend, I wonder if the list might have
reactions on the following two questions;
1) First, would the issue I note about NR having log(0) propagate into
other root finding algorithms? I realize bisection does not make use of
derivatives, and this occurs in the first derivative, but I'm not savvy
enough to see if this is an artifact of the function.
2) Second, is there a more thoughtful way to identify a search interval
that can be automatically programmed when iterating over hundreds of
thousands of cases such that universal values for the search interval
can be used?
Many thanks,
Harold
sessionInfo()
R version 2.8.1 (2008-12-22)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] MiscPsycho_1.4 statmod_1.3.8
__
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] Uniroot and Newton-Raphson Anomaly
Harold,
Here is an approach that might work well for all the students (it doesn't
use derivatives):
require(BB) # I use the dfsane() root-finder from the B package
f - function(a,R,c1,q1) sum((1 - (1-R)^a)^(1/a)) - c1 * q1 # I
have re-named the variables
dfsane(par=0.01, fn=f, R = c(.2,.4), c1=.5, q1=2) # Note it works
even for extreme initial values
The dfsane() algoithm is derivative free, and is robust to poor starting
values. Granted, this is bit of an overkill to use an algorithm designed
for large-scale systems to solve a univariate problem, but it seems to work.
Hope this helps,
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvarad...@jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Doran, Harold
Sent: Monday, March 16, 2009 10:34 AM
To: r-help@r-project.org
Subject: [R] Uniroot and Newton-Raphson Anomaly
I have the following function for which I need to find the root of a:
f - function(R,a,c,q) sum((1 - (1-R)^a)^(1/a)) - c * q
To give context for the problem, this is a psychometric issue where R is a
vector denoting the percentage of students scoring correct on test item i in
class j, c is the proportion correct on the test by student k, and q is the
number of items on the test in total.
I have programmed this using Newton-Raphson as follows:
numer - function(R,a,c,q){
result - sum((1-(1-R)^a)^(1/a))-c*q
result
}
denom - function(R,a,c,q){
result - sum(-((1 - (1 - R)^a)^(1/a) * (log((1 - (1 - R)^a)) *
(1/a^2)) +
(1 - (1 - R)^a)^((1/a) - 1) * ((1/a) * ((1 - R)^a * log((1 -
R))
result
}
aConst - function(R, c, q, con){
a - .5 # starting value for a
change - 1
while(abs(change) con) {
r1 - numer(R,a,c,q)
r2 - denom(R,a,c,q)
change - r1/r2
a - a - change
}
a
}
The function now works as follows. Assume there are two test items on the
test. Assume the proportion correct on the items in class j is .2 and .4,
and assume student k scored correctly on one item only.
aConst(R = c(.2,.4), c=.5, q=2, con=.001)
Now, one might notice that the first derivative of the function (in the
function denom) has in it log(1-R). In cases where all students in a class
answer the item correct, R = 1, and this creates an anomaly in that log(0)
is undefined. One cheap trick, I think, is to correct the vector R such that
any values of 1 become .999 and any values of 0 become .001 (0 is also an
anomaly).
I am accustomed to Newton-Raphson and don't use bisection or the uniroot
function. So, given the anomaly above, I am thinking a change of mindset
might be necessary, although I am not certain if the same issue that affects
NR will propagate to other root finding algorithms.
Now, to use uniroot, my understanding is that I need to start with two
guesses for a lower and upper limit of a such that:
f(x_l)*f(x_u) 0
Where f(x_l) and f(x_u) are the lower and upper limits, respectively.
With this, I can try:
f - function(R,a,c,q) sum((1 - (1-R)^a)^(1/a)) - c * q uniroot(f, c(.5,2),
R = c(.2,.4), c = .5, q=2)
Which gives the same root as my NR function. Now, the issue I am running
into is that this function is applied to each student in the data, which can
be in the hundreds of thousands, and the only value that is fixed is q. The
values of c vary by student and the values in the vector R vary by class.
So, as I have applied this to a larger dataset, I often can't find the root
because the values I use for the search interval are not universal to
maintain the necessary condition that f(x_l)*f(x_u) 0 for student k' in
class j'. As a result, the error that the endpoints are not of opposite sign
appears.
Now, it is not the error that confuses me, that I believe is rather
transparent in meaning. It is how to apply a search interval that can be
universally applied when implementing this function to many students when
the values of R and c vary. Obviously, with hundreds of thousands of
students I cannot toy around with the function for each student to find a
search interval to maintain f(x_l)*f(x_u) 0. So, after pondering this over
the weekend, I wonder if the list might have reactions on the following two
questions;
1) First, would the issue I note about NR having log(0) propagate into other
root finding algorithms? I realize bisection does not make use of
derivatives, and this occurs in the first derivative, but I'm not savvy
enough
Re: [R] uniroot() problem
Thanks for this reply. Here I was trying to calculate implied volatility
using BS formula. This is my code :
oo = 384.40 # traded option price
uu = 1563.25 # underlying price
tt = 0.656 # time to maturity in year
ii = 2.309/100 # interest rate, annualized
th.price = function(x)
{
d1 = (log(uu/K) + (ii + x^2/2)*tt) / (x*sqrt(tt)); d2 = d1 -
x*sqrt(tt)
option.price = uu * pnorm(d1) - K * exp(-ii*tt) * pnorm(d2)
return(option.price - oo)
}
uniroot(th.price, c(-20, 20), tol=1/10^12)
I got following result :
uniroot(th.price, c(-20, 20), tol=1/10^12)
$root
[1] 6.331672e-13
$f.root
[1] 36.28816
$iter
[1] 55
$estim.prec
[1] 7.385592e-13
Hence using implied volatility, difference between traded price and
theoretical price is coming as high as 36.28816, even I increse the
precision level. Any idea how to crack this problem?
Albyn Jones wrote:
One can't tell for sure without seeing the function, but I'd guess
that you have a numerical issue. Here is an example to reflect upon:
f=function(x) (exp(x)-exp(50))*(exp(x)+exp(50))
uniroot(f,c(0,100))
$root
[1] 49.7
$f.root
[1] -1.640646e+39
$iter
[1] 4
$estim.prec
[1] 6.103516e-05
.Machine$double.eps^0.25/2
[1] 6.103516e-05
uniroot thinks it has converged, at least in relative terms. Note
that the estimated precision is related to the machine epsilon, used
in the default value for tol. try fiddling with the tol argument.
uniroot(f,c(0,100),tol=1/10^12)
$root
[1] 50
$f.root
[1] 1.337393e+31
$iter
[1] 4
$estim.prec
[1] 5.186962e-13
albyn
Quoting megh megh700...@yahoo.com:
I have a strange problem with uniroot() function. Here is the result :
uniroot(th, c(-20, 20))
$root
[1] 4.216521e-05
$f.root
[1] 16.66423
$iter
[1] 27
$estim.prec
[1] 6.103516e-05
Pls forgive for not reproducing whole code, here my question is how
f.root
can be 16.66423? As it is finding root of a function, it must be near
Zero.
Am I missing something?
--
View this message in context:
http://www.nabble.com/uniroot%28%29-problem-tp21227702p21227702.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-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/uniroot%28%29-problem-tp21227702p21909423.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] uniroot() problem
Sorry, in my previous post I forgot to include strike price, which is K =
1160
megh wrote:
Thanks for this reply. Here I was trying to calculate implied volatility
using BS formula. This is my code :
oo = 384.40 # traded option price
uu = 1563.25 # underlying price
tt = 0.656 # time to maturity in year
ii = 2.309/100 # interest rate, annualized
th.price = function(x)
{
d1 = (log(uu/K) + (ii + x^2/2)*tt) / (x*sqrt(tt)); d2 = d1 -
x*sqrt(tt)
option.price = uu * pnorm(d1) - K * exp(-ii*tt) * pnorm(d2)
return(option.price - oo)
}
uniroot(th.price, c(-20, 20), tol=1/10^12)
I got following result :
uniroot(th.price, c(-20, 20), tol=1/10^12)
$root
[1] 6.331672e-13
$f.root
[1] 36.28816
$iter
[1] 55
$estim.prec
[1] 7.385592e-13
Hence using implied volatility, difference between traded price and
theoretical price is coming as high as 36.28816, even I increse the
precision level. Any idea how to crack this problem?
Albyn Jones wrote:
One can't tell for sure without seeing the function, but I'd guess
that you have a numerical issue. Here is an example to reflect upon:
f=function(x) (exp(x)-exp(50))*(exp(x)+exp(50))
uniroot(f,c(0,100))
$root
[1] 49.7
$f.root
[1] -1.640646e+39
$iter
[1] 4
$estim.prec
[1] 6.103516e-05
.Machine$double.eps^0.25/2
[1] 6.103516e-05
uniroot thinks it has converged, at least in relative terms. Note
that the estimated precision is related to the machine epsilon, used
in the default value for tol. try fiddling with the tol argument.
uniroot(f,c(0,100),tol=1/10^12)
$root
[1] 50
$f.root
[1] 1.337393e+31
$iter
[1] 4
$estim.prec
[1] 5.186962e-13
albyn
Quoting megh megh700...@yahoo.com:
I have a strange problem with uniroot() function. Here is the result :
uniroot(th, c(-20, 20))
$root
[1] 4.216521e-05
$f.root
[1] 16.66423
$iter
[1] 27
$estim.prec
[1] 6.103516e-05
Pls forgive for not reproducing whole code, here my question is how
f.root
can be 16.66423? As it is finding root of a function, it must be near
Zero.
Am I missing something?
--
View this message in context:
http://www.nabble.com/uniroot%28%29-problem-tp21227702p21227702.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-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/uniroot%28%29-problem-tp21227702p21909458.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] uniroot() problem
Megh,
The problem is due to jump discontinuity in your function at x=0. It is
always good practice to plot the function over the range of interest.
x - seq(-20, 20, by=0.01)
plot(x, th.price(x), type=l)
This will reveal the problem. The function value jumps from -384.4 to 36.29
at x=0.
If the singularity at x=0 is not an essential one, you may be able to
anayticallty remove this singularity.
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvarad...@jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of megh
Sent: Monday, February 09, 2009 4:27 AM
To: r-help@r-project.org
Subject: Re: [R] uniroot() problem
Thanks for this reply. Here I was trying to calculate implied volatility
using BS formula. This is my code :
oo = 384.40 # traded option price
uu = 1563.25 # underlying price
tt = 0.656 # time to maturity in year
ii = 2.309/100 # interest rate, annualized
th.price = function(x)
{
d1 = (log(uu/K) + (ii + x^2/2)*tt) / (x*sqrt(tt)); d2 = d1 -
x*sqrt(tt)
option.price = uu * pnorm(d1) - K * exp(-ii*tt) * pnorm(d2)
return(option.price - oo)
}
uniroot(th.price, c(-20, 20), tol=1/10^12)
I got following result :
uniroot(th.price, c(-20, 20), tol=1/10^12)
$root
[1] 6.331672e-13
$f.root
[1] 36.28816
$iter
[1] 55
$estim.prec
[1] 7.385592e-13
Hence using implied volatility, difference between traded price and
theoretical price is coming as high as 36.28816, even I increse the
precision level. Any idea how to crack this problem?
Albyn Jones wrote:
One can't tell for sure without seeing the function, but I'd guess
that you have a numerical issue. Here is an example to reflect upon:
f=function(x) (exp(x)-exp(50))*(exp(x)+exp(50))
uniroot(f,c(0,100))
$root
[1] 49.7
$f.root
[1] -1.640646e+39
$iter
[1] 4
$estim.prec
[1] 6.103516e-05
.Machine$double.eps^0.25/2
[1] 6.103516e-05
uniroot thinks it has converged, at least in relative terms. Note
that the estimated precision is related to the machine epsilon, used
in the default value for tol. try fiddling with the tol argument.
uniroot(f,c(0,100),tol=1/10^12)
$root
[1] 50
$f.root
[1] 1.337393e+31
$iter
[1] 4
$estim.prec
[1] 5.186962e-13
albyn
Quoting megh megh700...@yahoo.com:
I have a strange problem with uniroot() function. Here is the result :
uniroot(th, c(-20, 20))
$root
[1] 4.216521e-05
$f.root
[1] 16.66423
$iter
[1] 27
$estim.prec
[1] 6.103516e-05
Pls forgive for not reproducing whole code, here my question is how
f.root
can be 16.66423? As it is finding root of a function, it must be near
Zero.
Am I missing something?
--
View this message in context:
http://www.nabble.com/uniroot%28%29-problem-tp21227702p21227702.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-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/uniroot%28%29-problem-tp21227702p21909423.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-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] uniroot() problem
I have a strange problem with uniroot() function. Here is the result : uniroot(th, c(-20, 20)) $root [1] 4.216521e-05 $f.root [1] 16.66423 $iter [1] 27 $estim.prec [1] 6.103516e-05 Pls forgive for not reproducing whole code, here my question is how f.root can be 16.66423? As it is finding root of a function, it must be near Zero. Am I missing something? -- View this message in context: http://www.nabble.com/uniroot%28%29-problem-tp21227702p21227702.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] uniroot() problem
One can't tell for sure without seeing the function, but I'd guess that you have a numerical issue. Here is an example to reflect upon: f=function(x) (exp(x)-exp(50))*(exp(x)+exp(50)) uniroot(f,c(0,100)) $root [1] 49.7 $f.root [1] -1.640646e+39 $iter [1] 4 $estim.prec [1] 6.103516e-05 .Machine$double.eps^0.25/2 [1] 6.103516e-05 uniroot thinks it has converged, at least in relative terms. Note that the estimated precision is related to the machine epsilon, used in the default value for tol. try fiddling with the tol argument. uniroot(f,c(0,100),tol=1/10^12) $root [1] 50 $f.root [1] 1.337393e+31 $iter [1] 4 $estim.prec [1] 5.186962e-13 albyn Quoting megh megh700...@yahoo.com: I have a strange problem with uniroot() function. Here is the result : uniroot(th, c(-20, 20)) $root [1] 4.216521e-05 $f.root [1] 16.66423 $iter [1] 27 $estim.prec [1] 6.103516e-05 Pls forgive for not reproducing whole code, here my question is how f.root can be 16.66423? As it is finding root of a function, it must be near Zero. Am I missing something? -- View this message in context: http://www.nabble.com/uniroot%28%29-problem-tp21227702p21227702.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-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.
