餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Bloggerで独自ドメインを利用している場合のSSL/TLS化

ここ5年ぐらい代表的ブラウザーGoogle Chromeが非SSL/TLSであるhttp://のページに警告を出すようになっていて、Blogger独自ドメイン www.anlyznews.com で書いているブログのほうもLet's Encryptを使ってSSL/TLS化したのでメモ書きを。

ブログごとの設定のHTTPSの使用をオンにすれば、秘密鍵や証明書の作成と管理はBloggerがやってくれるのですが、オンにする前にDNSにCAAレコードを1つ追加する必要があります。ホスト名:anlyznews.com,VALUE:letsencrypt.org,FLAGS:0,TAG:issueで追加しました*1DNSのレコード追加後にHTTPの使用をオンにすると、数時間後にはSSL/TLS化が行えます。なお、証明局はLet's Encryptである必要はなく、例えばGoogle Trust Servicesを使うときはVALUEpki.googにすればよいそうです。

ところでBloggerは2010年から利用しているのですが、独自ドメインを使うときのDNSまわりの設定が当時と大きく変わっていて、今はばらばらAレコードを4つ追加する必要はなく、CNAMEレコードを2つ追加すればよくなっていました。ブログCNAMEとセキュリティCNAMEで、セキュリティCNAMEはカスタムドメインの利用開始時に表示されます。HTTPSの使用をオンにするときにエラーが出るので、もしやと思って今風のDNS設定に変えたのですが、カスタムドメインを一度削除したあと、再設定することになりました*2

*1: TAGがissuewildのレコードも追加したのですが、発行された証明書を眺める限り、ワイルドカード証明書になっていなかったので、こちらは参照されていません。

*2:今風のDNS設定後、無事にHTTPSの使用をオンにできましたが、関連は分かりません。

粒子フィルタ(particle filter)で状態とパラメーターを同時推定

ノイズが多い時系列データに対しては、観測値から観測誤差を取り除いて、分析者の興味関心である観察不能な状態を取り出すカルマンフィルターがよく使われています。しかし、カルマンフィルタは線形近似になるため、より一般的に用いることができるアルゴリズムとして粒子フィルタが開発されていて、カルマンフィルタほどではないですが様々な分野で利用されています。

矢野 (2014)「粒子フィルタの基礎と応用:フィルタ・平滑化・パラメータ推定」によると、粒子フィルタは1980年代に原案があり、1990年代に提案された相対的に新しい手法です。検索をすると、ここ5年間でも手法の改善を提案する論文があります。一方で、計算手順は簡素で実装しやすいものです。計算時間や記憶装置がカルマンフィルタよりも多く必要になりますが、イマドキの計算機にとっては誤差みたいなもんになっています*1

さて、プログラマ的にはコードを走らせて感覚を掴みたいわけですが、目についた例は簡素すぎました。粒子フィルタはシステム誤差を乱数で生成し、観測誤差を尤度関数で測って推定を行なうのですが、例ではこの2種類の乱数の確率密度関数のパラメーターが定数で与えられています。しかし、事前にシステム誤差と観測誤差のパラメーターを知ることはできないので、応用によっては実用的に思えないものになっています。

粒子フィルタではパラメーターの状態を考える事により、パラメーターの推定も同時に行なう事ができます。実用ではパッケージを使う方がよいですが、何をやっているか把握してみましょう。状態変数はコーシー分布か一様分布に従うシステムノイズによって変化し、観測値には正規分布に従う誤差が入る、状態変数が X, \tau, \thetaの3つ、観測値が Yの1つの推定モデルを考えて推定します。

1. データセット

Kitagawa (1998)の最初の例なのですが、不連続な状態遷移と正規分布する観測誤差を仮定してデータを生成します。なお、推定モデルのシステムノイズの仮定には沿っていません。

set.seed(406149)

# t: 観察不能な状態
# T: 分析期間の長さ
# y: 観測値
t <- rep(c(0, 1.5, -1, 0), each = 100)
T <- length(t)
y <- rnorm(T, mean = t, sd = 1)

2. 粒子フィルタで推定

計算自体は50行かからないぐらいで終わります。パッケージで実装されているものは、もっと複雑なことをしていると思いますが。なお、パラメーターの状態変数のシステムノイズが変則的に入れてありますが、パラメーターが0以下にならないようにです。

# 仮定するsystem noiseのνはコーシー分布
nu <- function(location, log_tau2) rcauchy(length(log_tau2), location = location, scale = sqrt(10^log_tau2))
# 観測誤差は正規分布
calcWeights <- function(y, x, sd){
    p <- dnorm(y, mean = x, sd = sd)
    p / sum(p)
}

# M: 各期でサンプリングする数
# x: 各期でサンプリングした値
# log_tau2: 各期でサンプリングに使う分布のパラメーター/x[i,j]ごとにある/log(tau^2, 10)
# w: 各期でリサンプリングに使うウェイト/尤度
# sd: 観測誤差の標準偏差

M <- 201

# 初期値を設定
w <- theta <- log_tau2 <- x <- matrix(NA, M, T)
log_tau2[, 1] <- runif(M, min = -8, max = -2)
theta[, 1] <- runif(M, min = 1e-6, max = 2)
x[, 1] <- rnorm(M, sd = 4)
# 始点のウェイトは雑につくる
w[, 1] <- calcWeights(y[1], x[, 1], sd = theta[, 1])

# パラメーターの誤差を大きくして、縮退しないようにする関数
g <- function(x) x*runif(M, min = 0.9, max = 1.1)

# フィルタリング分布を計算
for(j in seq(2, T)){
    # 予測分布をつくる
    x[, j] <- nu(x[, j - 1], log_tau2[, j - 1])

    # リサンプリング用のウェイトをつくる
    w[, j] <- calcWeights(y[j], x[, j], sd = theta[, j - 1])

    # リサンプリング
    i <- sample.int(M, replace = TRUE, prob = w[, j])
    x[, j] <- x[i, j]
    # パラメーターは縮退防止で粒子のばらつきを増している
    log_tau2[, j] <- g(log_tau2[i, j - 1])
    theta[, j] <- g(theta[i, j - 1])
}

3. プロット

推定結果を図示してみましょう。

3.1. フィルタリング分布

完璧とは言えませんが、観測誤差がだいぶ取れているのが分かります。

x_mean <- apply(x, 2, mean)
x_quantiles <- apply(x, 2, function(x) quantile(x, probs = c(0.025, 0.975)))

par(cex = 1, cex.lab = 1, cex.axis = 1, cex.main = 1, mar = c(5, 3, 1, 1))
plot(y, xlab = "t", ylab = "", ylim = c(-4, 4))
lines(x_mean, col = "blue", lwd = 2)
polygon(c(1:T, T:1), c(x_quantiles[1, ], rev(x_quantiles[2, ])), density=20, col="blue", border=NA)
lines(t, col = "red", lwd = 2)
legend("bottom", c("observations", "true states", "mean of the filtering", "95% confidence intervals"),
    ncol = 2, bty = "n", inset = 0.00, cex = 1.25,
    lty = c(NA, 1, 1, NA), lwd = c(1, 2, 2, NA), pch = c(21, NA, NA, NA), 
    fill = c(NA, NA, NA, "blue"), density = c(0, 0, 0, 20), border = NA,
    col = c("black", "red", "blue", "blue"))


3.2. 平滑化分布

構造的にはフィルタリング分布と同じなのですが、過去の粒子も推定に用いることで滑らかにします。

# 平滑化分布を計算
L <- 20
s <- matrix(NA, M*L, T)
for(j in 1:T){
    columns <- max(1, j - L + 1):j
    i <- sample.int(M, replace = TRUE, prob = w[, j])
    tmp <- sapply(columns, function(k){
        x[i, k]
    })
    s[1:length(tmp), j] <- tmp
}

s_mean <- apply(s, 2, mean, na.rm = TRUE)
s_quantiles <- apply(s, 2, function(x) quantile(x, probs = c(0.025, 0.975), na.rm = TRUE))

par(cex = 1, cex.lab = 1, cex.axis = 1, cex.main = 1, mar = c(5, 3, 1, 1))
plot(y, xlab = "t", ylab = "", ylim = c(-4, 4))
lines(s_mean, col = "darkgreen", lwd = 2)
polygon(c(1:T, T:1), c(s_quantiles[1, ], rev(s_quantiles[2, ])), density=20, col="darkgreen", border=NA)
lines(t, col = "red", lwd = 2)
legend("bottom", c("observations", "true states", "mean of the smoothing", "95% confidence intervals"),
    ncol = 2, bty = "n", inset = 0.00, cex = 1.25,
    lty = c(NA, 1, 1, NA), lwd = c(1, 2, 2, NA), pch = c(21, NA, NA, NA), 
    fill = c(NA, NA, NA, "darkgreen"), density = c(0, 0, 0, 20), border = NA,
    col = c("black", "red", "darkgreen", "darkgreen"))


4. カルマンフィルターとの比較

以下はKFASパッケージで計算したものですが、上で雑に計算した粒子フィルタの方が予測誤差は大きく、偏りが少なくなっています。この例だとカルマンフィルターに対する優位性は恐らくないです。

*1:組み込み用途ではつらいかも知れません。

Rで計量経済学の変量効果モデルを推定

だいぶ前にパッケージを使わず固定効果モデルで推定する例は書いていたのですが、変量効果モデルで推定する例も書いてみました。計量経済学で言う変量効果モデルは、観測値ごとの誤差 u_{it}に加えて、個体ごとの誤差c_{i}がある線形回帰モデル y_{it} = X\beta + c_{i} + u_{it}です。yは被説明変数、Xは説明変数、\betaは係数、iは個体を示す添字、tは時点を示す添字になります。実用ではパッケージを使うべきですが、説明でもパッケージを使うとブラックボックス感が出てしまうので。

1.データセット

定番テキストのWooldridgeではjtrainデータセット*1が例になっていたりするのですが、実データだと真のDGPとそのパラメーターが分からないので、乱数から生成します。

mkdata <- function(n_of_t, n_of_i){
    # n_of_t: 時点の数, n_of_i: 個体の数

    # 個体番号をつくる
    id <- rep(1:n_of_i, each = n_of_t, min=-0.5, max=0.5)

    # 時点変数をつくる
    t <- rep(1:n_of_t, n_of_i)

    # 説明変数x,zをつくる
    x <- runif(n_of_i*n_of_t)
    z <- runif(n_of_i*n_of_t)

    # 誤差項をつくる
    e_ind <- rep(rnorm(n_of_i, sd = 2), each = n_of_t)
    e_ids <- rnorm(n_of_i*n_of_t)
    e <- e_ind + e_ids

    # 従属変数yをつくる
    y <- 1 + 2*x - 3*z + e

    data.frame(id = as.factor(id), t = as.factor(t), y, x, z)
}

set.seed(1303)
df01 <- mkdata(5, 100)

2. 行列演算などで推定

定番テキストのWooldridgeのGSLとRandom Effect Modelの説明に沿って計算します。OLSの推定結果から、個体の分散 \sigma^2_cと観測ごとの分散 \sigma^2_uを推定して、ウェイト行列 \SigmaをつくってGLSにかけます。

 \Sigma = \begin{bmatrix}
\sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c \\
\sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c \\
\sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c \\
\sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c \\
\sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c \\
\end{bmatrix}

 \hat{\beta} = (^t\!X (I_n \otimes \Sigma^{-1}) X)^{-1} {^t\!X} (I_n \otimes \Sigma^{-1}) y

自由度がややこしいので注意してください。

frml <- y ~ x + z
X <- model.matrix(frml, data = df01)
mf <- model.frame(frml, data = df01)
y <- model.response(mf)
beta_ols <- solve(t(X)%*%X) %*% t(X) %*% y
residuals <- X %*% beta_ols - y
df <- nrow(X) - ncol(X)
sigma2_ols <- sum(residuals^2)/df
T <- length(levels(df01$t))
n <- length(levels(df01$id))

# total error variance
sigma2_v <- sum(residuals^2)/df

# individual error variance
sigma2_c <- with(df01, {
    sum(sapply(levels(id), function(i){
        r <- residuals[id == i]
        cp <- outer(r, r)
        sum(cp[upper.tri(cp, diag = FALSE)])
        # 以下と同じ結果
        # sum <- 0
        # for(t in 1:(T-1)){
        #     for(k in (t+1):T){
        #         sum <- sum + r[t]*r[k]
        #     }
        # }
        # sum
    }))/((n*T*(T-1)/2) - ncol(X))
})

# idiosyncratic error variance
sigma2_u <- sigma2_v - sigma2_c

omega <- matrix(sigma2_c, T, T) + sigma2_u*diag(T)
iomega <- diag(n) %x% solve(omega)
beta_gls <- solve(t(X) %*% iomega %*% X) %*% t(X) %*% iomega %*% y

# Usual Variance Matrix Estimator
usual_vcov <- solve(t(X) %*% iomega %*% X)

# Robust Variance Matrix Estimator
rbst_vcov <- with(df01, {
    dfs <- split(df01, id)
    so <- solve(omega)
    A <- matrix(0, ncol(X), ncol(X))
    B <- matrix(0, ncol(X), ncol(X))
    n <- length(dfs)
    for(i in 1:n){
        X <- model.matrix(frml, data = dfs[[i]])
        y <- dfs[[i]]$y
        u <- y - X %*% beta_gls
        A <- A + t(X) %*% so %*% X
        B <- B + t(X) %*% so %*% u %*% t(u) %*% so %*% X
    }
    solve(A) %*% B %*% solve(A)
})

3. パッケージで推定

比較のために、plmパッケージで推定します。

library(plm)
pd01 <- pdata.frame(df01, index=c("id"))
r_plm <- plm(frml, model = "random", data = pd01)

4. 推定結果を比較

係数の推定量と標準誤差を比較してみましょう。

beta <- matrix(
    c(1, 2, -3, beta_ols, beta_gls, coef(summary(r_plm))[, 1]),
    ncol(X), dimnames = list(rownames(beta_gls), c("true", "OLS", "GLS", "plm/random")))

sigma <- matrix(
    c(sqrt(sigma2_ols*diag(solve(t(X) %*% X))), sqrt(diag(usual_vcov)), sqrt(diag(rbst_vcov)), coef(summary(r_plm))[, 2]),
    ncol(X), dimnames = list(rownames(beta_gls), c("OLS", "GLS/Usual", "GLS/Robust", "plm/random")))

係数が

true OLS FGLS plm/random
(Intercept)  1.0000  0.7393  1.1459  1.1471
x  2.0000  2.6266  2.0862  2.0845
z -3.0000 -2.7875 -3.0411 -3.0419

標準誤差が

OLS FGLS/Usual FGLS/Robust plm/random
(Intercept)  0.28772  0.25141  0.25710  0.25087
x  0.36308  0.18589  0.20435  0.18237
z  0.35117  0.17542  0.15863  0.17208

となります。FGLSで推定した係数は、OLSよりも真の値に近くなっており、係数も標準誤差もplmパッケージとほとんど同じです。

ほとんど同じと言うところに引っかかると思いますが、plmパッケージでは以下の連立方程式を解いて個体内変化と個体間相違を計算して、idiosyncratic error varianceとindividual error varianceにしており、分散分解の方法が教科書のものと異なっているためです。

  • {withinの残差平方和} = {個体内変化}*((T-1)*N-K+1) + {個体間相違}*0
  • {betweenの残差平方和} = {個体内変化}*(N-K)/T + {個体間相違}*(N-K)

*1:Rだとwooldridgeパッケージの中にあります。

miceの並列処理で時間短縮

欠測処理でRのmiceパッケージのmice()を使うとして、処理時間の問題皆さんどうしてますかね。」と言う質問に、「多重代入法はM個の代入を独立して実行するのが望ましいので、マルチコアを使える演算環境ならば、並列処理することができます…Flexible Imputation of Missing Data第2版の5.5節Parallel computationが参考になると思います。」と言うリプライがあって、参照先を見ると

  • doParallelパッケージのcomplete関数とibind関数を組み合わせて何とかする
  • parlMICE関数を使う
  • micemdパッケージのpar.mice関数を使う(追記:同様に使えるfuturemice関数の方が使いやすいそうです)

方法が示唆されていたのですが、一番手軽そうなparlMICEを使う方法の検索したら最初に出てきた記事の内容が古くなっていたので、例をつくってみました。

library(mvtnorm)
library(mice)

# 練習用のデータセットをつくる
mdata <- with(list(n = 100,
    mu = rep(0, 4),
    sigma = matrix(rep(c(1.0, rep(0.5, 4)), 4)[1:16], 4, 4),
    beta = matrix(c(1, -2, 3, -4), 4, 1)), {
    data <- rmvnorm(n, mean = mu, sigma = sigma)
    colnames(data) <- paste("x", 1:ncol(data), sep = "")
    y <- -1 + data %*% beta + rnorm(n, sd = 2)
    nhanes <- ampute(data, prop = 0.8)
    cbind(y, nhanes$amp)
})

# 並行処理のための準備
# n.core: 平行処理で利用するコア数
# n.imp.core: コアごとの生成データセット数
# methodやpredictorMatrixはmice関数と同様(blocksは上手く動かない?)
imp <- parlmice(mdata, n.imp.core = 20, cluster.seed = 223, n.core = 3)
# もしくは
# imp <- futuremice(mdata, n.imp.core = 20, cluster.seed = 223, n.core = 3)
imp$method

# 欠損値補定
fit <- with(imp, lm(y ~ x1 + x2 + x3 + x4))
pool(fit)

普段、miceを使わないので上手く動いているか良く分かりません。

コーシー分布も扱えるランク基準線形モデル推定(Rank-based Estimation for Linear Models)

Kloke and McKean (2012) "Rfit: Rank-based Estimation for Linear Models," The R Journal, Vol.4(2)にアルゴリズムの概要と利用方法の詳細が書いてあるので説明する必要はないのですが、U検定やBrunner-Munzel検定を使ってきた外れ値がある観測値がコーシー分布らしいデータに線形回帰や分散分析をかけるのに便利なパッケージがありました。

ランク基準線形モデル推定は、ランク*1を用いた線形回帰で一致推定量を計算する手法で、OLSでは一致推定量が得られない誤差項がコーシー分布である場合なども、頑強な結果を示してくれます。ノンパラ手法の紹介ページで存在を知ったのですが、カーネル回帰のように大きなサンプルサイズが必要になるわけではなく、U検定やBrunner-Munzel検定を使う状況で、統制変数を加えてみたいときなどに気軽に使えそうです。

誤差項がコーシー分布

誤差項がコーシー分布の場合は、OLSよりも真のパラメーター(i.e. 切片項が1,xの係数が2,zの係数が-3)に近い推定量を出してきます。

# データセットを作成
df01 <- with(list(n = 90), {
	set.seed(520)
	x <- rnorm(n, sd = 2)
	z <- runif(n, min = -1, max = 1)
	g <- 1:n %% 3
	gcoef <- c(0, 0.05, -0.05)
	y_c <- 1 + 2*x - 3*z  + rcauchy(n, scale = 1)
	y_n <- 1 + 2*x - 3*z  + rnorm(n)
	data.frame(y_c, y_n, x, z, g = as.factor(g))
})

# 誤差項がコーシー分布のときの推定
library("Rfit")
r_rfit_c <- rfit(y_c ~ x + z + g, data = df01) # 推定
summary(r_rfit_c) # 推定結果の確認
summary(lm(y_c ~ x + z + g, data = df01)) # OLSの場合


誤差項が正規分布

誤差項が正規分布の場合は、OLSとだいたい同じ結果です。

r_rfit_n <- rfit(y_n ~ x + z + g, data = df01) # 推定
summary(r_rfit_n) # 推定結果の確認
summary(lm(y_n ~ x + z + g, data = df01)) # OLSの場合


モデル選択

モデル選択はlmとanovと同様にできます。

drop.test(r_rfit_c, rfit(y_c ~ x + z, data = df01)) # モデル選択


分散分析

一元配置分散分析はoneway.rfitで、n元配置分散分析はraovでできます。

r_oneway_anova <- with(df01, oneway.rfit(y, g)) # 一元配置分散分析
summary(r_oneway_anova, method = "tukey") # 多重比較の補正をして結果表示

ぽちっとな感がある実用的なソリューションですね。

*1:誤差項のベクトルをRのrank関数で計算した値。

尖度で自由度を補正したF検定による等分散性の検証

教科書の分散比のF検定をつかった等分散性の検証は、F検定が正規分布に依存しているので利用できる場合が限られるわけですが、それは分布の裾野の厚さ、尖度が問題であって、尖度が分かればそれで自由度を補正してより広い分布に使えることが知られています。

理屈

O'Neill (2014)のp.294の S^2_n/\sigma^2 \sim \chi^2(DF_n)から、 \frac{S^2_m/\sigma^2}{S^2_n/\sigma^2}  \sim \frac{\chi^2(DF_m)}{\chi^2(DF_n)}=F(DF_m,DF_n)であることが直ちに分かります。ここで、 n,mはサブサンプルの観測数で、 DF_n, DF_mは尖度による補正をした自由度です。 DF_n = 2n / (\kappa - (n-3)/(n-1))であり、 DF_mも同様です。 \sigma, \kappa帰無仮説の上で比較する2つのサンプルにとって共通の母集団の標準偏差と尖度になります。

尖度 \kappaの推定が問題になりますが、比較する2つの標本の平均がそれぞれゼロになるように、それぞれの標本の観測値をシフトした上で、2つの標本をプールした上で尖度の不偏推定量を計算します。こうすることで、平均がシフトしているだけで分散が同じ分布の比較で、帰無仮説を棄却しないようになります。

なお、参考にした元ネタではO'Neill (2014)のResult 15を参照しろと言う話になっていたのですが、Result 15はサンプルとサブサンプルの分散を比較する話で、サンプルとサンプルの分散を比較する話ではないので、そのまま利用しませんでした。

コード

kurtosis <- function(...){
    # 尖度の不偏推定量を計算(注意:正規分布の尖度が0ではなく3になる方)
    args <- list(...)
    kappa <- numeric(length(args))
    for(i in 1:length(args)){
        a <- args[[i]]
        n <- length(a)
        m <- mean(a)
        kappa[i] <- n * (n + 1) / (n - 1) * sum((a - m)^4) / sum((a - m)^2)^2
    }
    kappa
}

on2014ss2ss <- function(a, b, alternative = NULL){
    # サブサンプルとサブサンプルの分散比の統計量を、自由度補正してF検定
    n <- c(length(a), length(b))
    s2 <- c(var(a), var(b))
    kappa <- kurtosis(c(a - mean(a), b - mean(b)))
    DF_n <- function(n) 2*n / (kappa - (n - 3)/(n - 1))
    cpf <- function(i, j) pf(s2[i]/s2[j], DF_n(n[i]), DF_n(n[j]))

    if(is.null(alternative)) alternative <- "two.sided"
    (switch(alternative, "two.sided" = function(){
        2*min(cpf(1, 2), cpf(2, 1))
    }, "less" = function(){
        cpf(1, 2)
    }, "greater" = function(){
        cpf(2, 1)
    }))()
}

使い方は簡単で、以下のように二つのサンプルのベクトルを引数に与えるだけです。

n <- 30
a <- rgamma(n, shape = 2, scale = 3)
b <- rgamma(n, shape = 2, scale = 3)
on2014ss2ss(a, b)

シミュレーション

この方法(O’Neill)がどの程度有用かシミュレーションして確かめてみました。教科書的なF検定と、非正規分布でも信頼性が高いといわれるLevene検定と比較します。

まずは偽陽性が過剰にでないか、有意水準5%で同一の正規分布、一様分布、ラプラス分布、ガンマ分布、平均値だけシフトしたラプラス分布、平均値だけシフトした一様分布から生成した2つの乱数列で検定を行なってみました。

norm laplace laplace/shift unif unif/shift gamma
F-test 0.041 0.132 0.137 0.006 0.006 0.114
O’Neill 0.041 0.058 0.053 0.043 0.052 0.046
Levene 0.043 0.046 0.043 0.043 0.044 0.044

0.05が望ましい値ですが、教科書的なF-testは裾野が厚い一様分布では検出力が過小で、裾野が薄いラプラス分布では件出力が過剰になっています。今回紹介した方法とLeveneはこのような問題はありません。

2つの乱数列を生成した分布のパラメーター(と分散)が異なる場合も試してみましょう。

norm laplace unif gamma
F-test 0.833 0.681 0.553 0.966
O’Neill 0.824 0.536 0.797 0.914
Levene 0.661 0.440 0.449 0.895

こちらは高い方が望ましいと言うことになりますが、今回紹介した方法はLevene法よりも高くなっています。なお、教科書的なF-testは高いのですが、偽陽性も高いので頼りにはなりません。

なお、シミュレーションに使ったRのコードは以下になります。

library(VGAM) # ラプラス分布用
library(car) # Levene検定用
set.seed(547)
nr <- 1000 # 母集団分布×手法ごとの試行回数
n <- c(70, 30) # 試行ごとのサンプルサイズ
names_of_methods <- c("F-test", "O’Neill", "levene")

calcPs <- function(a, b){
    pv <- numeric(length(names_of_methods))
    n <- c(length(a), length(b))

    df01 <- data.frame(y = c(a, b), g = as.factor(c(rep("a", n[1]), rep("b", n[2]))))
    r <- leveneTest(y ~ g, data=df01)
    pv[3] <- r[[3]][1]

    # Old fashion F-test
    s2 <- c(var(a), var(b))
    df <- n - 1
    qf <- s2[1]/s2[2]
    pv[1] <- pf(qf, df[1], df[2])

    # Newbie based on O'Neill (2014)
    pv[2] <- on2014ss2ss(a, b)
    pv
}

simulate <- function(dgp1, dgp2, n, nr = 1000, alpha = 0.05){
    pv <- matrix(NA, nr, length(names_of_methods), dimnames = list(NULL, names_of_methods))
    for(i in 1:nr){
        pv[i, ] <- calcPs(dgp1(n[1]), dgp2(n[2]))
    }
    apply(pv, 2, function(x) sum(x < alpha))
}

r_fp <- sapply(c(
    function(n) rnorm(n, mean = 0, sd = 1),
    function(n) runif(n, min = 0, max = 1),
    function(n) rlaplace(n, location = 0, scale = 1),
    function(n) rgamma(n, shape = 2, scale = 3)
), function(dgp, n, nr){
    simulate(dgp, dgp, n, nr)
}, n = n, nr = nr)

r_fp <- cbind(r_fp, simulate(
    function(n) rlaplace(n, location = 0, scale = 1),
    function(n) rlaplace(n, location = 10, scale = 1),
    n = n, nr = nr
), simulate(
    function(n) runif(n, min = 0, max = 1),
    function(n) runif(n, min = 20, max = 21),
    n = n, nr = nr
))

r_fp <-  r_fp / nr
colnames(r_fp) <- c("norm", "unif", "laplace", "gamma", "laplace/shift", "unif/shift")
r_fp <- r_fp[, c(1, 3, 5, 2, 6, 4)]

r_pp <- cbind(simulate(
    function(n) rnorm(n, mean = 0, sd = 1),
    function(n) rnorm(n, mean = 0, sd = 1.5),
    n = n, nr = 1000
),simulate(
    function(n) rlaplace(n, location = 0, scale = 1),
    function(n) rlaplace(n, location = 0, scale = 1.5),
    n = n, nr = 1000
), simulate(
    function(n) rgamma(n, shape = 2, scale = 3),
    function(n) rgamma(n, shape = 2, scale = 6),
    n = n, nr = 1000
), simulate(
    function(n) runif(n, min = 0, max = 1),
    function(n) runif(n, min = -0.15, max = 1.15),
    n = n, nr = 1000
))

r_pp <- r_pp / nr
colnames(r_pp) <- c("norm", "laplace", "gamma", "unif")
rownames(r_pp) <- names_of_methods
r_pp <- r_pp[, c(1, 2, 4, 3)]

print("偽陽性率")
print(r_fp)

print("真陽性率(違う分布の値は比較不可)")
print(r_pp)

RでBEEP音を鳴らす方法

BEEP音を鳴らしたかったのですが、検索するとbeeprパッケージを使え、alarm関数を使えと言うコレじゃないソリューションや、system関数でOSのコマンドを叩けと言うベタな手が引っかかります。どうやらRでBEEP音を鳴らしたい人々はほとんどいないようです。

Cの標準関数でBEEP音がサポートされていないせいか、Rは標準ではBEEP音を鳴らせません。その他の音も鳴らせません。パッケージの力を借りる事になるのですが、audioパッケージのplay関数が正弦波を音に変えてくれるので、BEEP音の代わりに使えます。

早速、どれみふぁそらしど…とベタに鳴らしてみましょう。

library(audio)

# 12平均律の音階周波数を生成
mkHz <- function(h = 4){
    Hz <- 440*(2^(1/12))^{(-12*4):(12*h)}
    n <- length(Hz)
    scale <- c("C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B")
    j <- 1
    lst <- list()
    while(0 < n){
        lst[[j]] <- paste(scale[1:min(n, 12)], j, sep="")
        j <- j + 1
        n <- n - 12
    }
    names(Hz) <- unlist(lst)
    Hz
}

# 音階と音程と音量から正弦波を生成
mkWave <- function(..., length, volume, rate){
    if(0 >= ...length()) stop("no scale inputed!")
    Hz <- mkHz()
    lst <- list()
    for(i in 1:length(...elt(1))){
        sinc <- 0
        for(j in 1:...length()){
            s <- ...elt(j)[i]
            if(!is.na(Hz[s])){
                sinc <- sinc + sin(1:rate/rate*(2*pi)*Hz[s])
            }
        }
        sinc <- sinc / ...length()
        lst[[i]] <- (sinc * volume[i])[1:(length[i])]
    }
    unlist(lst)
}

rate <- 44100 # サンプリングレート
scale <- c("C5", "D5", "E5", "F5", "G5", "A5", "B5", "C6")
length <- rep(rate/4, length(scale))
volume <- rep(1/4, length(scale))

# playで音をならし、waitで鳴らし終わるのを待つ
wait(play(mkWave(scale, length = length, volume = volume, rate = rate), rate))

正弦波を加算すれば和音(コード)になりますし、原理的には色々とできるわけですが、MIDIの再発明みたいな事になるのでこの辺でやめておこうかと思っていたのですが、正解/不正解/挨拶の効果音を用意しました。

correct <- function(rate = 44100){
    scale <- c("E5", "C5", "E5", "C5", "E5", "C5")
    length <- c(rep(rate/8, length(scale)-1), rate/4)
    volume <- rep(1/4, length(scale))
    wait(play(mkWave(scale, length = length, volume = volume, rate = rate), rate))
}

wrong <- function(rate = 44100){
    scale1 <- c("F#3", "F#3")
    scale2 <- c("G2", "G2")
    length <- c(rate/3, rate/1.5)
    volume <- rep(2/3, length(scale1))
    wait(play(mkWave(scale1, scale2, length = length, volume = volume, rate = rate), rate))
}

bow <- function(rate = 44100){
    length <- c(rate, rate, rate)
    volume <- rep(2/3, 3)
    wait(play(mkWave(
        c("C3", "G2", "C2"),
        c("C4", "G3", "C3"),
        c("E4", "D4", "E4"),
        c("G4", "F4", "G4"),
        c("C5", "G4", "C5"),
        c("", "B3", ""),
        length = length, volume = volume, rate = rate), rate))
}

fanfare <- function(rate = 44100){
    scale <- c("F4", "A4", "D5", "F5", "", "D5", "F5")
    length <- c(rep(1/5, 6), 1/2)*rate
    volume <- rep(1/2, length(scale))
    wait(play(mkWave(scale, length = length, volume = volume, rate = rate), rate))
}

bigben <- function(rate = 44100){
    scale <- c("F4", "A4", "G4", "C4", "F4", "G4", "A4", "F4", "A4", "F4", "G4", "C4", "C4", "G4", "A4", "F4")
    length <- rep(c(1, 1, 1, 2), 4) * rate
    volume <- rep(1/2, length(scale))
    wait(play(mkWave(scale, length = length, volume = volume, rate = rate), rate))
}

bow()
fanfare()
correct()
wrong()
bigben()