ポアソン分布に従うカウントデータの平均値の差の検定

はじめに

業務でポアソン分布に従うカウントデータの平均値の差の検定を行う必要があったのですが、日本語の情報がなく、下記の論文を参考に実装を行ったので備忘録として残しておきます。

[A more powerful test for comparing two Poisson means]
(http://www.sciencedirect.com/science/article/pii/S0378375802004081)

検定の方法には、the conditional test (C-test) と呼ばれる方法と、この論文で提案されているP値を使った検定 (E-test) があるとのことでしたが、今回は検出力の高い E-test の方を実装しました。

簡単な説明

$n1, n2$：単位時間の経過回数
$k1, k2$：事象の発生回数の全期間の合計
$d$：検定したい平均の差

としたとき、

帰無仮説 $H_0: \lambda_1 = \lambda_2 + d$
対立仮説 $H_1: \lambda_1 \neq \lambda_2 + d$

の仮説検定を行う場合のP値は下記の式で求められる。

\begin{align}
& \hat{\lambda}_{2k} = \frac{k_1 + k_2}{n_1 + n_2} - \frac{d n_1}{n_1 + n_2}\\[2mm]

& p = \sum_{x_1=0}^{k_1} \sum_{x_2=0}^{k_2}
\frac{e^{-n_1(\hat{\lambda}_{2k}+d)}(n_1(\hat{\lambda}_{2k} + d))^{x_1}}{x_1!}
\frac{e^{-n_2\hat{\lambda}_{2k}}(n_2(\hat{\lambda}_{2k}))^{x_2}}{x_2!}
\end{align}

実装

使い方

上記の式をPython と Ruby で実装したものが下記になります。例えば、比較したいカウント数がそれぞれ、40と65だった場合、２つの平均値に有意な差があるかどうか知りたい場合は、

pois_mean_diff_test(40, 65)
=> 0.04921332025427114

とするとP値を計算することができます。今回は、P値が 0.05 を下回ることがわかるので、5%水準で有意に差があるということがわかります。

注意点

数式をそのまま実装すると、浮動小数点の計算でオーバーフローするので、階乗の計算を下記のように分解して実装しました。

\frac{e^{-n_1(\hat{\lambda}_{2k}+d)}(n_1(\hat{\lambda}_{2k} + d))^{x_1}}{x_1!}
= e^{-n_1(\hat{\lambda}_{2k}+d)} \times
\frac{n_1(\hat{\lambda}_{2k} + d)}{x_1} \times \frac{n_1(\hat{\lambda}_{2k} + d)}{x_1 -1} \times \dots \times \frac{n_1(\hat{\lambda}_{2k} + d)}{1}

Pythonでの実装

Python

import math
import numpy as np

def pois_mean_diff_test(k1, k2, n1=1, n2=1, d=0.0):

    x1_seq = range(0, k1 + 1)
    x2_seq = range(0, k2 + 1)
    l2k = (k1+k2)/(n1+n2) - d*n1/(n1+n2)

    p_value = sum([math.exp(-n1*(l2k+d)) * np.prod([n1*(l2k+d)/i for i in range(1, x1+1)]) * \
                   math.exp(-n2*l2k) * np.prod([n2*l2k/j for j in range(1, x2+1)]) \
                   for x2 in x2_seq for x1 in x1_seq])

    return p_value

if __name__ == '__main__':
    print "P-value is " + str(pois_mean_diff_test(40, 65))

Rubyでの実装

Ruby

def pois_mean_diff_test(k1, k2, n1=1, n2=1, d=0.0, alpha=0.05)

  x1_seq = Array(0..k1)
  x2_seq = Array(0..k2)

  l2k = (k1+k2)/(n1+n2) - d*n1/(n1+n2)

  p_value = x1_seq.product(x2_seq).map{|x| x1=x[0]; x2=x[1];
    Math.exp(-n1*(l2k+d)) * Array(1..x1).map{|i| n1*(l2k+d)/i }.inject(:*).to_f *
    Math.exp(-n1*l2k) * Array(1..x2).map{|j| n2*l2k/j }.inject(:*).to_f
  }.inject(:+)

  return p_value
end

out = "P-value is " + pois_mean_diff_test(20,10).to_s
puts out

参考文献

[K. Krishnamoorthy, Jessica Thomson, A more powerful test for comparing two Poisson means, Journal of Statistical Planning and Inference, Volume 119, Issue 1, 15 January 2004, Pages 23-35]
(http://www.sciencedirect.com/science/article/pii/S0378375802004081)