# Peacock Test Implementation in R

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?

• 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?
– whuber
Jun 24 '19 at 19:42
• I'm using the formula on page 622 (8 of the paper linked in the post.) Jun 25 '19 at 12:20
• 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.
– whuber
Jun 25 '19 at 13:32
• @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
– mjs
Nov 28 '19 at 14:41
• @user1329307 see also discussion on the C -> R translation topic stackoverflow.com/questions/58816535/…
– mjs
Nov 28 '19 at 14:48

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)
(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))