Note that warning message. Its trying to evaluate qnorm outside of its allowable domain.
On 3/12/06, Tolga Uzuner <[EMAIL PROTECTED]> wrote: > 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 > >>> > >>>______________________________________________ > >>>R-help@stat.math.ethz.ch 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 > >> > >>______________________________________________ > >>R-help@stat.math.ethz.ch mailing list > >>https://stat.ethz.ch/mailman/listinfo/r-help > >>PLEASE do read the posting guide! > >>http://www.R-project.org/posting-guide.html > >> > >> > >> > > > > > > > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html