# Calculating p-value from an arbitrary distribution

I hope this isn't a silly question. Let's say I have some arbitrary continuous distribution. I also have a statistic, and I'd like to use this arbitrary distribution to get a p-value for this statistic.

I realize that in R it's easy to do this as long as your distribution fits one of the built-in ones, like if it's normal. But is there an easy way to do this with any given distribution, without making that kind of assumption?

-
It's a silly question if I can answer it. –  bill_080 Feb 6 '11 at 22:04
@bill_080 As you can see, it is more complex in practice. –  mbq Feb 7 '11 at 10:49
@mbq In theory, even I could be correct –  bill_080 Feb 7 '11 at 18:21

If you have a cumulative distribution function $F$, then calculating the $p$-value for given statistic $T$ is simply $1-F(T)$. This is straightforward in R. If you have probability density function on the other hand, then $F(x)=\int_{-\infty}^xp(t)dt$. You can find this integral analytically or numerically. In R this will look like this:

dF <- function(x)dnorm(x)
pF <- function(q)integrate(dF,-Inf,q)$value > pF(1) [1] 0.8413448 > pnorm(1) [1] 0.8413447  You can tune integrate for better accuracy. This of course may fail for specific cases, when the integral does not behave well, but it should work for majority of the density functions. You can of course pass parameters into pF, if you have several parameter values to try-out and do not want to redefine dF each time. dF <- function(x,mean=0,sd=1)dnorm(x,mean=mean,sd=sd) pF <- function(q,mean=0,sd=1)integrate(dF,-Inf,q,mean=mean,sd=sd)$value

> pF(1,1,1)
[1] 0.5
> pnorm(1,1,1)
[1] 0.5


Of course you can also use Monte-Carlo methods as detailed by @suncoolsu, this would be just another numerical method for integration.

-
I think your method is way easier than what I suggested, especially if there are no restriction on the function that you are integrating. I am not aware of the numerical tech. in R. –  suncoolsu Feb 7 '11 at 22:12
Yes, I think this is more in line with my current abilities. Thanks! –  Alan H. Feb 7 '11 at 23:19
Actually, I'm not sure I quite follow how these functions work. The examples give results for a normal distribution, but where do I plug in my probability density function? –  Alan H. Feb 12 '11 at 23:19
(I've done tests, and my data doesn't seem to be remotely normal.) –  Alan H. Feb 12 '11 at 23:40
@Alan H., plug your density function in to dF. That is dF should return density function value at given argument. –  mpiktas Feb 13 '11 at 4:39

Yes, it is possible to use any arbitrary distribution to get a p-value for any statistic. Theoretically and practically you can calculate (one-sided) p-value by this formula.

$$\mathrm{p-value} = P[T > T_{observed} | H_0 \quad \mathrm{holds}]$$

Where $T$ is the test-statistic of interest and $T_{observed}$ is the value that you have calculated for the observed data.

If you know the theoretical distribution of $T$ under $H_0$, great! Otherwise, you can use MCMC simulation to generate from the null distribution of $T$ and calculate the Monte Carlo integral to obtain p-value. Numerical integration techniques will also work in case you don't want to use (may be) easier Monte Carlo methods (especially in R; in Mathematica integration may be easier, but I have no experience using it)

The only assumption you are making here is -- you know the null distribution of T (which may not be in the standard R random number generator formats). That's it -- as long as you know the null distribution, p-value can be calculated.

-
I must remark -- this is one the reasons why p-values are so popular and easy to be misunderstood. (IMHO) –  suncoolsu Feb 6 '11 at 18:47
Ok, this makes sense. I do have what I believe is a good estimation of the null distribution. Any hints on how to implement this in R? Thanks! –  Alan H. Feb 6 '11 at 19:07
@Alan - Do you know how to generate random values from your Null Distribution? If yes, suppose - T= c(T1, ..., TN) are draws from the null distribution - p-value = sum(T > T_obs)/N. If you don't know how to generate, may you need to use Metropolis Sampling or Gibbs Sampling to get T1...TN, but it is very do-able. –  suncoolsu Feb 6 '11 at 19:30
Ah! I get it. Thanks again. –  Alan H. Feb 6 '11 at 19:37