ポアソン分布に従うカウントデータの平均値の差の検定 ref: http://qiita.com/kndt84/items/1c855df473d6a256a191

file0.txt

\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}

file1.txt

pois_mean_diff_test(40, 65)
=> 0.04921332025427114

file2.txt

\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

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

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