[R] using "integrate" in a function definition
theo borm
theo_rstt at borm.org
Fri Feb 23 19:42:50 CET 2007
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.
More information about the R-help
mailing list