6

I would like to run some tests for the null hypothesis that the times of events I have was created from a homogeneous Poisson process (see e.g. http://en.wikipedia.org/wiki/Poisson_process ). For a fixed number of events the times should therefore look like a sorted version of a uniform distribution in the appropriate range. There is an implementation of the Kolmogorov-Smirnov test at http://docs.scipy.org/doc/scipy-0.7.x/reference/generated/scipy.stats.kstest.html but I can't see how to use it here as scipy.stats doesn't seem to know about Poisson processes.

As a simple example, this sample data should give a high p-value for any such test.

import random
nopoints = 100
max = 1000

points = sorted([random.randint(0,max) for j in xrange(nopoints)])

How can I make a sensible test for this problem?

From www.stat.wmich.edu/wang/667/classnotes/pp/pp.pdf‎ I see

" REMARK 6.3 ( TESTING POISSON ) The above theorem may also be used to test the hypothesis that a given counting process is a Poisson process. This may be done by observing the process for a fixed time t. If in this time period we observed n occurrences and if the process is Poisson, then the unordered occurrence times would be independently and uniformly distributed on (0, t]. Hence, we may test if the process is Poisson by testing the hypothesis that the n occurrence times come from a uniform (0, t] population. This may be done by standard statistical procedures such as the Kolmogorov-Smirov test."

3
  • scipy.stats does know a poisson distribution: >>> import scipy.stats as stats >>> stats.poisson.pmf(1, mu=3) 0.14936120510359185
    – ev-br
    Sep 17, 2013 at 0:14
  • The problem with using the Kolmogorov-Smirnov test is that the distribution of the test statistic is not the same for discrete and continuous distributions, and the same if you estimate parameters. The p-values from kstest assume that it's a continuous distribution with no estimated parameters, and the p-values will not be correct if that's not the case in your problem. (applies also to Hooked's answer)
    – Josef
    Sep 17, 2013 at 11:35
  • If you have the continuous times for Remark 6.3, then you can still use kstest with the uniform as Null distribution. If you have discrete events or event times, then you can use the chisquare distribution. There was just the same question on the statsmodels mailing list with other tests for homogeneity of a Poisson process (for example with binning into time intervals), but I don't have any ready to use tests for this.
    – Josef
    Sep 17, 2013 at 17:08

3 Answers 3

4

Warning: quickly written, some details are not verified

what's the appropriate estimator for exponential, degrees of freedom for chisquare test

based on lecture notes

The implications of homogeneity is not rejected with any of the three tests. Illustrates how to use kstest and chisquare test of scipy.stats

# -*- coding: utf-8 -*-
"""Tests for homogeneity of Poissson Process

Created on Tue Sep 17 13:50:25 2013

Author: Josef Perktold
"""

import numpy as np
from scipy import stats

# create an example dataset
nobs = 100
times_ia = stats.expon.rvs(size=nobs) # inter-arrival times
times_a = np.cumsum(times_ia) # arrival times
t_total = times_a.max()

# not used
#times_as = np.sorted(times_a)
#times_ia = np.diff(times_as)

bin_limits = np.array([ 0. ,  0.5,  1. ,  1.5,  2. ,  np.inf])
nfreq_ia, bins_ia = np.histogram(times_ia, bin_limits)


# implication: arrival times are uniform for fixed interval
# using times.max() means we don't really have a fixed interval
print stats.kstest(times_a, stats.uniform(0, t_total).cdf)

# implication: inter-arrival times are exponential
lambd = nobs * 1. / t_total
scale = 1. / lambd

expected_ia = np.diff(stats.expon.cdf(bin_limits, scale=scale)) * nobs
print stats.chisquare(nfreq_ia, expected_ia, ddof=1)

# implication: given total number of events, distribution of times is uniform
# binned version
n_mean_bin = 10
n_bins_a = nobs // 10
bin_limits_a = np.linspace(0, t_total+1e-7, n_bins_a + 1)
nfreq_a, bin_limits_a = np.histogram(times_a, bin_limits_a)
# expect uniform distributed over every subinterval
expected_a = np.ones(n_bins_a) / n_bins_a * nobs
print stats.chisquare(nfreq_a, expected_a, ddof=1)
2

The KS test, when determining if two distributions differ, is simply the largest difference between them:

enter image description here

This is simple enough to calculate yourself. The program below computes the KS statistic for two Poisson processes with different parameter sets:

import numpy as np

N = 10**6
X  = np.random.poisson(10, size=N)
X2 = np.random.poisson(7, size=N)

bins = np.arange(0, 30,1)
H1,_ = np.histogram(X , bins=bins, normed=True)
H2,_ = np.histogram(X2, bins=bins, normed=True)

D = np.abs(H1-H2)

idx = np.argmax(D)
KS = D[idx]

# Plot the results
import pylab as plt
plt.plot(H1, lw=2,label="$F_1$")
plt.plot(H2, lw=2,label="$F_2$")
text = r"KS statistic, $\sup_x |F_1(x) - F_2(x)| = {KS:.4f}$"
plt.plot(D, '--k', label=text.format(KS=KS),alpha=.8)
plt.scatter([bins[idx],],[D[idx],],s=200,lw=0,alpha=.8,color='k')
plt.axis('tight')
plt.legend()

enter image description here

3
  • @Anush The KS test assumes you have a distribution you are comparing it against. It sounds like you have a model in mind - the Poisson distribution. In this case you would change F1 to your data above. When looking at a distribution, I'm not sure what you mean by the ordering, normally I think of the CDF as resulting from an average. The ordering should be irrelevant, in fact a Poisson process assumes idd, something stronger.
    – Hooked
    Sep 16, 2013 at 20:07
  • @Anush to go from a series of observations (your exponentially distributed times) to a distribution you can simply histogram your observations into bins (as shown in the example) and normalize.
    – Hooked
    Sep 16, 2013 at 20:12
  • Thanks. I think my confusion is just how to choose the bin sizes or whether I really have to bin at all.
    – Simd
    Sep 17, 2013 at 8:14
0

The problem is, as the doc you linked to suggests: "The KS test is only valid for continuous distributions." while a poisson distribution is discrete.

I would suggest to you maybe use the example in this link: http://nbviewer.ipython.org/urls/raw.github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/master/Chapter1_Introduction/Chapter1_Introduction.ipynb (look for "##### Example: Inferring behaviour from text-message data")

In that link, they check for the appropriate lambda(s) for a particular dataset which they assume distributes according to a poisson process.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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