Actually, I did implement this using richardson extrapolation, but am
having trouble vectorising it. For some reason, it fails within integrate...
Anyone willing to look over the below and let me know what I am doing
wrong, helps much appreciated. You can cut paste the below into the
console..
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
richardson.grad <- function(func, x10, d=0.01, eps=1e-4, r=6, show=F){
sapply(x10,function(x){
v <- 2 # reduction factor.
n <- length(x) # Integer, number of variables.
a.mtr <- matrix(1, r, n)
b.mtr <- matrix(1, (r - 1), n)
h <- abs(d*x)+eps*(x==0.0)
for(k in 1:r) { # successively reduce h
for(i in 1:n) {
x1.vct <- x2.vct <- x
x1.vct[i] <- x[i] + h[i]
x2.vct[i] <- x[i] - h[i]
if(k == 1) a.mtr[k,i] <- (func(x1.vct) - func(x2.vct))/(2*h[i])
else{
if(abs(a.mtr[(k-1),i])>1e-20)
# some functions are unstable near 0.0
a.mtr[k,i] <- (func(x1.vct)-func(x2.vct))/(2*h[i])
else a.mtr[k, i] <- 0
}
}
h <- h/v # Reduced h by 1/v.
}
if(show) {
cat("\n","first order approximations", "\n")
print(a.mtr, 12)
}
for(m in 1:(r - 1)) {
for(i in 1:(r - m)) b.mtr[i,]<- (a.mtr[(i+1),]*(4^m)-a.mtr[i,])/(4^m-1)
if(show & m!=(r-1) ) {
cat("\n","Richarson improvement group No. ", m, "\n")
print(a.mtr[1:(r-m),], 12)
}
}
a.mtr[length(a.mtr)]})
}
## try it out
richardson.grad(function(x){x^3},2)
#works fine... should return 12.
# now try integrating something simple
integrate(function(i) richardson.grad(function(x) x^2,i),0,1)
#also works fine, but instead try this:
CDFLHP <-function(x,D,B)
pnorm((sqrt(1-B^2)*qnorm(x)-D)/B)
integrate(function(i) richardson.grad(function(x) CDFLHP(x,-2,0.1),i),0,1)
# fails, for some annoying reason, even tho richardson.grad is vectorised...
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
This is the error message:
Error in if (abs(a.mtr[(k - 1), i]) > 1e-20) a.mtr[k, i] <-
(func(x1.vct) - :
missing value where TRUE/FALSE needed
In addition: Warning message:
NaNs produced in: qnorm(p, mean, sd, lower.tail, log.p)
Peter Dalgaard wrote:
>"Gray Calhoun" <[EMAIL PROTECTED]> writes:
>
>
>
>>Tolga,
>>
>>Look at numericDeriv.
>>
>>
>>
>>>arbfun <- function(x) x^2
>>>x <- 3
>>>numericDeriv(quote(arbfun(x)), "x")
>>>
>>>
>>[1] 9
>>attr(,"gradient")
>> [,1]
>>[1,] 6
>>
>>
>
>However, numericDeriv is not particularly intelligent. It is
>effectively doing what Tolga was trying not to. A more refined
>function could be a good idea, e.g. implementing higher order
>approximations, a tunable stepsize, box constraints...
>
>
>
>>--Gray
>>
>>On 3/12/06, Tolga Uzuner <[EMAIL PROTECTED]> wrote:
>>
>>
>>>Hi,
>>>
>>>Suppose I have an arbitrary function:
>>>
>>>arbfun<-function(x) {...}
>>>
>>>Is there a robust implementation of a numerical derivative routine in R
>>>which I can use to take it's derivative ? Something a bit more than
>>>simple division by delta of the difference of evaluating the function at
>>>x and x+delta...
>>>
>>>Perhaps there is a way to do this using D or deriv but I could not
>>>figure it out. Trying:
>>>
>>>eval(deriv(function(x) arbfun(x),"x"),1)
>>>
>>>does not seem to work.
>>>
>>>Thanks,
>>>Tolga
>>>
>>>______________________________________________
>>>[email protected] mailing list
>>>https://stat.ethz.ch/mailman/listinfo/r-help
>>>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>>>
>>>
>>>
>>--
>>Gray Calhoun
>>
>>Economics Department
>>UC San Diego
>>
>>______________________________________________
>>[email protected] mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>>
>>
>>
>
>
>
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
