strptime が面白い

strptime が面白い

> strptime("2014年12月27日 14時27分34秒", "%Y年%m月%d日 %H時%M分%S秒")
[1] "2014-12-27 14:27:34 JST"

> strptime("2014/12/27 14-27-34", "%Y/%m/%d %H-%M-%S")
[1] "2014-12-27 14:27:34 JST"

> strptime("2014.12.27 14.27.34", "%Y.%m.%d %H.%M.%S")
[1] "2014-12-27 14:27:34 JST"

> strptime("14-27-34 2014/12/27", "%H-%M-%S %Y/%m/%d")
[1] "2014-12-27 14:27:34 JST"

> strptime("14-27-34 2014/12/27", "%H-%M-%S %Y/%m/%d")
[1] "2014-12-27 14:27:34 JST"

> strptime("14-27 2014/12/27", "%H-%M %Y/%m/%d")
[1] "2014-12-27 14:27:00 JST"

> strptime("2014/12/27", "%Y/%m/%d ")
[1] "2014-12-27 JST"

> strptime("14-27-34", "%H-%M-%S")
[1] "2015-01-08 14:27:34 JST"

適合度の検定について（2）

「適合度の検定について」に対し，中澤さんからコメントをいただいた。

> vcd パッケージの goodfit が便利です

はじめて使ってみました。

> goodfit(as.table(tab1))

Observed and fitted values for poisson distribution
with parameters estimated by `ML'

 count observed       fitted
     0        0 2.150092e-03
     1        0 2.805870e-02
     2        0 1.830830e-01
     3        1 7.964111e-01
     4        1 2.598291e+00
     5        3 6.781540e+00
     6       12 1.474985e+01
     7       27 2.749793e+01
     8       42 4.485600e+01
     9       63 6.504121e+01
    10       80 8.487877e+01
    11      108 1.006971e+02
    12      129 1.095081e+02
    13      108 1.099293e+02
    14      106 1.024698e+02
    15       87 8.914871e+01
    16       79 7.271192e+01
    17       45 5.581709e+01
    18       42 4.046739e+01
    19       31 2.779471e+01
    20       18 1.813605e+01
    21        8 1.127026e+01
    22        5 6.685312e+00
    23        1 3.793188e+00
    24        1 2.062546e+00
    25        2 1.076649e+00
    26        1 5.403950e-01

> summary(goodfit(as.table(tab1)))

     Goodness-of-fit test for poisson distribution

                      X^2 df  P(> X^2)
Likelihood Ratio 19.73707 22 0.5994826

fitted の合計が 999.5218 となるのがちょっといやな感じです（x >= 27 に対する確率が考慮されていない）。

goodfit がやっていることは，以下の通り（冗長な部分を削除）。

x = 3:26
tab1 = c(1, 1, 3, 12, 27, 42, 63, 80, 108, 129, 108, 106, 87, 79, 45, 42, 31, 18, 8, 5, 1, 1, 2, 1)
names(tab1) = x
lambda = sum(x * tab1) / sum(tab1)
p = dpois(x, lambda)
cbind(x, tab1, p, p*sum(tab1))
ln = function(o, e) sum(ifelse(o == 0, 0, o*log(o/e)))
( G2 = sum(ln(tab1, p*sum(tab1)))*2 )
( df = length(x) - 1 -1 )
( p.value = pchisq(G2, df, lower.tail=FALSE) )

結果

> ( G2 = sum(ln(tab1, p*sum(tab1)))*2 )
[1] 19.73707
> ( df = length(x) - 1 -1 )
[1] 22
> ( p.value = pchisq(G2, df, lower.tail=FALSE)　)
[1] 0.5994826

G2（G スクエアー）を使う検定（独立性の検定でも，G2 を使う検定が定義されている）。

この検定では，観察値が 0 のセルは検定統計量に反映されないので，観察値が 0 のセルがあっても良いし，「期待値が小さいセルについての考慮」も不要ということにはなっている。

適合度の検定について

「交代再生過程のある時刻における状態がポアソン分布に従うか否かシミュレーションしてみた.(非公開にされました)」において，

> 検定の仕方まちがえてるのだろうか？

ということですが。そのとおり，残念ながら，間違えています。

著者が最初にやった chisq.test(tab1) は，一様分布に従うかどうかの検定です。一様とはとてもいえないので，帰無仮説は棄却されてしまいますね。

やるべき検定は，適合度検定です。chisq.test 関数を使います。

一様分布に従うかどうかも，適合度検定ですが，chisq.test の引数に p というのがあり，分布が理論比 p に従うかどうかの適合度の検定を行うのです。n 個のカテゴリーがある場合の一様分布の適合度検定は p = 1/n がデフォルトになっているのです。

ですから，ポアソン分布などへの適合度検定の場合には p を明示的に指定してやらなければなりません。これが，一筋縄ではいかないことがあります。

まず，最初にポアソン分布に従うという帰無仮説のもとで適合度検定を行うわけですが，ポアソン定数が分かっているならそれを使えばよいのですが，もし標本から推定したポアソン定数（平均値）を使うならば，後述するカイ二乗分布の自由度が余分に 1 減ります。

確率変数（例の場合は 3 ～ 26）をとる確率は dpois で求まります（x1 = 3:26); dpois(x1, l1)）が，ポアソン分布の確率変数は 0 から ∞ ですが，0 近辺や 大きな値をとる確率は小さくなるので，小さい方は 0 から ある値までの確率を合計します（プールします）。同様にある値以上の確率も合計します。「ある値」は，期待値が 1 以上になるような値です。幸いなことに例題の場合は，最小と最大（3 と 26）になります。期待値が 1 以上になるようにするためにはもっと大きい，あるいは，小さい値までを含めるなければならないこともあるでしょう。この場合には，度数分布表もまとめる必要があるでしょう（そして，そのような場合には結果として自由度が余分に減少します）。
ちょっとややこしいですが，丁寧に説明してあるページなどを参照しましょう。

で，最終的には以下に示すデータフレームのような数値を使って検定を行います。

         x1 tab1       Dpois Expectation
0 thru 3  3    1 0.001009703    1.009703
4         4    1 0.002598291    2.598291
5         5    3 0.006781540    6.781540
6         6   12 0.014749849   14.749849
7         7   27 0.027497933   27.497933
8         8   42 0.044856004   44.856004
9         9   63 0.065041206   65.041206
10       10   80 0.084878773   84.878773
11       11  108 0.100697090  100.697090
12       12  129 0.109508086  109.508086
13       13  108 0.109929271  109.929271
14       14  106 0.102469784  102.469784
15       15   87 0.089148712   89.148712
16       16   79 0.072711919   72.711919
17       17   45 0.055817090   55.817090
18       18   42 0.040467391   40.467391
19       19   31 0.027794708   27.794708
20       20   18 0.018136047   18.136047
21       21    8 0.011270258   11.270258
22       22    5 0.006685312    6.685312
23       23    1 0.003793188    3.793188
24       24    1 0.002062546    2.062546
25       25    2 0.001076649    1.076649
26 over  26    1 0.001018650    1.018650

自由度は，カテゴリー数（階級数）マイナス 2 ですが，マイナス 1 は度数の合計がサンプルサイズになるという制約，もう一つのマイナス 1 は，標本からポアソン定数を決めたための制約です。

ということで，結果は，カイ二乗値=16.18529，自由度=22，P値=0.8065844ということになり，「ポアソン分布に従っていないとはいえない」という結論になります。

なお，二項分布への適合度検定も同じように行えますが，その場合も，帰無仮説は棄却できないようです。

x1 = 3:26
tab1 = c(1, 1, 3, 12, 27, 42, 63, 80, 108, 129, 108, 106, 87, 79, 45, 42, 31, 18, 8, 5, 1, 1, 2, 1)
n = sum(tab1)

# chisq.test(d$freq) # これは一様分布の検定なので間違い

( l1 = sum(x1 * tab1) / n ) # lambda
Dpois = dpois(x1, l1)
Dpois[1] = sum(dpois(0:3, l1)) # = ppois(3, l1)
Dpois[24] = 1-sum(Dpois[1:23]) # = ppois(25, l1, lower.tail=FALSE)
Expectation = Dpois * n
d = data.frame(x1, tab1, Dpois, Expectation)
rownames(d) = c("0 thru 3", 4:25, "26 over")
d
( chisq = sum((tab1-Expectation)^2/Expectation) )
( df = length(tab1) - 1 - 1 )
( p.value = pchisq(chisq, df, lower.tail=FALSE) )

結果

> ( chisq = sum((tab1-Expectation)^2/Expectation) )
[1] 16.18529
> ( df = length(tab1) - 1 - 1 )
[1] 22
> ( p.value = pchisq(chisq, df, lower.tail=FALSE) )
[1] 0.8065844