"+" を文字列連結演算子とする場合
"+" = function(e1, e2) {
if (is.numeric(e1) && is.numeric(e2)) {
base::"+"(e1, e2)
} else {
paste0(e1, e2)
}
}
> "abc" + "12345" + "あいうえお"
[1] "abc12345あいうえお"
> 123456 + "numeric"
[1] "123456numeric"
> 456 + 100
[1] 556
> 456L + 100L
[1] 556
> 456 + 0
[1] 456
> TRUE + FALSE
[1] "TRUEFALSE"
> "123" + TRUE + "asd"
[1] "123TRUEasd"
> (1 == 1) + (2 != 3)
[1] "TRUETRUE"
文字列と文字列の連結は，言語によって異なる記号が演算子として使われているが，R では二項演算子ではなく paste 関数，paste0 関数が使われる。
> paste("asd", "poi", sep="")
[1] "asdpoi"
> paste0("asd", "poi")
[1] "asdpoi"
> paste0(123, "asd")
[1] "123asd"
> paste0(123, TRUE)
[1] "123TRUE"
> paste0(123, T)
[1] "123TRUE"
> paste0(TRUE, FALSE)
[1] "TRUEFALSE"
TRUE と FALSE を連結するなんてことは普通は考えないこと。
& 記号を文字列としての連結演算子と定義する。ただし，2 項ともに論理型の場合には本来の & 演算子として使うというような関数を書いてみる。
"&" を "|" に変えても同じように動く。お好きな方を。
"&" = function(e1, e2) {
if (is.logical(e1) && is.logical(e2)) {
base::"&"(e1, e2)
} else {
paste0(e1, e2)
}
}
> "abc" & "12345" & "あいうえお"
[1] "abc12345あいうえお"
> 123456 & "numeric"
[1] "123456numeric"
> 456 & 100
[1] "456100"
> 456 & 0
[1] "4560"
> TRUE & FALSE
[1] FALSE
> "123" & TRUE & "asd"
[1] "123TRUEasd"
> (1 == 1) & (2 != 3)
[1] TRUE
超幾何分布
要因 a, b が独立な場合，以下のような分割表が得られる確率を求める
B not(B) sum
A x o m
not(A) (k-x) o n
sum k o (m+n)
超幾何分布により，求める確率は xCm * xCk-x / m+nCk
R の dhyper の引数はちょっと変で（？），x が生じる確率は dhyper(x, m, n, k)
B not(B) sum
A 2 o 10
not(A) (1) o 5
sum 3 o (15)
> x = 2; m = 10; n = 5; k = 3
> dhyper(x, m, n, k)
[1] 0.494505494505494
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] 0.494505494505495
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.494505494505495
しかし，
> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] NaN
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.00104601828356506
となり，普通のやり方では精度が足りない。
Rmpfr パッケージを用いて，単純に計算すると以下のようになる。
library(Rmpfr)
precBits=13
dhyper2 = function(x, m, n, k) {
choose2 = function(a, b) {
prod(mpfr(1:a, precBits)) / prod(mpfr(1:b, precBits)) / prod(mpfr(1:(a-b), precBits))
}
as.numeric(choose2(m, x) * choose2(n, k-x) / choose2(m+n, k))
}
> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> dhyper2(x, m, n, k)
[1] 0.00104601828356522
dhyper2 の実行時間は 0.15 秒
二項分布 B(x, n, p) = nCx p^x (1-p)^(n-x) の計算
R では，dbinom(x, n, p) で速く，簡単に求まる
これを，自分で計算してみるとどうなるかをやってみる。
式をそのまま記述すると以下のようになる。
options(digits = 15)
f = function(x, n, p) {
return(choose(n, x) * p ^ x * (1 - p) ^ (n - x))
}
n が小さいうちは，これでも十分である。
> x = 0
> n = 10
> p = 0.5
> dbinom(x, n, p)
[1] 0.0009765625
> f(x, n, p)
[1] 0.0009765625
しかし，n が大きくなると，これじゃダメだ。
> x = 4000
> n = 10000
> p = 0.4
> dbinom(x, n, p)
[1] 0.00814316030659453
> f(x, n, p)
[1] NaN
==============================================
一つのやりかたは，掛算・割り算を log の足し算・引き算で計算し，結果の exp をとる。
g = function(x, n, p) {
g0 = function(m) {
ifelse(m <= 1, 0, sum(log(2:m))) # R の仕様からくる問題を避けるために関数を定義
}
return(exp(g0(n) - g0(x) - g0(n - x) + x * log(p) + (n - x) * log(1 - p)))
}
> dbinom(x, n, p)
[1] 0.00814316030659453
> g(x, n, p)
[1] 0.00814316030660582
このやり方では，有効数字は 10 桁ほどである。
==============================================
もう一つのやり方は，n の階乗は n+1 に対するΓ関数なので，それを利用して対数Γ関数を使う。
h = function(x, n, p) {
return(exp(lgamma(n + 1) - lgamma(x + 1) - lgamma(n - x + 1) +
x * log(p) + (n - x) * log(1 - p)))
}
> dbinom(x, n, p)
[1] 0.00814316030659453
> h(x, n, p)
[1] 0.00814316030654657
このやり方では，有効数字は 11 桁ほどである。
==============================================
R では nCx の対数は lchoose 関数で計算できるので，少し簡単に書ける。
v = function(x, n, p) {
return(exp(lchoose(n, x) + x * log(p) + (n - x) * log(1 - p)))
}
> dbinom(x, n, p)
[1] 0.00814316030659453
> v(x, n, p)
[1] 0.00814316030660582
このやり方では，有効数字は 10 桁ほどである。
==============================================
nCx は整数値（約分すると分母は 1 になる）ということで，約分後の分子の log の和を求めて lchoose と同じ答えを目指す。
gcd = function(m, n) {
while ((temp <- n %% m) != 0) {
n <- m
m <- temp
}
return(m)
}
lchoose2 = function(n, k) {
k = min(k, n - k)
if (k == 0) {
return(0)
}
numerator = (n - k + 1):n
denominator = 1:k
for (den in seq_along(denominator)) {
for (num in seq_along(numerator)) {
GCD = gcd(denominator[den], numerator[num])
if (GCD != 1) {
denominator[den] = denominator[den] / GCD
numerator[num] = numerator[num] / GCD
if (denominator[den] == 1) {
break
}
}
}
numerator = numerator[numerator != 1]
}
return(sum(log(numerator)))
}
w = function(x, n, p) {
return(exp(lchoose2(n, x) + x * log(p) + (n - x) * log(1 - p)))
}
> dbinom(x, n, p)
[1] 0.00814316030659453
> w(x, n, p)
[1] 0.00814316030659842
このやり方では，有効数字は 12 桁に増加するが，計算時間が 5 秒近くかかる。
dbinom の結果の違いは，計算途中での計算精度が不足しているからである。
==============================================
R には多倍長演算をサポートするパッケージがいくつかあるが，そのうちの Rmfpr を使ってみる。
計算精度を適当（適切！）に指定して以下のような結果になる。
積・商を対数の和・差にする場合。
w2 = function(x, n, p) {
library(Rmpfr)
precBits = 70
P = mpfr(p, precBits)
lchoose = sum(log(mpfr((n-x+1):n, precBits))) - sum(log(mpfr(1:x, precBits)))
as.numeric(exp(lchoose + x * log(P) + (n - x) * log(1 - P)))
}
積・商をそのまま計算する場合。
w3 = function(x, n, p) {
library(Rmpfr)
precBits = 70
P = mpfr(p, precBits)
choose = prod(mpfr((n-x+1):n, precBits)) / prod(mpfr(1:x, precBits))
as.numeric(choose * (P ^ x) * (1 - P) ^ (n - x))
}
> dbinom(x, n, p)
[1] 0.00814316030659453
> w2(x, n, p)
[1] 0.00814316030659453
> w3(x, n, p)
[1] 0.00814316030659453
いずれでも， bdinom と同じ精度の結果が得られる。
計算時間は，積・商をそのまま計算する w3 の方が速く 0.14 秒ほどである（dbinom は瞬時）。
==============================================
二項分布 B(x, n, p) = nCx p^x (1-p)^(n-x) の計算の場合は dbinom があるので，以上のような事をやる必要はさらさらないが，「ものすごく大きな数とものすごく小さな数の積がほどほどの値になる場合」や「ものすごく大きな 2 つの数の商がほどほどの値になる場合」には，「積と商を対数の和と対数の差をとって結果の逆対数をとる」というのが必要になることもある。
そういうときには Rmfpr を使うのがもっとも効果的のようである。
... と，そういうこと。
factor 変数の場合
> fac = factor(c("c", "a", "b", "b", "a", "c", "b"))
> z = levels(fac)
> t(sapply(fac, function(x) z %in% x))[, -1] + 0
[,1] [,2]
[1,] 0 1
[2,] 0 0
[3,] 1 0
[4,] 1 0
[5,] 0 0
[6,] 0 1
[7,] 1 0
数値変数の場合
> num = c(1, 3, 2, 1, 3, 3, 2)
> z = as.numeric(names(table(num)))
> t(sapply(num, function(x) z %in% x))[, -1] + 0 # factor 変数の場合と同じ
[,1] [,2]
[1,] 0 0
[2,] 0 1
[3,] 1 0
[4,] 0 0
[5,] 0 1
[6,] 0 1
[7,] 1 0
