4
$\begingroup$

I have two empirical distribution (without access to the observations themselves). So, in order to see if the empirical distributions are "significantly" different, I took two large random samples from the empirical distributions and tried a few different tests in R -- Cramer test, non-parametric density equality test, and Peacock's test.

#library(cramer)
#cramer::cramer.test(x=dat1, y=dat2, sim="permutation")
#library(np)
#npdeneqtest(x=data.frame(dat1), y=data.frame(dat2))
#library(Peacock.test)
#D <- peacock2(dat1, dat2)
#Z <- sqrt(n1*n2/(n1+n2))*D
#2*exp(-2*(Z-0.5)^2)

I used this code and the result of peacock2 did not match the other two tests. This led me to start playing around with the peacock test to see what was going on.

> mat.unif <- matrix(data=rep(1/100, 100), ncol = 10, nrow=10)
> n1 <- 1000
> n2 <- n1
> data1 <- arrayInd(sample(size=n1, x=length(mat.unif), 
    prob=c(mat.unif), replace=T), dim(mat.unif))
> (D <- peacock2(data1, data1))
[1] 0.119
> (Z<-sqrt(n1*n2/(n1+n2))*(D))
[1] 2.549117
> (p.val<-2*exp(-2*(Z-0.5)^2))
[1] 0.0004507409

This doesn't make any sense to me because both datasets are the same -- the p-value should be very close to 1. The function returns the test statistic and I think I correctly translated this paper to r code.

Is the R-package broken, or am I calculating the p-value incorrectly?

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Could you explain your calculation of the p-value? The formula $2\exp(-2(Z-0.5)^2)$ produces results between $0$ and $2$ and therefore does not correspond to any possible valid calculation, so it's wrong; but what calculation are you trying to carry out? $\endgroup$
    – whuber ♦
    Jun 24, 2019 at 19:42
  • $\begingroup$ I'm using the formula on page 622 (8 of the paper linked in the post.) $\endgroup$
    – user1329307
    Jun 25, 2019 at 12:20
  • $\begingroup$ Thank you for pointing to that. It's an asymptotic formula for small p-values, which explains why it might work. However, in comparing the paper to your code I find that your code does not execute what you claim: in particular, the line (D <- peacock2(data1, data1)) compares the dataset data1 to itself. Note, too--unless I missed some special provision--the paper's results do not apply to your data, anyway: it assumes the data are governed by a continuous probability law, ruling out the chance of any ties in the datasets. $\endgroup$
    – whuber ♦
    Jun 25, 2019 at 13:32
  • $\begingroup$ @user1329307, have you tried 'Kolmogorov-Smirnov Test for Two-Dimensional Data' by William H. Press, and Saul A. Teukolsky? aip.scitation.org/doi/abs/10.1063/1.4822753 The C code is relatively easy to translate to R or can be used direclty. It is a simplified version of the Peacock test as simplified initially by Fasano and G. & Franceschini, A. See adsabs.harvard.edu/full/1987MNRAS.225..155F $\endgroup$
    – mjs
    Nov 28, 2019 at 14:41
  • $\begingroup$ @user1329307 see also discussion on the C -> R translation topic stackoverflow.com/questions/58816535/… $\endgroup$
    – mjs
    Nov 28, 2019 at 14:48

1 Answer 1

Reset to default
0
$\begingroup$

if you jitter the data, it works although p>1 when I ran it - the paper acknowledges this possibility

data11<-data1+0.01*rnorm(2000)
data12<-data1+0.01*rnorm(2000)
head(data11)
head(data12)
(D<-peacock2(data11,data12))

there's also a small adjustment to use Zinf - see the summary on p12 of the pdf - but almost no difference with your example

n=500
(Zn<-n^0.5*D)
(Zinf<-Zn*(1-0.53*n^-0.9)^-1)
(p<-2*exp(-2*(Zinf-0.5)^2))

data11 and data12 are extremely close, leading to p>1. If you resample to generate data2, it is well behaved

data2<-arrayInd(sample(size=n1,x=length(mat.unif),prob=c(mat.unif),
replace=T),dim(mat.unif))
data21<-data2+0.01*rnorm(2000)
(D<-peacock2(data11,data21))
(Zn<-n^0.5*D)
(Zinf<-Zn*(1-0.53*n^-0.9)^-1)
(p<-2*exp(-2*(Zinf-0.5)^2))
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.