Your function "jjj" is not vectorized. Try this:
jjj <- function(www) sapply(www, function(x)2*integrate(dnorm,0,x)$value) plot(jjj, 0, 5) It should work. 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: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of theo borm Sent: Friday, February 23, 2007 1:43 PM To: r-help@stat.math.ethz.ch Subject: [R] using "integrate" in a function definition Dear list members, I'm quite new to R, and though I tried to find the answer to my probably very basic question through the available resources (website, mailing list archives, docs, google), I've not found it. If I try to use the "integrate" function from within my own functions, my functions seem to misbehave in some contexts. The following example is a bit silly, but illustrates what's happening: The following behaves as expected: > jjj<-function(www) {2*integrate(dnorm,0,www)$value} > kkk<-function(www) {2*(pnorm(www)-0.5)} > jjj(2) [1] 0.9544997 > kkk(2) [1] 0.9544997 However, if I do: > plot(jjj,0,5) Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ whereas > plot(kkk,0,5) produces a nice plot. > xxx<-seq(0:5) > yyy<-jjj(xxx) > zzz<-kkk(xxx) produces no errors, but: > yyy [1] 0.6826895 > zzz [1] 0.6826895 0.9544997 0.9973002 0.9999367 0.9999994 1.0000000 Why is this? Is this some R language feature that I've completely missed? Ultimately my problem is that I have a mathematical function describing a distribution, and I want to use this in a Kolmogorov-Smirnov test where I need a cumulative distribution function. If I use the following (synthetic) dataset with ks.test with either the "jjj" or "kkk" function, I get: > ppp<-c(1.74865955,1.12220426,0.24760427,0.24351439,0.10870853, 0.72023272,0.40245392,0.16002948,0.24203950,0.44029851, 0.34830446,1.66967152,1.71609574,1.17540830,0.22306379, 0.64382628,0.37382795,0.84361547,1.92138362,0.02844235, 0.11238473,0.21237557,1.01058073,2.33108243,1.36034473, 1.88951679,0.18230647,0.96571916,0.91239760,2.05574766, 0.33681130,2.76006257,0.83952491,1.24038831,1.46199671, 0.24694163,0.01565159,0.94816108,1.04642673,0.36556444, 2.20760578,1.59518590,0.83951030,0.07113652,0.97422886, 0.59835793,0.08383890,1.09180853,0.43508943,0.35368637, 0.75805978,0.12790161,1.56798563,1.68669770,0.56168021) > ks.test(ppp,kkk) One-sample Kolmogorov-Smirnov test data: ppp D = 0.1013, p-value = 0.5895 alternative hypothesis: two.sided [ which seems correct as the samples come from abs(rnorm()) ] > ks.test(ppp,jjj) One-sample Kolmogorov-Smirnov test data: ppp D = 0.9875, p-value < 2.2e-16 alternative hypothesis: two.sided [ which is *incorrect* ] Note: This example is artificial as I have used a function for which I know the integral; my real problem involves a distribution that I can only integrate numerically. How would I go about to solve this problem? With kind regards, Theo. ______________________________________________ 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 and provide commented, minimal, self-contained, reproducible code. ______________________________________________ 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 and provide commented, minimal, self-contained, reproducible code.