ぎるばーとの日記

もっともっと遠くへ行きたい 空が広く見える場所まで

2×2分割表の検定

 久しぶりにこちらのブログに。何ヶ月ぶりだろ、、4ヶ月ぶりか。

 奥村先生のサイトに分割表の検定方法を比較しているページがありました。
Fisherの正確検定かカイ2乗検定か
 フィッシャーの正確検定か? ピアソンのカイ二乗か?
 巷の解説では、「フィッシャーの正確検定はパターンの数え上げに基づいているから正確だ。使用できる場合はいつでも使用すべき。」のような文言も見られますが、上のページを見れば正しくないと分かります!

 上のページでは、総度数だけ固定される場合と周辺度数が完全に固定される場合を、それぞれunconditional/conditionalとして検討しています。
 分割表の生成ルールは、どんな現象からどんな実験計画でデータをとったかによって変わります。

  • 度数が非固定の場合
  • 総度数が固定の場合
  • 各行の度数が固定の場合
  • 各行・各列の度数が固定の場合

 この4つについてシミュレーションを追試してみます。ソースの全体は記事の最後にあります。

度数が非固定の場合

 総度数含め度数が一切固定されない。

列1 列2
行1 行1合計
行2 行2合計
列1合計 列2合計 全合計

 一定時間内の観測データなどがあたりそうです。実験のたびに標本サイズnが揺らぐのが特徴です。
 セルの度数は独立なポアソン分布×4に従うと考えます。セルごとに平均(期待度数)のパラメータ \lambda_1 \cdots \lambda_4があります。

 \lambda_1  \lambda_2
 \lambda_3  \lambda_4

 パラメータ間に制約は設けないので、自由度は4。

行と列の独立条件

 独立とは、セルの度数に対して行と列が別個に効いてくるという意味です。
 ポアソン分布の平均パラメータのオッズ比が1になるように制限するのが独立モデルです。
 (「行列式が0」と等価なので、パラメータの行列が一次従属なら「独立」モデル、一次独立なら「非独立」モデルに……!)

 \frac{\lambda_1 \lambda_4}{\lambda_2 \lambda_3} = 1

 制約が入った分、自由に動かせるパラメータが減って、自由度は3。

総度数が固定の場合

 全セルの度数の合計が決まっている。

列1 列2
行1 行1合計
行2 行2合計
列1合計 列2合計 全合計

 決められた数の観測を行ったり、被験者数が事前に決まっているときはこれです。
 セルの度数は多項分布に従うと考えます。セルごとに確率パラメータがあって、n回試行した結果の分布です。

 \alpha_1  \alpha_2
 \alpha_3  \alpha_4

(ただし、 \alpha_1 + \alpha_2 + \alpha_3 + \alpha_4 = 1

 パラメータは見た目4つでも、3つ決まれば残り1つは自動的に決まり、自由度は3。

行と列の独立条件

 総度数が固定の場合の独立モデルは、多項分布の確率パラメータのオッズ比が1になるように制限します。
 セルの確率が行と列の確率の積の形に表せるということです。

 \frac{\alpha_1 \alpha_4}{\alpha_2 \alpha_3} = 1

 この制約で自由に動かせるパラメータが減って、自由度は2。

各行の度数が固定の場合

 各行にあるセルの度数の合計が決まっている。

列1 列2
行1 行1合計
行2 行2合計
列1合計 列2合計 全合計

 新薬の臨床試験で、新薬を投与するグループと従来薬を投与するグループの人数が決まっているときとか……。
 列を「改善なし/あり」とすると、各列の度数は固定されないと考えられます。

 セルの度数は行ごとに二項分布に従うとし、確率パラメータは同じとは限らない設定です。

  • 1行目
 p_1  1 - p_1
  • 2行目
 p_2  1 - p_2

 パラメータは p_1 p_2の2つ。パラメータ間の制約はないので、自由度は2。

行と列の独立条件

 列の確率がどちらの行でも同じなら独立といえます。
 つまり、1行目と2行目の二項分布の確率パラメータが等しいのが条件です。

 p_1 = p_2

 この制約で自由に動かせるパラメータが減って、自由度は1。

各行・各列の度数が固定の場合

 周辺度数(各行の度数、各列の度数、総度数)が完全に決まっている。

列1 列2
行1 行1合計
行2 行2合計
列1合計 列2合計 全合計

 奥村先生のページでconditionalと呼ばれているモデルです。
 セルの度数は左上のセルが超幾何分布に従って、残りのセルは自動的に決まります。

 推定が必要なパラメータはなく、自由度は0……と思いきや、黒木玄先生のツイートによると1ですか。。うむむ。

 ここから、シミュレーション開始です!

ピアソンのカイ二乗検定(補正なし)

 まず一番ポピュラーなピアソンのカイ二乗
 独立条件(下表の期待度数)で分割表を生成して、p値の一様性を評価します。点が斜め45度の線上に並ぶのが理想的な検定です。

45 15
30 10
# 1000回繰り返す
trials = 1000
pearson_p_values = numeric(trials)
for (i in 1:trials) {
  data = matrix(rtable_2x2(45,15,30,10,conditions_on="total"),
    ncol=2,byrow=T)
  pearson_p_values[i] = chisq.test(data,correct=F)$p.value
}

plot(ecdf(pearson_p_values),asp=1,xlab="", ylab="",
  xlim=c(0,1),ylim=c(0,1),verticals=TRUE,
  main="ECDF for chisq.test(correct=F)\n-- conditions_on: total")
abline(0,1)

f:id:zwxadz:20170924075122p:plain
f:id:zwxadz:20170924075134p:plain
f:id:zwxadz:20170924075143p:plain
f:id:zwxadz:20170924075151p:plain

 度数が揺らぐ場合や、各行の度数が固定される場合でも、ピアソンのカイ二乗は問題なく使用できそうです。
 周辺度数が完全に固定される場合だけ、上にずれています。棄却域を設定する0の近くで傾きが大きくなっていて危険率が上がっています。

イェーツの連続性補正あり

# 1000回繰り返す
trials = 1000
yates_p_values = numeric(trials)
for (i in 1:trials) {
  data = matrix(rtable_2x2(45,15,30,10,conditions_on="total"),
    ncol=2,byrow=T)
  yates_p_values[i] = chisq.test(data,correct=T)$p.value
}

plot(ecdf(yates_p_values),asp=1,xlab="", ylab="",
  xlim=c(0,1),ylim=c(0,1),verticals=TRUE,
  main="ECDF for chisq.test(correct=T)\n-- conditions_on: total")
abline(0,1)

f:id:zwxadz:20170924075309p:plain
f:id:zwxadz:20170924075314p:plain
f:id:zwxadz:20170924075317p:plain
f:id:zwxadz:20170924075321p:plain

 周辺度数が完全に固定される場合に線上に並んでいます。それ以外では、下にずれているので保守的(鈍感)に。
 「連続性補正(correction for continuity)」の呼び方は歴史的な誤解かもしれません。むしろ小標本の方が、ずれがひどい。。

フィッシャーの正確検定

# 1000回繰り返す
trials = 1000
fisher_p_values = numeric(trials)
for (i in 1:trials) {
  data = matrix(rtable_2x2(45,15,30,10,conditions_on="total"),
    ncol=2,byrow=T)
  fisher_p_values[i] = fisher.test(data)$p.value
}

plot(ecdf(fisher_p_values),asp=1,xlab="", ylab="",
  xlim=c(0,1),ylim=c(0,1),verticals=TRUE,
  main="ECDF for fisher.test()\n-- conditions_on: total")
abline(0,1)

f:id:zwxadz:20170924075352p:plain
f:id:zwxadz:20170924075355p:plain
f:id:zwxadz:20170924075358p:plain
f:id:zwxadz:20170924075401p:plain

 傾向はイェーツのカイ二乗と同じ感じです。この検定が正確なのは、周辺度数が完全に固定される場合ですね。この場合は、近似が入らないフィッシャーの方法が最適と考えられます。
 奥村先生のページに「(Fisherの方法を)近似するものとして連続性の補正をしたカイ2乗検定があります。」と書かれていました。

G検定

 対数線形モデルの尤度比検定だから、glm関数でいけそう……?(マニュアルにあたらないと、、)

ソース

 最後に、2×2分割表を生成するRのソースです。

# 2×2分割表
# 各セルの期待度数:a, b, c, d
# 度数の固定:conditions_on
#   ="none"	度数を固定しない
#   ="total"	総度数を固定(初期値)
#   ="rows"	各行の度数を固定(総度数も固定)
#   ="columns"	各列の度数を固定(総度数も固定)
#   ="rows.and.columns"	各行・各列の度数を固定(総度数も固定)

rtable_2x2 = function(a, b, c, d, conditions_on = "total") {
  if (min(a, b, c, d) < 0) {
    # 期待度数は非負でないとだめ
    print("Illegal argument!")
    return()
  }
  
  switch(conditions_on,
    "none" = return(rtable_2x2_n(a, b, c, d)),
    "total" = return(rtable_2x2_t(a, b, c, d)),
    "rows" = return(rtable_2x2_r(a, b, c, d)),
    "columns" = return(rtable_2x2_c(a, b, c, d)),
    "rows.and.columns" = return(rtable_2x2_rc(a, b, c, d)),
    {
      # その他の値はエラー
      print("Illegal argument!")
      return()
    }
  )
}

rtable_2x2_n = function(a, b, c, d) {
  # 乱数(ポアソン×4)
  counts = integer(4)
  counts[1] = rpois(1, lambda = a)
  counts[2] = rpois(1, lambda = b)
  counts[3] = rpois(1, lambda = c)
  counts[4] = rpois(1, lambda = d)
  return(counts)
}

rtable_2x2_t = function(a, b, c, d) {
  n = a + b + c + d	# 総度数
  # 乱数(多項)
  counts = rmultinom(1, size = n, prob = c(a, b, c, d) / n)
  return(counts)
}

rtable_2x2_r = function(a, b, c, d) {
  r1 = a + b	# 1行目の度数
  r2 = c + d	# 2行目の度数
  # 乱数(二項×2)
  counts = integer(4)
  counts[1] = rbinom(1, size = r1, prob = a / r1)
  counts[2] = as.integer(r1 - counts[1])
  counts[3] = rbinom(1, size = r2, prob = c / r2)
  counts[4] = as.integer(r2 - counts[3])
  return(counts)
}

rtable_2x2_c = function(a, b, c, d) {
  c1 = a + c	# 1列目の度数
  c2 = b + d	# 2列目の度数
  # 乱数(二項×2)
  counts = integer(4)
  counts[1] = rbinom(1, size = c1, prob = a / c1)
  counts[2] = rbinom(1, size = c2, prob = b / c2)
  counts[3] = as.integer(c1 - counts[1])
  counts[4] = as.integer(c2 - counts[2])
  return(counts)
}

rtable_2x2_rc = function(a, b, c, d) {
  r1 = a + b	# 1行目の度数
  r2 = c + d	# 2行目の度数
  c1 = a + c	# 1列目の度数
  # 乱数(超幾何)
  counts = integer(4)
  counts[1] = rhyper(1, m = r1, n = r2, k = c1)
  counts[2] = as.integer(r1 - counts[1])
  counts[3] = as.integer(c1 - counts[1])
  counts[4] = as.integer(r2 - counts[3])
  return(counts)
}
# 1000回繰り返す
trials = 1000
pearson_p_values = numeric(trials)
yates_p_values = numeric(trials)
fisher_p_values = numeric(trials)
for (i in 1:trials) {
  data = matrix(rtable_2x2(45,15,30,10,conditions_on="total"),
    ncol=2,byrow=T)
  pearson_p_values[i] = chisq.test(data,correct=F)$p.value
  yates_p_values[i] = chisq.test(data,correct=T)$p.value
  fisher_p_values[i] = fisher.test(data)$p.value
}

色温度からRGB(sRGB)に変換(完成版)

何をするの?

 こういうの見たことありませんか? 色温度に対応するRGB(sRGB)値を求めるコードを実装します。

f:id:zwxadz:20170502053830p:plain

計算した値の一部

 相関色温度4000K〜12000KのCIE昼光がどんな色をしているか計算しました。
 100K刻みで表にしましたが、実装コードではその間の値も計算できます。また、最高25000Kまで入力OKです。
 特定の色温度のカラーコード値を知るのに役立つかもしれません。

4000K: #FFD6A1 6000K: #FFF9F0 8000K: #E0ECFF 10000K: #CADDFF
4100K: #FFD9A6 6100K: #FFFBF3 8100K: #DFEBFF 10100K: #C9DCFF
4200K: #FFDBAB 6200K: #FFFCF6 8200K: #DDEAFF 10200K: #C8DCFF
4300K: #FFDDAF 6300K: #FFFDF9 8300K: #DCE9FF 10300K: #C7DBFF
4400K: #FFDFB4 6400K: #FFFEFC 8400K: #DBE8FF 10400K: #C6DBFF
4500K: #FFE2B9 6500K: #FFFFFF 8500K: #D9E7FF 10500K: #C6DAFF
4600K: #FFE4BD 6600K: #FCFEFF 8600K: #D8E6FF 10600K: #C5DAFF
4700K: #FFE6C1 6700K: #FAFCFF 8700K: #D7E5FF 10700K: #C4D9FF
4800K: #FFE7C5 6800K: #F7FAFF 8800K: #D5E5FF 10800K: #C3D9FF
4900K: #FFE9C9 6900K: #F5F9FF 8900K: #D4E4FF 10900K: #C3D8FF
5000K: #FFEBCD 7000K: #F3F7FF 9000K: #D3E3FF 11000K: #C2D8FF
5100K: #FFEDD1 7100K: #F1F6FF 9100K: #D2E2FF 11100K: #C1D7FF
5200K: #FFEED5 7200K: #EEF5FF 9200K: #D1E2FF 11200K: #C1D7FF
5300K: #FFF0D9 7300K: #ECF4FF 9300K: #D0E1FF 11300K: #C0D6FF
5400K: #FFF1DC 7400K: #EBF2FF 9400K: #CFE0FF 11400K: #C0D6FF
5500K: #FFF3E0 7500K: #E9F1FF 9500K: #CEE0FF 11500K: #BFD6FF
5600K: #FFF4E3 7600K: #E7F0FF 9600K: #CDDFFF 11600K: #BED5FF
5700K: #FFF6E7 7700K: #E5EFFF 9700K: #CCDEFF 11700K: #BED5FF
5800K: #FFF7EA 7800K: #E3EEFF 9800K: #CBDEFF 11800K: #BDD4FF
5900K: #FFF8ED 7900K: #E2EDFF 9900K: #CADDFF 11900K: #BDD4FF

参考にしたWeb上の情報

 各温度の黒体の色を求めています。等色関数を介してXYZ色空間にマッピングするので手間がかかりますが、広い範囲の温度について計算できます。

 色の用語、特にCIE関連について詳しい解説があります。

  • JIS Z 8720:2012

 表題は「測色用の標準イルミナント(標準の光)及び標準光源」。ググると法的に灰色(?)っぽいサイトがヒットして見られます。リンクは貼りません。

 sRGBへの変換はWikipediaを参考にしました。

前提条件

CIE昼光(daylight)

 黒体(black body)ではなく、CIE昼光です。昼光とほぼ等しい光と考えてください。
 JISでは、「主として、北空についての多くの分光測定値から統計的手法によって各々の相関色温度における昼光の代表としてCIEが定めた分光分布」です。昼光の測定データに基づいて定めたスペクトルの光ということですね。
 色温度というとき、普通は黒体、つまりプランクの法則に従う放射を考えます。しかし、黒体と異なるスペクトルをもつ光についても、知覚的に色の近い黒体の温度で表すことができ、これを相関色温度といいます。

白色点はD65

 約6500Kが白になります。より正確には約6504Kなのですが、最終的に8bit/channelで整数に丸めるので差が出ません。
 JIS Z 8720にあるxy色度座標への変換式を使用します。定義域の4000K〜25000Kから外れると計算できません。(というか、定義されない……?)

sRGB符号化

 だいぶ前からsRGBが標準になってると思うので。
 Webでよくある"#XXXXXX"表記はsRGB前提ですよね?(たぶん)

処理の大筋

  1. 色温度の範囲チェック
  2. 色温度→xy色度座標
  3. xyz→線形sRGB
  4. 線形sRGBを規格化
  5. 線形sRGB→非線形sRGB
  6. 非線形sRGBを量子化(0〜255の整数)
入力チェックとxy色度座標への変換(1と2)

 JIS Z 8720にある、CIE昼光のxy色度座標を求める式は次の通り。

 相関色温度T_{\rm cp}が4000K以上7000K以下の場合の色度座標x_{\rm D}

x_{\rm D} = -4.6070\,\frac{10^9}{{T_{\rm cp}}^3} + 2.9678\,\frac{10^6}{{T_{\rm cp}}^2} + 0.09911\,\frac{10^3}{T_{\rm cp}} + 0.244063

 相関色温度T_{\rm cp}が7000Kを超え25000K以下の場合の色度座標x_{\rm D}

x_{\rm D} = -2.0064\,\frac{10^9}{{T_{\rm cp}}^3} + 1.9018\,\frac{10^6}{{T_{\rm cp}}^2} + 0.24748\,\frac{10^3}{T_{\rm cp}} + 0.237040

 色度座標y_{\rm D}

y_{\rm D} = -3.000{x_{\rm D}}^2 + 2.870x_{\rm D} - 0.275

 4000K未満、25000K超えは、とりあえず例外投げておけばいいかな……。

線形sRGBへの変換(3)

 まず、色度座標z_{\rm D}を求めます。

z_{\rm D} = 1 - x_{\rm D} - y_{\rm D}

 ここでは色度(色合い)だけに興味があって、輝度は気にしないでいいので、xyzをそのままXYZとして計算します。(xyzはx+y+z = 1となるよう、XYZをスケーリングしたもの。)

\begin{bmatrix}
R_\text{linear} \\ G_\text{linear} \\ B_\text{linear}
\end{bmatrix} = \begin{bmatrix}
3.2406 & -1.5372 & -0.4986 \\ -0.9689 & 1.8758 & 0.0415 \\ 0.0557 & -0.2040 & 1.0570
\end{bmatrix} \begin{bmatrix}
X \\ Y \\ Z
\end{bmatrix}

線形sRGBを規格化(4)

 赤緑青の最大値が1となるようにスケーリングします。

\begin{bmatrix}
R_\text{scaled} \\ G_\text{scaled} \\ B_\text{scaled}
\end{bmatrix} = \frac{1}{\max \{R_\text{linear}, G_\text{linear}, B_\text{linear}\}} \begin{bmatrix}
R_\text{linear} \\ G_\text{linear} \\ B_\text{linear}
\end{bmatrix}

 なお、xyzからY = 1と仮定してXYZに変換、さらに線形sRGBに変換して1を超えた値をクリッピングする方法も考えられます。クリッピングによって色度が変わってしまうのを避けるため、スケーリングを選択しました。

非線形sRGBへの変換(5)

 いわゆるガンマ補正です。次式でチャンネル毎に変換します。

C_\text{sRGB} = \begin{cases}
12.92 C_\text{linear} & \text{if $C_\text{linear} \leq 0.0031308$,} \\
1.055 {C_\text{linear}}^{1/2.4} - 0.055 & \text{otherwise.}
\end{cases}

8bit/channelで量子化(6)

 最後に、0〜1の浮動小数点数を0〜255の整数に丸めます。四捨五入なり偶数丸めなりお好みの関数で。(コードでは四捨五入にしました。)

C_\text{int-sRGB} = \text{round} (255 C_\text{sRGB})

 実装コードはQiitaの記事へどうぞ!

qiita.com

三角形内の一様乱数

 頂点ABCの三角形内で一様な乱数はどうすればいいでしょうか?

A new method to simulate the triangular distribution
 前回と同じ文献。「Random vectors in the plane」の部分で触れられてます。この論文より先に知られていた方法らしく。

u \sim \text{Uniform}(0,1),\ v \sim \text{Uniform}(0,1);
\text{MIN} = \text{MIN}(u,v),\ \text{MAX} = \text{MAX}(u,v)として、
\vec{a} \text{MIN}+ \vec{b} (1 - \text{MAX}) + \vec{c} (\text{MAX} - \text{MIN})

 なぜこれでOKなのかよく分かりませんが、2次元の場合でプロットしてみると確かに三角形になる……。ふむ。

Point2D randomPoint(double x1, double y1, double x2, double y2, double x3, double y3) {
    double u = Math.random(), v = Math.random();
    
    double m1 = Math.min(u, v);
    double m2 = 1.0 - Math.max(u, v);
    double m3 = Math.max(u, v) - Math.min(u, v);
    double x = m1 * x1 + m2 * x2 + m3 * x3;
    double y = m1 * y1 + m2 * y2 + m3 * y3;
    return new Point2D.Double(x, y);
}

三角乱数

 マニアックな乱数記事を徐々に増やしていくキャンペーン。

 前回・前々回と、なぜ数学関数の使用を回避したいかを書いてなかったので、その補足から、、
 プログラムの数学関数、具体的にはサイン、コサイン、対数、平方根などですが、処理が重いことが多いです。除法も加・減・乗法と比べれば若干重いとはいえ、これの回避は難しい……。
 数学関数を使わないかわりに、単位円内の一様乱数を作るための繰返しが入りました。その繰返し回数を見積もってみます。
 円と外接正方形の面積の比は\pi/4 \approx 0.785 \ldotsで、円内に点が入ったときに脱出。その逆数が繰返し回数の平均になります。幾何分布なので数回程度かかる場合もありますが、平均は1回ちょっとです。
 ちなみに、実際にサインやコサインを使うより速いかは環境によりけりかも……?

 ここからが今日の本文。
 三角分布をPERTで使うという話は本当なのかどうか知りません。

Triangular distribution - Wikipedia
 三角乱数、平方根を使う方法がWikipediaに書いてあります。が、もっといいやり方(感動的!)があるので紹介します。

A new method to simulate the triangular distribution - ScienceDirect
u \sim \text{Uniform}(0,1),\ v \sim \text{Uniform}(0,1)として、
(1-c)\text{MIN}(u,v) + c\text{MAX}(u,v)
 これで、三角乱数(a=0, b=1, cは最頻値で0以上1以下)になります。MINとMAXは比較だけなので、平方根の計算よりずっと軽い!

double randomTriangular(double c) {
    double u = Math.random(), v = Math.random();
    
    return (1.0 - c) * Math.min(u, v) + c * Math.max(u, v);
}

 この方法が紹介されているのをあまり見ませんが、美しいアルゴリズムの一つだと思います。

コーシー乱数

 連日の更新なんていつぶり……?

 コーシー分布について。ちなみにコーシーと孔子は別人です。(豆知識
 コーシー分布は自由度1のt分布、平均も分散も発散するたちが悪い分布として名を知られてますね。
 和に対して安定(和の分布もコーシー分布)なのは正規分布と同じです。が、正規分布はスケールの二乗(分散)が加法的で、コーシー分布はスケールが加法的になります。

 分散1の独立な正規乱数の和は分散2の正規乱数になって、スケール(標準偏差)で見ると√2倍。さらに標本平均の分散は、
\text{Var}(\bar X) = \text{Var}(\frac{\sum X_i}{n}) = \frac{1}{n^2}\text{Var}(\sum X_i) = \frac{n\sigma^2}{n^2} = \frac{\sigma^2}{n}
 標本サイズが大きくなることで標本平均の分散が小さくなり、母平均の推定の信頼性が高まるのでした。

 片や、スケール1の独立なコーシー乱数の和はスケール2のコーシー乱数になります。n個の和をとって、スケールnのコーシー乱数。nで割るとスケール1のコーシー乱数。と、標本平均で元の分布に戻ってしまいました。母確率分布がコーシー分布のとき、標本サイズを大きくしても母中心の推定の信頼性は高まりません。

Cauchy distribution - Wikipedia
 コーシー乱数の生成には、逆関数法とか標準正規乱数の比をとるとかいろいろ見つかりますが、ここでは数学関数を使わない方法を、、

double randomCauchy() {
    double u, v;
    do {
        u = 2.0 * Math.random() - 1.0;
        v = Math.random();
    } while (u * u + v * v >= 1.0 || v == 0.0);
    
    return u / v;
}

 単位円内(上半分)からランダムに選んだ点で、X座標をY座標で割るとコーシー分布(0,1)の乱数が得られます。説明は割愛、Box–Mullerとかで検索するといいんじゃないかな……。