【R/Stan】項目反応理論でM-1を分析してみる
RとStanを使って、過去三回のM-1についての分析を項目反応理論を応用して分析してみました。目次は以下の通りです。
はじめに
こんにちは。最近Stan力上げないとヤバいなと思いネットサーフィンをしていたら、Stanを用いてM-1について分析する記事をいくつか発見しました（といっても2個目は去年のゼミで話題になりました）。
めちゃくちゃ面白くないですか？Stanで柔軟なモデリングをすれば何か漫才についての色々なことが見えてきそうな気がします。
そこで、私自身もStanの練習を兼ねてM-1について分析してみようと思いました。
何を分析するか
この記事では、過去三回のM-1の成績データを使って、
をStanでモデルを組んで推定しようと思います。
では、なぜ最近3年のデータしか使わないといえば、最近の動向が気になるのはもちろんのことですが、この3年では審査員が固定されているので分析が楽ちんという点があります。この点については後ほど説明します。
どうやって分析するのか
潜在的なものを推定したいのでベイズを使うのですが、ここでは項目反応理論（Item Response Theory, IRT）という手法を応用しようと思います。私自身は政治学の畑の人ですが、卒論で心理学や教育学系の手法であるIRTを使うので、これは格好の機会だなと思い、やってみることにします（実際のところ、政治学でのIRT応用例は結構あります）。
項目反応理論を応用する
IRTについては、いろんなところで解説されつくされていると思うのでWikipediaでもいいので適当に見てください。
私は主にこの本で勉強しています。
通常のIRTのモデルは、以下のようなロジスティック回帰の時に用いられる形で書かれます（またはプロビット）。
$$
\text{Pr}(Y_i=1) = \frac{\exp(\beta_j\theta_i - \alpha_j)}{1 + \exp(\beta_j \theta_i - \alpha_j)}
$$
一般に、$i$をあるテストの被験者とし、$j$をあるテストの項目（問題）とします。$\theta_i$は個人$i$の潜在的な能力とされ、$\alpha_j$*1は項目$j$の困難度、$\beta_j$は識別パラメータと呼ばれ、能力$\theta_i$が高い人とそうでない人とを分けている度合いとして捉えられます（要は曲線の傾き）。
しかし、今回扱うM-1のデータは従属変数$Y_i$がバイナリーでもなければ多値選択変数でもなく、得点というカウント変数です。そこで、カウントといえばポアソン分布の登場です。「上記のロジットをポアソンへと拡張すれば、カウント変数でもIRTできるんじゃね？」と考えました。実際に、Slapin & Proksch (2008) が提案したWordfishというモデルも同じような発想だと思います。つまり、
\begin{align}
Y_i &\sim Poisson(\lambda_i),\\
\lambda_i &= \exp[\beta_j(\theta_i - \alpha_j)].
\end{align}
というモデルを考えます。$\beta_j(\theta_i - \alpha_j)$をlinear aggregator内に投入した形です。この発想が正しいかどうかは微妙ですが、一応テストしてみました。
上が$\alpha$の値を動かしたときのもので、下が$\beta$の値を動かしたものです。通常の場合と同じような曲線の動きが見て取れると思います。なので、不安ですが概ね大丈夫かと。。。
それでは、モデルの説明をしたいと思います。まず、ここではM-1の分析をしたいので項目$j$を審査員$j$であるとして、困難度である$\alpha_j$をその審査員$j$の採点に関する厳しめ度合いということにします。また、M-1の文脈で$\beta$をどう表そうかと考えてましたが、良い日本語が浮かばなかったので、良い出演者$i$とそうでないものを分けるパラメータとして考えます。つまり、良い出演者にはしっかりといい点数が付き、逆にそうでない出演者には厳しい値が付くというような識別度合いですね。また、$\theta_i$を出演者$i$の潜在的な漫才力を表すとします。
さらに今回推定するデータは、時間が3年間に渡っているため、時間の影響考慮をするためにモデルに時間のパラメータ$\gamma_{t}$を加えることにします。
\begin{align}
Y_i &\sim Poisson(\lambda_i),\\
\lambda_i &= \exp[\beta_j (\theta_i - \alpha_j) + \gamma_{t}].
\end{align}
とは言っても、時間に関するパラメータは実際のIRTにはほとんど入ることはないと思います。
それでは、早速分析していきましょう。
パッケージの読み込み
以下のパッケージを読み込みます。{tidybayes}はstan回した後の事後処理に非常に便利です。
library(rstan) library(tidyverse) library(tidylog) library(tidybayes)
また、以下のStanに関する設定をします。これやると速くなる（？）。
options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE)
データの読み込み
今回用いるデータは、過去三回のM-1（2018 ~ 2020）の成績データです。M-1のサイトやら何やらを参照して作成しました（ミスありましたらスミマセン）。ダウンロードはこちら。
# そのまま読み込むと文字化けするのでencodingをセットしてます。 df <- read.csv("data/m1_data.csv", encoding = "shift_jis") %>% as_tibble() # tibbleが好きなので
中身を確認しましょう。
dim(df) head(df) tail(df)
[1] 161 4
# A tibble: 6 x 4 time group referee score <int> <chr> <chr> <int> 1 14 ギャロップ オール巨人 87 2 14 ギャロップ 松本人志 86 3 14 ギャロップ 上沼恵美子 89 4 14 ギャロップ 塙宜之 89 5 14 ギャロップ 富澤たけし 87 6 14 ギャロップ 立川志らく 86
# A tibble: 6 x 4 time group referee score <int> <chr> <chr> <int> 1 16 東京ホテイソン 松本人志 86 2 16 東京ホテイソン 上沼恵美子 92 3 16 東京ホテイソン 塙宜之 85 4 16 東京ホテイソン 富澤たけし 91 5 16 東京ホテイソン 立川志らく 89 6 16 東京ホテイソン 礼二 88
ここでは、2回連続で出場している組も何組かあり（e.g. 和牛）、分析の際に重複すると面倒くさいので、その組については一番新しい出場データのみを適用しています。今回やりたいことは、審査員$j$をIRTでいうところのItemとして捉えることで各出演者$i$の潜在的な能力$\theta_i$（潜在的な漫才力）をIRTで推定することです。この時、時系列の分析にしてしまうと欠損値の扱いが大変になってしまったりするし、単純にデータを集めてしまうと重複する出演者が出てきてしまうので、それを防ぐために上記の収集ルールにしています。
それでは、このデータを用いてまずは簡単な可視化をします。まず、出演者は全部で何組いるのか確認します。
length(unique(df$group))
[1] 23
通常、M-1決勝は10組ずつなので23組はちょっとおかしいですが、上記のルールを適用したからこうなっているだけなので気にしないようにします。また、審査員も確認します。
length(unique(df$referee))
[1] 7
先ほども説明した通り、過去三回は審査員が全員同じなので、データとしては非常に分析しやすいです。
出演者ごとの総得点を見てみます。ここで、いきなりボンッと得点に関するグラフを貼ってしまうと、どの組も600点以上を取っているのであまり得点の差が見えないようなグラフができてしまいます。そこで、全出演者の総得点から最低総得点を引いて+1した値を採用します。
df %>% group_by(group) %>% summarise(sum = sum(score)) %>% mutate(sum = sum - min(sum) + 1) %>% ggplot(aes(x = fct_reorder(group, sum), y = sum)) + # fct_reorder(factor, value)で大きい順に並び変える geom_bar(stat = "identity", fill = "#00b3b3") + coord_flip() + ylab("Total Score (Adjusted)") + xlab("") + theme_light()
14 ~ 15回の最終決戦チームが結構高得点を取ってますが、前回は何か微妙な感じですね（笑）。次に、審査員ごとの総得点を見てみます。
df %>% group_by(referee) %>% summarise(sum = sum(score)) %>% mutate(sum = sum - min(sum) + 1) %>% ggplot(aes(x = fct_reorder(referee, sum), y = sum)) + geom_bar(stat = "identity", fill = "#00b3b3") + coord_flip() + ylab("Total Score (Adjusted)") + xlab("") + theme_light()
ここから単純に見れば、巨人さんは厳しめに得点を付けていて上沼さんは甘めに付けているということが何となくですがわかります。しかし、本当にそうなのかはわかりません。
項目反応理論（IRT）をやってみる
モデル
それではIRTの実装に入りたいと思います。先述の通り、ここで行うIRTは従属変数がバイナリーというわけではなく、また多値の選択変数でもないただのカウント変数ため、いわゆるロジットやプロビットで行われる通常の2Pモデルではなく、ポアソンへと応用したモデルを扱います。generative modelは、
\begin{align}
Y_i &\sim Poisson(\lambda_i),\\
\lambda_i &= \exp[\beta_j (\theta_i - \alpha_j) + \gamma_{t}].
\end{align}
です。
データの整形
IRTを行う際には、IRTのための行列（以下、IRT行列）を作るのが慣例というか恒例というか義務（？）なので、早速作りたいと思います。IRT行列とは、行に回答者$i$・列に項目$j$、要素は正解or不正解のような形です。ここでの事例に合わせるならば、行は出演者$i$・列は審査員$j$、要素は得点です。以下で加工してみます。
df %>% select(group, referee, score) %>% mutate(score = score - min(df$score)) %>% pivot_wider(names_from = referee, values_from = score) %>% # wider型へ変換 as.matrix() -> m_irt
得点に関してですが、90点とかだと大きすぎるので、各得点から最低点を引いた値を採用しています。諸々残りをやります。
gname <- m_irt[, 1] # 1列目（出演者名）を抽出 m_irt <- m_irt[, -1] # 1列目を削除 rownames(m_irt) <- gname # 行ラベルとして抽出した出演者名を適用 rname <- colnames(m_irt) # 後ほど使うので審査員名も抽出
できた行列を見てみます。
m_irt
| オール巨人 | 松本人志 | 上沼恵美子 | 塙宜之 | 富澤たけし | 立川志らく | 礼二 | |
|---|---|---|---|---|---|---|---|
| ギャロップ | 7 | 6 | 9 | 9 | 7 | 6 | 10 |
| ジャルジャル | 13 | 12 | 8 | 13 | 10 | 19 | 13 |
| スーパーマラドーナ | 7 | 5 | 9 | 9 | 9 | 8 | 10 |
| トムブラウン | 7 | 11 | 6 | 13 | 9 | 17 | 10 |
| ミキ | 10 | 8 | 18 | 10 | 10 | 9 | 13 |
| ゆにばーす | 4 | 0 | 4 | 2 | 6 | 7 | 11 |
| 霜降り明星 | 13 | 14 | 17 | 18 | 11 | 13 | 16 |
| かまいたち | 13 | 15 | 15 | 15 | 13 | 15 | 14 |
| からし蓮根 | 13 | 10 | 14 | 10 | 10 | 9 | 13 |
| すゑひろがりず | 12 | 9 | 12 | 11 | 10 | 12 | 11 |
| ぺこぱ | 13 | 14 | 16 | 14 | 14 | 11 | 12 |
| ミルクボーイ | 17 | 17 | 18 | 19 | 17 | 17 | 16 |
| 和牛 | 12 | 12 | 12 | 16 | 11 | 16 | 13 |
| アキナ | 9 | 5 | 12 | 7 | 8 | 10 | 11 |
| インディアンス | 9 | 10 | 13 | 5 | 9 | 9 | 10 |
| ウェストランド | 8 | 10 | 12 | 5 | 11 | 6 | 10 |
| おいでやすこが | 12 | 15 | 14 | 13 | 13 | 16 | 15 |
| オズワルド | 8 | 8 | 12 | 15 | 11 | 13 | 15 |
| ニューヨーク | 8 | 12 | 14 | 13 | 13 | 11 | 11 |
| マヂカルラブリー | 8 | 13 | 14 | 14 | 14 | 10 | 16 |
| 錦鯉 | 7 | 9 | 13 | 15 | 12 | 15 | 12 |
| 見取り図 | 11 | 11 | 15 | 13 | 12 | 13 | 13 |
| 東京ホテイソン | 6 | 6 | 12 | 5 | 11 | 9 | 8 |
出演者の順番はあまり気にしないでください。
次に、Stanへ渡すためのデータを作らなくてはなりません。Stanくんはlist型のデータではないと読み込んでくれないのがちょっとダルい点です。
Y <- as.numeric(m_irt) # 得点 N <- length(Y) # observation i <- rep(1:dim(m_irt)[1], times = dim(m_irt)[2]) # 出演者の観察数 j <- rep(1:dim(m_irt)[2], each = dim(m_irt)[1]) #審査員の観察数 I <- max(i) # 出演者数 J <- max(j) # 審査員数 # 時間の影響考慮もしたいので抜き出します df %>% arrange(referee) -> df2 t <- df2$time - min(df2$time) + 1 # 時間の観察数 T <- max(t) # 時間数 data <- list(I = I, J = J, N = N, T = T, i = i, j = j, t = t, Y = Y)
Stanの実装
次に、Stanを書きます。今回のStanコードは以下の通りです（pois_irt.stanという名前で保存しました）。
また、priorとして
$$
\begin{align}
\alpha &\sim \mathcal{N}(0, 1), \\\
\beta &\sim \mathcal{N}(0, 10), \\\
\gamma &\sim \mathcal{N}(0, 10), \\\
\theta &\sim \mathcal{N}(0, 1).
\end{align}
$$
を与えます。
data{
int<lower=1> N; // # of observations
int<lower=1> I; // # of groups
int<lower=1> J; // # of referees
int<lower=1> T; // # of time
int<lower=1, upper=I> i[N]; // groups observation
int<lower=1, upper=J> j[N]; // referees observation
int<lower=1, upper=T> t[N]; // time observation
int<lower=0> Y[N]; // score
}
parameters {
vector[T] gamma; //time parameter
vector[J] alpha; //difficulty
vector[J] beta; //discrimination
vector[I] theta; //latent ability
}
model {
for (n in 1:N){
Y[n] ~ poisson(exp(beta[j[n]] * (theta[i[n]] - alpha[j[n]]) + gamma[t[n]]));
}
theta ~ std_normal(); // for standardizing
alpha ~ std_normal();
beta ~ normal(0, 10);
gamma ~ normal(0, 10);
}
できたので早速Stanを回しましょう。最初にモデルを読み込みます。
model <- stan_model("stan/pois_irt.stan")
警告が出る時がありますが無視してかまいません。キック用のコードは以下の通りです。
fit <- sampling(model, data = data, iter = 20000, warmup = 5000, chains = 1, seed = 1, thin = 10, control = list(adapt_delta = 0.99, # controlなしでやったら警告が出たので入れてます max_treedepth = 15))
収束しているかチェックしましょう。
all(summary(fit)$summary[, "Rhat"] <= 1.10, na.rm = TRUE)
[1] TRUE
TRUEが返されているので、すべてのパラメータに関するRhatが1.1以下であったことがわかり、収束していると判断します。上記の便利な判定コードは、以下のslideshareを参照しました。めちゃくちゃ便利ですね。
Stanの便利な事後処理関数 from daiki hojo
www.slideshare.net
結果
次に、各パラメータの事後分布を可視化します。まずは、各出演者の潜在的な漫才力$\theta_i$を見てみましょう。この値が高ければ高いほど潜在的な漫才力が高いことを表します。また、推定値の線は細いものから順に95%、90%、80%の信用区間を示し、黒丸は事後分布の中央値です。
group <- tibble(group = gname, n = 1:length(gname)) # 出演者の名前を投入するためにデータフレームを作ります fit %>% spread_draws(theta[n]) %>% # stanの結果をtidyに変換 median_qi(.width = c(.95, .9, .8)) %>% # 信用区間やら諸々 left_join(group, by = "n") %>% ggplot(aes(x = theta, y = fct_reorder(group, theta), xmin = .lower, xmax = .upper, size = -.width)) + geom_pointinterval(interval_size_range = c(0.4, 1.5)) + theme_light() + xlab("Estimated Latent Score") + ylab("")
霜降り明星とミルクボーイ強いですねぇ～。得点だけで見たらミルクボーイが圧倒的でしたが、いざ分析してみると霜降りの方が強いとは。個人的にはトムブラウンが割と上位の方にいて嬉しいですね～。でんじゃらすじーさんを見ている感じがしてとても良い。
個人$i$だけでなく項目$j$の推定も見てみましょう。まずは、審査員の厳しめ度合いを表す$\alpha_j$から見てみます。
referee <- tibble(referee = rname, n = 1:length(rname)) # 審査員の名前のデータフレーム fit %>% spread_draws(alpha[n]) %>% median_qi(.width = c(.95, .9, .8)) %>% left_join(referee, by = "n") %>% ggplot(aes(x = alpha, y = fct_reorder(referee, alpha), xmin = .lower, xmax = .upper, size = -.width)) + geom_pointinterval(interval_size_range = c(0.4, 1.5)) + theme_light() + xlab("Estimated Difficulties") + ylab("")
結果としては言われている通り巨人さんでが厳しめで、上沼さんは甘めなのでしょうか。しかし、信用区間がかなり広く不確実性が大きいので何とも言えないです。
最後に$\beta$の値を見てみます。解釈が難しいパラメータですが、、、、
fit %>% spread_draws(beta[n]) %>% median_qi(.width = c(.95, .9, .8)) %>% left_join(referee, by = "n") %>% ggplot(aes(x = beta, y = fct_reorder(referee, beta), xmin = .lower, xmax = .upper, size = -.width)) + geom_pointinterval(interval_size_range = c(0.4, 1.5)) + theme_light() + xlab("Estimated Saliency") + ylab("")
これ、何と呼んでいいかわからず横軸をSaliency（政治学では$\beta$はsailiency parameterと呼ばれることが多い）にしてしまいました。塙さんが一番大きな値を獲得してますね。能力のある芸人にはしっかりと点数を付けて、そうでない人には厳しいといった評価を与えているという風に解釈できますね。そう考えると、礼二さんが一番可もなく不可もなくのような採点をしているのでしょうか。
初めてお笑いの分析をしてみましたが結構面白いですね。今回やったことがいろんな方面に活かせたら結構面白いですね～。
【R】Bayesian heteroskedastic linear regresssionをRで実装してみる（実装編）
はじめに
前回は以下の記事にて、Bayesian heteroskedastic linear regressionの理論について整理しました。
sunaninattahito.hatenablog.com
この記事では、前回の記事で導出したposteriorを基に、Rで自力でMCMCを実装してみることにします。また、この記事はGreenberg (2011)の教科書の練習問題がベースです。
改善点や間違い等あればぜひ指摘してください。
パッケージの読み込み
まず最初に以下のパッケージを読み込みましょう。{tidylog}は{tidyverse}で行ったデータ加工の結果を表示してくれる便利な代物ですが、必要がなければ読み込まなくても全く問題ありません。
library(tidyverse) library(tidylog) library(mvtnorm) #多変量正規分布を発生させるやつ set.seed(1)
データの読み込み
データは、Greenberg (2011) 8.4 Excersisesの問題8.3に基づいて著者のサイト(https://www.stat.berkeley.edu/users/statlabs/labs.html) にある``Birth weight II"というデータを使います。このデータは生まれてきた赤ちゃんの体重とそれに対して影響を与えるであろう変数が盛り込まれています。早速、読み込みましょう。
df <- read.table("https://www.stat.berkeley.edu/users/statlabs/data/babies.data", header = TRUE) %>% as_tibble() # tibbleが好きなので
リンク先が.dataファイルなのでread.table関数を使用しました。少しデータを確認してみましょう。
dim(df) head(df) tail(df)
[1] 1236 7
# A tibble: 6 x 3 variable value iterations <chr> <dbl> <int> 1 intercept -103. 1 2 intercept -76.8 2 3 intercept -83.7 3 4 intercept -63.4 4 5 intercept -115. 5 6 intercept -68.5 6
# A tibble: 6 x 3 variable value iterations <chr> <dbl> <int> 1 smoke -8.07 995 2 smoke -7.58 996 3 smoke -8.25 997 4 smoke -9.25 998 5 smoke -7.94 999 6 smoke -7.05 1000
また、記録されている変数は
names(df)
[1] "bwt" "gestation" "parity" "age" "height" "weight" [7] "smoke"
- bwt: 赤ちゃんの出生時の体重 (ounces)
- gestation: 妊娠期間
- parity: 第一子かどうか (0 = first born)
- age: 出産時の年齢
- height: 母親の身長（inch）
- weight: 母親の妊娠時の体重 (pounds)
- smoke: 母親が喫煙するかどうか (0 = not now, 1 = yes now)
です。次に、データから欠損値がある観察iを削除していいとの事のなので、消しちゃいましょう。ちなみに、gestation・weightでは999、age・heightでは99、smokeでは9が欠損値とされています。
df %>% filter(gestation != 999 & age != 99 & height != 99 & weight != 999 & smoke != 9) -> df
filter: removed 62 rows (5%), 1,174 rows remaining
上記のように{tidylog}は、filterなどの結果を教えてくれます。
記述統計をして、データがどのようになっているか確認します。
summary(df)
bwt gestation parity age Min. : 55.0 Min. :148.0 Min. :0.0000 Min. :15.00 1st Qu.:108.0 1st Qu.:272.0 1st Qu.:0.0000 1st Qu.:23.00 Median :120.0 Median :280.0 Median :0.0000 Median :26.00 Mean :119.5 Mean :279.1 Mean :0.2624 Mean :27.23 3rd Qu.:131.0 3rd Qu.:288.0 3rd Qu.:1.0000 3rd Qu.:31.00 Max. :176.0 Max. :353.0 Max. :1.0000 Max. :45.00 height weight smoke Min. :53.00 Min. : 87.0 Min. :0.000 1st Qu.:62.00 1st Qu.:114.2 1st Qu.:0.000 Median :64.00 Median :125.0 Median :0.000 Mean :64.05 Mean :128.5 Mean :0.391 3rd Qu.:66.00 3rd Qu.:139.0 3rd Qu.:1.000 Max. :72.00 Max. :250.0 Max. :1.000
実装
データも準備できたので、早速実装してみようと思います。また、簡単のために予めデータをattachしておきます。
attach(df)
この問題の要旨としては、赤ちゃんの出生時体重であるbwtを従属変数として、その他の変数を独立変数とする形になります。要は赤ちゃんの出生時体重の決定要因を探ろうみたいな感じですね。
まず初めに、分析で用いるデータを行列に変換します。
# 行列X (N × K)を作成 X <- cbind(1, gestation, parity, age, height, weight, smoke) # Yを N × 1 行列に変換 Y <- as.matrix(bwt) # 後ほど使用するので変数名を保存 X_name <- colnames(X) # 先頭をinterceptという名前にする X_name[1] <- "intercept"
Posteriors
次に各パラメータのposteriorを計算するための関数を以下で定義していきます。
βのposterior
ここで$\beta$のposteriorは、
where
$$
\begin{align}
B_1 &= [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1},\\
\tilde{\beta} &= B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y].
\end{align}
$$
でした。これに基づいて、$\beta$のposteriorを計算する関数beta_postを以下で定義します。
beta_post <- function(Y, X, b0, B0_inv, sigma2, Lambda){ B1 <- solve(B0_inv + t(X) %*% Lambda %*% X / sigma2) beta_mean <- B1 %*% (B0_inv %*% b0 + t(X) %*% Lambda %*% Y / sigma2) beta <- matrix(rmvnorm(n = 1, mean = beta_mean, sigma = B1), , 1) return(beta) }
σ2のposterior
次に、σ2のposteriorを導出する関数を定義します。posteriorは、
$$
σ^2|\beta, \lambda, Y \sim \mathcal{IG}\left(\frac{\alpha_1}{2}, \frac{\delta_1}{2}\right),
$$
where
$$
\begin{align}
\alpha_1 &= \alpha_0 + n, \\\
\delta_1 &= \delta_0 + (Y-X\beta)^{T} \Lambda (Y-X\beta).
\end{align}
$$
です。よって、以下のsigma2_postを定義します。
sigma2_post <- function(Y, X, alpha0, delta0, beta, Lambda){ alpha1 <- alpha0 + nrow(X) delta1 <- delta0 + t(Y - X %*% beta) %*% Lambda %*%(Y - X %*% beta) sigma2 <- MCMCpack::rinvgamma(n = 1, alpha1 / 2, delta1 / 2) return(sigma2) }
λiのposterior
最後に$\lambda_i$のposteriorを見ていきます。posteriorは以下のように表せます。
$$
\begin{align}
\lambda_i | \beta, σ^{2}, Y \sim Ga\left(\frac{\nu_1}{2}, \frac{\nu_{2, i}}{2}\right),\\
\end{align}
$$
where
$$
\begin{align}
\nu_1 &= \nu + 1,\\
\nu_{2, i} &= \nu + σ^{-2}(Y_i - X_i\beta).
\end{align}
$$
また、ここで自由度$\nu$は既知として、関数を以下で定義します。
lambda_i_post <- function(Y, X, b0, B0_inv, nu, sigma2, Lambda) { nu1 <- nu + 1 nu2_i <- matrix(NA, nrow(X), 1) lambda_i <- matrix(NA, nrow(X), 1) b <- beta_post(X = X, Y = Y, b0 = b0, B0_inv = B0_inv, sigma2 = sigma2, Lambda = Lambda) nu2_i <- nu + ((Y - X %*% b) / sigma2) for (i in 1:nrow(nu2_i)) { lambda_i[i] <- rgamma(n = 1, scale = (nu1 / 2), shape = (nu2_i[i] / 2)) } return(lambda_i) }
Gibbs Sampler
次に、Gibbs samplingをするために関数を定義します。
gibbs <- function(iter, warmup, thin, verbose = FALSE, Y, X, b0, B0_inv, alpha0, delta0, nu){ set.seed(1) # 時間計測スタート stime <- Sys.time() # storage store_beta <- matrix(NA, iter / thin, ncol(X)) ### storage for beta store_sigma2 <- matrix(NA, iter / thin, 1) ### storage for sigma2 store_lambda <- matrix(NA, iter / thin, nrow(X)) ### storage for lambda # initital value sigma2 <- 1 # initial value fot sigma2 Lambda <- diag(1, nrow(X)) # initital value for Lambda total_iter <- iter + warmup for (i in 1:total_iter){ beta <- beta_post(Y, X, b0, B0_inv, sigma2, Lambda) sigma2 <- sigma2_post(Y, X, alpha0, delta0, beta, Lambda) lambda_i <- lambda_i_post(Y, X, b0, B0_inv, nu, sigma2, Lambda) Lambda <- diag(as.numeric(lambda_i)) # save posteriors... if (i > warmup & ((i - warmup) %% thin == 0)){ store_beta[(i - warmup) / thin, ] <- beta # save beta store_sigma2[(i - warmup) / thin, ] <- sigma2 # save sigma2 store_lambda[(i - warmup) / thin, ] <- t(lambda_i) # save lambda } if (verbose == TRUE & i <= warmup & i %% 100 == 0){ cat("Warmup", i,"DONE.", crayon::yellow("Total time:", Sys.time() - stime), "\n") } if (verbose == TRUE & i > warmup & i %% 100 == 0){ cat("Iterations", i - warmup,"DONE.", crayon::yellow("Total time:", Sys.time() - stime), "\n") } } # naming colnames(store_beta) <- X_name colnames(store_sigma2) <- "sigma2" colnames(store_lambda) <- paste0("observation", 1:nrow(X)) # Convert into data frame store_beta <- as_tibble(store_beta) %>% mutate(iterations = 1:(iter / thin), .before = X_name[1]) store_sigma2 <- as_tibble(store_sigma2) %>% mutate(iterations = 1:(iter / thin), .before = "sigma2") store_lambda <- as_tibble(store_lambda) %>% mutate(iterations = 1:(iter / thin), .before = "observation1") return(list(beta_posterior = store_beta, sigma2_posterior = store_sigma2, lambda_posterior = store_lambda))
Greenberg (2011)の問題8.3に基づいて自由度$\nu$ = 5とします。また、このモデルをspecifyするために$\beta_0$, $B_0^{-1}$, $\alpha_0$, $\delta_0$にpriorを置きます。今回は適当にnoninformative priorを置くこととします。ここでは、
$$
\begin{align}
\beta_0 &= 0, \\
B_0^{-1} &= 0, \\
\alpha_0 &= 0.001, \\
\delta_0 &= 0.001.
\end{align}
$$
とします。それでは、早速MCMCしていきましょう！
fit <- gibbs(Y = Y, X = X, iter = 3000, warmup = 1000, thin = 3, verbose = TRUE, b0 = matrix(0, ncol(X), 1), B0_inv = diag(0, ncol(X), ncol(X)), alpha0 = 0.001, delta0 = 0.001, nu = 5)
Warmup 100 DONE. Total time: 9.5207188129425 Warmup 200 DONE. Total time: 16.6507158279419 Warmup 300 DONE. Total time: 22.6647148132324 Warmup 400 DONE. Total time: 28.8763778209686 Warmup 500 DONE. Total time: 34.9007859230042 Warmup 600 DONE. Total time: 41.0568199157715 Warmup 700 DONE. Total time: 46.8590228557587 Warmup 800 DONE. Total time: 52.8343329429626 Warmup 900 DONE. Total time: 58.8880109786987 Warmup 1000 DONE. Total time: 1.08210526704788 Iterations 100 DONE. Total time: 1.17846973339717 Iterations 200 DONE. Total time: 1.27412329912186 Iterations 300 DONE. Total time: 1.37088464895884 Iterations 400 DONE. Total time: 1.46519904931386 Iterations 500 DONE. Total time: 1.56177023251851 Iterations 600 DONE. Total time: 1.65910194714864 Iterations 700 DONE. Total time: 1.75826054811478 Iterations 800 DONE. Total time: 1.85463916460673 Iterations 900 DONE. Total time: 1.94738284746806 Iterations 1000 DONE. Total time: 2.04625216325124 Iterations 1100 DONE. Total time: 2.14010365009308 Iterations 1200 DONE. Total time: 2.23587293227514 Iterations 1300 DONE. Total time: 2.33102181355158 Iterations 1400 DONE. Total time: 2.42716313203176 Iterations 1500 DONE. Total time: 2.52142798105876 Iterations 1600 DONE. Total time: 2.61657724777857 Iterations 1700 DONE. Total time: 2.71183398167292 Iterations 1800 DONE. Total time: 2.80728929837545 Iterations 1900 DONE. Total time: 2.90187638203303 Iterations 2000 DONE. Total time: 2.99862943092982 Iterations 2100 DONE. Total time: 3.09441526333491 Iterations 2200 DONE. Total time: 3.18699366648992 Iterations 2300 DONE. Total time: 3.28234141667684 Iterations 2400 DONE. Total time: 3.37668879826864 Iterations 2500 DONE. Total time: 3.47107758124669 Iterations 2600 DONE. Total time: 3.56890536546707 Iterations 2700 DONE. Total time: 3.66383141676585 Iterations 2800 DONE. Total time: 3.7627916653951 Iterations 2900 DONE. Total time: 3.85598178307215 Iterations 3000 DONE. Total time: 3.95359458128611
mutate: new variable 'iterations' (integer) with 1,000 unique values and 0% NA mutate: new variable 'iterations' (integer) with 1,000 unique values and 0% NA mutate: new variable 'iterations' (integer) with 1,000 unique values and 0% NA
結構遅くてイライラしますが気長に待ちましょう。
結果
結果を確認してみましょう。まずはconvergenceをチェックします。$\lambda_i$については各$i$ごとにチェックしなければならず、結構面倒くさいので割愛します。
βのconvergence
{ggplot2}に含まれるfacet_wrap関数を活用して、変数ごとにconvergenceを見ていきたいと思います。まずは、プロット用に脳筋的にデータを加工していきます。
# βのposteriorを抽出 fit_df <- fit$beta_posterior # 各変数の名前を抽出 v_name <- colnames(fit_df[, -1]) # 各変数の名前を1000回ずつ入れた縦ベクトルを作る v <- rep(v_name, each = nrow(fit_df)) %>% as_tibble() # 各変数の推定値を縦にどんどん追加していく for (i in 2:8) { if (exists("value") == FALSE) { value <- fit_df[, i] %>% as.matrix() } else{ value2 <- fit_df[, i] %>% as.matrix() value <- rbind(value, value2) } } # tibble化 value %>% as_tibble() -> value # 各変数の名前が入った縦ベクトルにvalue追加する v %>% bind_cols(value) %>% rename(variable = value, value = intercept) %>% mutate(iterations = rep(1:1000, times = 7)) -> plot_df # 確認 head(plot_df) tail(plot_df)
rename: renamed one variable (variable) mutate: new variable 'iterations' (integer) with 1,000 unique values and 0% NA
# A tibble: 6 x 3 variable value iterations <chr> <dbl> <int> 1 intercept -103. 1 2 intercept -76.8 2 3 intercept -83.7 3 4 intercept -63.4 4 5 intercept -115. 5 6 intercept -68.5 6
# A tibble: 6 x 3 variable value iterations <chr> <dbl> <int> 1 smoke -8.07 995 2 smoke -7.58 996 3 smoke -8.25 997 4 smoke -9.25 998 5 smoke -7.94 999 6 smoke -7.05 1000
早速、プロットしてみましょう。
plot_df %>% ggplot(aes(x = iterations, y = value)) + geom_line() + facet_wrap(~ variable, scales = "free") + theme_light() + xlab("Iterations") + ylab("Estimated") + ggtitle("Convergence check for β")
大体収束してそうな気がします。
σ2のconvergence
次に、σ2についてのconvergenceを見てみましょう。
fit$sigma2_posterior %>% ggplot(aes(x = iterations, y = sigma2)) + geom_line() + theme_light() + xlab("Iterations") + ylab("Estimated") + ggtitle(expression(paste("Convergence check for ", sigma^2)))
こちらもよさそうです。
Posterior distributions
posteriorsも確認しましょう。先ほど作ったplot_dfというデータフレームを用いて作図します。また、各変数ごとにposterior
meanを図示するために、数ステップ踏みます。
# 変数ごとにposterior meanを計算 plot_df %>% group_by(variable) %>% summarise(mean = mean(value)) -> mean_df
group_by: one grouping variable (variable) summarise: now 7 rows and 2 columns, ungrouped
plot_df %>% ggplot(aes(x = value)) + geom_density(alpha = .3, fill = "black") + facet_wrap(~ variable, scales = "free") + geom_vline(data = mean_df, aes(xintercept = mean), linetype = "dashed") + theme_light() + xlab("Estimated") + ylab("Density") + ggtitle("Posterior distributions with posterior mean")
density中央付近にある破線はposterior meanを表しています。interceptの値がデカすぎますが、OLSでもほとんど同じ位の推定値だったので見逃します。
また、通常のBayesian linear regressionの結果との比較をするために、以下のようなプロットを作ってみました。
ここから見ると、Heteroskedasticityを考慮したモデルの方がposteriorの裾がやや長いことが確認できます。上手く実装できているか不安ですが（（（（（（
【R】Bayesian heteroskedastic linear regresssionをRで実装してみる（理論編）
はじめに
先日、Twitterを見ていたら、「StanとかJAGSは、MCMCを自力で書けない人が使うべきでない」というような指摘をするツイートを見つけました。 それを見て、私自身「ギクッーー」となりました。この指摘は本当におっしゃる通りであり、Stanなんかは事前分布と尤度を書くだけで勝手にベイズ推定ができちゃいます。その他背後でどのような計算が行われているかは、やはり自分の手で書いてみないと何もわかりません。何が起こっているかわからないものを使うのは個人的にちょっと気持ち悪い感じがします。何でもかんでもStanにぶっこんでベイズベイズ～とか言ってましたが、確かに自分の手でMCMCを実装するべきだなと痛感しました。
そこで、この記事では前学期に受講したベイズの授業で扱った教科書である``Introduction to Bayesian econometrics" (Greenberg, 2011) の練習問題を基にして、自力でMCMCをしてみよう！という記事になります。
ここでは何を扱うかといえば、通常のlinear regressionだけではあまり面白くないので、``Heteroskedasticity"を考慮したHeteroskedastic linear regressionをRで実装します。ちなみに、Greenberg (2011)の 8.4 Excersisesから8.1と8.3の練習問題に基づいています。今回は理論編です。実装編は以下から飛んでください。また、間違いや改善点等ございましたら是非指摘してください。
sunaninattahito.hatenablog.com
Heteroskedasticityについて
まず、そもそもHeteroskedasticityとは何かについて説明していきます。通常のlinear regressionにおいては、誤差項の分散である
がすべての
について一定であることを仮定します。しかし、実際に扱う（現実の）データは、すべての観察
において
が一定であるとは限りません*1。
ここで、出てくるのがHeteroskedasticityの問題です。もし、分散が均一でないのなら、通常のlinear regressionによる推定はバイアスを持つことになります。そこで、このような問題に対処しようとするモデルが、今回実装しようとするHeteroskedastisic linear regressionです。
通常のlinear regressionは
といったgenerative modelですが、heteroskedastic linear regressionでは
というような形で、表すことになります。ここでは、分散を
ごとのパラメータである
で補正するような形になってます。この
の導入によって、いわゆるheteroskedasticityの問題に対処します。ここで、
とは自由度（degree of freedom）を表し、また簡単のために自由度
は既知とします。
Posteriorの導出
次に、それぞれのパラメーターについてのposteriorを導出していきたいと思います。このセクションはあくまで私の理解度確認テストのようなものなので、「早くRで実装したいンだが」という方は実装編に飛んじゃってください。
sunaninattahito.hatenablog.com
Setting
まず、以下のようなモデルを仮定します。
\begin{align*}
Y_i &= X_i\beta + \epsilon_i, \\
\text{where } &\epsilon_i \sim \mathcal{N}(0, σ^{2}).
\end{align*}
ここで、 は
かつ 1 列目が 切片である（=1）行列であり、
は
の係数の行列であるとします。
さらに、 を
となるような
行列とし、
を観察$i$の全数と対応するように
の行列とします。また、
を
となるような
行列とします。そうした場合、上式を以下のように書き換えることができます。
ここで、Heteroskedastic errors を得るために、以下を仮定します。
Likelihood
次に、Likelihood function (尤度関数) を見ていきます。ここで、$\Lambda$を$\lambda_i$についての対角行列であるとして
$$
\Lambda = \left(
\begin{array}{cccc}
\lambda_1& 0 & \ldots & 0 \\
0 & \lambda_2 & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & \lambda_n
\end{array}
\right)
$$
とします。また、$Y=(Y_1,\cdots, Y_n)$とすると、$\beta$, σ2, $\lambda_i$を所与とした場合の尤度関数は
$$
f(Y|\beta, σ^2, \lambda_i) = \prod\limits_{i=1}^{n}\left(\frac{\lambda_i}{2\pi σ^2}\right)^{1/2} \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right].
$$
となります。
Priors
次にpriorを見ていきます。
βについて
ここで、$\beta$に以下の多変量正規分布のpriorを置きます。
$$
\beta \sim \mathcal{N}_K(\beta_0, B_0).
$$
この時、$\beta$に関する事前分布の確率密度 $\pi (\beta)$は以下のように表すことができます。
$$
\pi(\beta) = \frac{1}{(2\pi)^{k/2}|B_0|^{1/2}} \exp\left[-\frac{1}{2}(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0)\right].
$$
σ2について
次に、σ2のpriorを見ていきます。ここでは、σ2に以下の逆ガンマ分布のpriorを置きます。
$$
σ^2 \sim \mathcal{IG}\left(\frac{\alpha_0}{2}, \frac{\delta_0}{2}\right).
$$
この時、σ2に関する事前分布の確率密度 $\pi (σ^{2})$は以下のように表すことができます。
$$
\pi(σ^2) =\frac{(\delta_0 / 2)^{\alpha_0/2}}{\Gamma(\alpha_0/2)} \left(\frac{1}{σ^2}\right)^{(\alpha_0/2) + 1} \exp\left[-\frac{\delta_0}{2σ^2}\right].
$$
λiについて
ここでは、$\lambda_i$に以下のガンマ分布のpriorを置きます。
$$
\lambda_i \sim Ga\left(\frac{\nu}{2}, \frac{\nu}{2}\right), \ \ \nu \text{ known}.
$$
ここで、自由度$\nu$は既知とします。この時、$\lambda_i$に関する事前分布の確率密度$\pi(\lambda_i)$は
$$
\pi(\lambda_i) = \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)}\lambda_i^{(v/2)-1} \exp\left[-\sum\limits_{i=1}^{n}\frac{\nu\lambda_i}{2}\right].
$$
となります。また、$\pi(\lambda) = \pi(\lambda_1, \cdots, \lambda_n)$と定めると
$$
\begin{align}
\pi(\lambda) &= \prod\limits_{i=1}^{n}\pi (\lambda_i)\\
&= \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)} \prod\limits_{i=1}^{n}\lambda_i^{(v/2)-1} \exp\left[-\sum\limits_{i=1}^{n}\frac{\nu\lambda_i}{2}\right].
\end{align}
$$
として表します。
Posterior
最後にPosteriorを導出していきます。これはめちゃくちゃダルいので心が折れそうになりますが頑張ります。
Joint posterior
Joint posterior はベイズの定理より以下のように書き表すことができます。
$$
\begin{align}
\pi(\beta, σ^2, \lambda|Y) &\propto \pi(\beta)\pi(σ^2)\pi(\lambda)f(Y|\beta, σ^2, \lambda)\\
&= \frac{1}{(2\pi)^{k/2}|B_0|^{1/2}} \exp\left[-\frac{1}{2}(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0)\right]\\
&\times \frac{(\delta_0 / 2)^{\alpha_0/2}}{\Gamma(\alpha_0/2)} \left(\frac{1}{σ^2}\right)^{(\alpha_0/2) + 1} \exp\left[-\frac{\delta_0}{2σ^2}\right]\\
&\times \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)} \prod\limits_{i=1}^{n}\lambda_i^{(v/2)-1} \exp\left[-\sum\limits_{i=1}^{n}\frac{\nu\lambda_i}{2}\right]\\
&\times \prod\limits_{i=1}^{n}\left(\frac{\lambda_i}{2\pi σ^2}\right)^{1/2} \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right]
\end{align}
$$
Conditional posterior of σ2
ここから先は、上記のjoint posteriorからそれぞれのconditional posteriorを求めることになります。まずは、導出のし易いσ2から見ていきます。σ2のconditonal posteriorを求めたいときは、上記のjoint posteriorからσ2に関する項のみに着目します。この時
$$
\begin{align}
\pi(σ^2|\beta, \lambda, Y) &\propto \left(\frac{1}{σ^2}\right)^{(\alpha_0/2) + 1} \exp\left[-\frac{\delta_0}{2σ^2}\right]\\
&\times \left(\frac{1}{σ^2}\right)^{n/2} \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right]\\
&= \left(\frac{1}{σ^2}\right)^{[(\alpha_0/ + n)/ 2] + 1}\exp\left[-\frac{1}{2σ^2}\{(\delta_0 + Y-X\beta)^{T}\Lambda(Y-X\beta)\}\right].
\end{align}
$$
となります。また、
$$
\begin{align}
\alpha_1 &= \alpha_0 + n, \\
\delta_1 &= \delta_0 + (Y-X\beta)^{T}\Lambda(Y-X\beta).
\end{align}
$$
とすると、上式は
$$
\begin{align}
\pi(σ^2|\beta, \lambda, Y) &\propto \left(\frac{1}{σ^2}\right)^{(\alpha_1/2) + 1} \exp\left[-\frac{\delta_1}{2σ^2}\right]
\end{align}
$$
と書き換えられます。この時、σ2のconditional posteriorは逆ガンマ分布からのサンプリングとなります。よって、
$$
\begin{align}
&σ^2 | \beta, \lambda, Y \sim \mathcal{IG} \left(\frac{\alpha_1}{2}, \frac{\delta_1}{2}\right)\\
\text{where }\ &\alpha_1 = \alpha_0 + n, \\
&\delta_1 = \delta_0 + (Y-X\beta)^{T}\Lambda(Y-X\beta).
\end{align}
$$
Conditional posterior of λi
Joint posteriorの式より、$\lambda_i$に関する式を抜き出します。$\lambda_i$のconditional posteriorは
$$
\begin{align}
\pi(\lambda_i|\beta, σ^2, Y_i) &\propto \lambda_i^{(v/2)-1} \exp\left[-\frac{\nu\lambda_i}{2}\right]\\
&\times \lambda_i^{1/2}\exp\left[-\frac{\lambda_i}{2σ^2}(Y_i - X_i\beta)^2\right]\ \ \cdots\ \ \text{Likelihoodの}Y, X\text{を}i\text{についての式に書き換えるとこうなる。} \\
&= \lambda_i^{[(\nu+1)/ 2]-1} \exp\left[-\frac{\lambda_i}{2}\{\nu + σ^{-2} (Y_i - X_i\beta\}\right].
\end{align}
$$
また、
$$
\begin{align}
\nu_1 &= \nu + 1\\
\nu_{2, i} &= \nu + σ^{-2}(Y_i - X_i\beta)
\end{align}
$$
とすると、上式は
$$
\pi(\lambda_i|\beta, σ^2, Y_i) \propto \lambda_i^{[(\nu_1)/ 2]-1} \exp\left[-\frac{\nu_{2, i}\lambda_i}{2}\right]
$$
と書き換えられます。この時、$\lambda_i$のconditional posteriorはガンマ分布からのサンプリングとなります。よって
$$
\begin{align}
&\lambda_i | \beta, σ^2, Y_i \sim Ga\left(\frac{\nu_1}{2}, \frac{\nu_{2, i}}{2}\right), \\
\text{where }\ \ &\nu_1 = \nu + 1,\\
&\nu_{2, i} = \nu + σ^{-2}(Y_i - X_i\beta).
\end{align}
$$
Conditional posterior of β
最後です。丁寧に解説するのはかなり面倒なので、ところどころ端折ってやります。Joint posteriorより、$\beta$のconditional posteriorは
$$
\begin{align}
\pi(\beta| \lambda, σ^2, Y) &\propto \exp\left[-\frac{1}{2}(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0)\right]\\
&\times \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right]\\
&= \exp \left[-\frac{1}{2}\{(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0) + σ^{-2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\}\right]
\end{align}
$$
となります。また、上式の$\exp$内を色々やって整理すると
$$
\begin{align}
\pi(\beta| \lambda, σ^2, Y) &\propto \exp\left[-\frac{1}{2}\{(\beta-\tilde{\beta})^{T}B_1^{-1}(\beta-\tilde{\beta})\}\right]\\
\text{where}\ \ B_1 &= [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1}, \\
\tilde{\beta} &= B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y].
\end{align}
$$
となる。この時、$\beta$のconditional posteriorは正規分布からのサンプリングとなります。よって
$$
\begin{align}
&\beta | \lambda, σ^2, Y \sim \mathcal{N}_K(\tilde{\beta}, B_1), \\
\text{where}\ \ &B_1 = [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1}, \\
&\tilde{\beta} = B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y].
\end{align}
$$
まとめ
ここまでで求められたposteriorを以下で並べてみましょう。
$$
\begin{align}
\beta|\lambda, σ^2, Y &\sim \mathcal{N}_K(\tilde{\beta}, B_1), \\
σ^2|\beta, \lambda, Y &\sim \mathcal{IG}\left(\frac{\alpha_1}{2}, \frac{\delta_1}{2}\right), \\
\lambda_i|\beta, σ^2, Y_i &\sim \textit{Ga}\left(\frac{\nu_1}{2}, \frac{\nu_{2i}}{2}\right), \\
\\
\text{where}\ \ B_1 &= [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1}, \\
\tilde{\beta} &= B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y],\\
\alpha_1 &= \alpha_0 + n, \\
\delta_1 &= \delta_0 + (Y-X\beta)^{T}\Lambda (Y-X\beta), \\
\nu_1 &= \nu + 1, \\
\nu_{2,i} &= \nu + σ^{-2}(Y_i-X_i\beta).
\end{align}
$$
どうでも良い話ですが、はてなブログ上でtex記法を使って\sigmaを打とうとすると、なぞにキーワードにリンクされてしまい、数式がうまく表示されないのでクソすぎると思いました。
*1:実際には、Breusch-Pagan testやWhite testを用いてこのHeteroskedasticityを確認します。
最新のグランジミュージックの動向 ― Modern Grungeについて
突然に音楽の話を。
最近の気になる動向として、海外音楽シーンにおける「グランジリバイバル」があります。文字通りですが、Pearl JamやNirvanaなどのレジェンド的バンドに代表される'Grunge'が、現行の音楽シーンにおいて再解釈・再構築され、新たな形で復活を遂げているです。この動向を語るに際して'Grunge'とは何かを考えなくてはならない気がしますが、こんなニッチ中のニッチのその下レベルに該当する弱小ブログのグランジに関する記事にたどり着く読者の方は、'Grunge'を良く理解していると勝手に思い込んでいますのでスルーします（そもそも音楽のジャンルに厳格な定義を当てはめて議論していくことが反吐が出るレベルで嫌いですが）。
勝手ですが、それらに該当するバンドを'Modern Grunge'と呼んでいます。個々のサウンドを一緒くたにしたり、ジャンルという枠で括る無意味さ、もとい愚かさを感じますが、今はそう呼んでいます。
これらのバンドに一貫して観察される特徴として、Grungeの渋さやアグレッシブさとShoegazerの浮遊感を兼ね備えているバンドが非常に多いという点が挙げられます。その点から、'Grunge Gaze'や'Space Grunge'といういかにも造語チックな呼称で表現することもできるかもしれません（Bandcampに一応タグは存在します）。
Modern Grungeバンドの多くは、多種多様な音楽的文脈的視点からグランジサウンドに対してアプローチを繰り広げているイメージがあります。例えば、先述の通りShoegazer的な視点からグランジを再構築しているバンド、また、前身がHardcoreであるためにHardcore的なアプローチを執っているバンドもあります。中には、Indie Rockの要素も…と、例を挙げてるとキリがありませんので、この辺りで留めておきましょう。このように、多文脈を包含した、というよりはジャンルのクロスオーバーを織りなしているのがこの'Modern Grunge'なのです。
小中学生の頃は狂ったようにPuddle of MuddやLifehouseなどの2000年代のPost-Grungeを聴きまくっていたことがあり、懐かしさからこのグランジブームに強く惹かれるという点もあるかと思いますが、何が私をここまで惹きつけるのかといえば、やはりModern Grungeサウンドが抱擁する多文脈性だと思います。住む世界を異にする（完全に異かと言われると、そうではないと思いますが）サウンド同士が交わることによる危険な化学反応が、やはり何とも言えぬ魅力なのだと思います。
言葉を羅列するだけではあまり伝わらないと思うので、Modern Grungeに直接触れていただきたいと思います。以下に7バンドほど、個人的な独断と偏見で紹介します。長いとは思いますが、Quarantine Playlist だと考えてぜひお付き合いください。
① Soul Blind
NY出身。ハードコア等に立脚するであろうスタイルと，昔懐かしのサウンド，さらには独特の浮遊感というかSpace感というかがクロスオーバーしていて非常に面白い。個人的最重要バンドである。Twitter上でリプライをもらったのですが，Failureに強く影響を受けているとのこと。確かにあのSpace感とベースの唸りはFailureゆずりだなと思う。
Soul Blind - "Crawling Into You"
②Narrow Head
TX出身。最近新譜が出たのだが，これがなかなかというより非常に良すぎる。DeftonesがOleanderと出会った結果みたいなサウンドを醸し出していて非常にユニークである。懐かしサウンドすぎて要所要所で鳥肌の嵐に見舞われる。GodfleshのTeeを着てるのも非常に好感を持てる。しかも，安心安全のRun For Coverである。
Narrow Head - "Stuttering Stanley"
③My Ticket Home
OH出身。このバンド面白いのが，2010年あたりにゴリゴリのTHE RISE RECORD METALCOREをやっていたのである。今，ベースボーカルをしているNickはかつてピンのスクリームボーカルをやっていたのである。そのため，Run For Cover界隈では語られることもないのである。
彼らのLatest Album（といっても2017年だが）は強烈なポストグランジアルバムである。彼らのことを"Nu Metalcore"と形容する人が多いが，どう考えてもPuddle of Mudd直系サウンドのポストグランジなのである。さらに，スクリームや疾走感あふれるリフの嵐とのクロスオーバーで聴いていて非常に気持ちいい。さらに，ベースボーカルNickのルックスがイケすぎており，どこかカート・コバーンっぽさがある。
余談であるが，彼らを聴いていると2006年あたりのHyde（666のアルバムあたり？）のサウンドが想起される。Countdownは特にMy Ticket Homeっぽい。いや，彼らがHydeっぽいのか。
My Ticket Home - "Flypaper"
④NOTHING
PA出身。多分，この記事にたどり着く人には説明不要であろう。シューゲなのかドゥームゆずりなのかは上手く言い表せないが，重厚な音像と多幸感あふれるサウンドがミックスされ，さらにグランジテイストがむんむんにあふれる何ともいいとこどりなエッチなバンド。言うまでもなく，Jesus PieceのAaronもベースとして参加しており，さらにRelapseから出ているという怖すぎるバンドである。
NOTHING - "Vertigo Flowers"
⑤Superheaven
PA出身。安心安全のRun For Cover。大学の空きコマで彼らを聴きまくっているせいで，大学周辺を歩いていると脳内再生されるレベルで頭が洗脳されている。ちなみに，JarのLPもあるので家でも狂ったように聴きまくっている。
サウンドはPure Post Grungeである。が，演奏がかなり無機質かつ謎の寂寥感にあふれており，かつてのポストグランジの文脈とはまた異なったサウンドで興味深い。イントロのリフで使われているチョーキングは，最高のグランジリフである。
Superheaven - "Life In a Jar"
⑥Teenage Wrist
LA出身。GTAVで有名なThe Chain Gang of 1974の人がボーカルをやっていた。今は，確か在籍メンバーが2人しかいない。ギターがボーカルを兼任するようになったらしく，新譜『Silverspoon』を聴くと，今のボーカルの声の方が好きだなと思ってしまった。
雑談は置いておいて，Teenage Wristは名門Epitaphからのバンドである。Epitaphらしからぬサウンドなのも面白いが，ボーカルのウィスパーボイスと声の抜き方が非常に独特でとにかくエロい。ギターもwet全開のリバーブが聴いていて浮遊感があって非常に心地よい。
Teenage Wrist - "Dweeb"
⑦Fleshwater
MA出身。Vein.FMのサイドプロジェクト的な扱いらしい。曲調からは，ドゥームを軸に据えた音像にグランジやエモの香りがする。シューゲイザー要素も強く多幸感あふれているが，ギターのリフがかなりヘヴィーかつアグレッシブでコントラストが非常に面白い。個人的には，Nothingの一歩先のヘヴィネスを据えたバンドだと勝手に解釈している。
以上で切り上げますが，後日また適当に音楽紹介でもしようかと思います。
イントロダクション的な何か。
初めまして。
都内で大学生をしています。
自分が学んだ勉強に関する事項や趣味である音楽や服、読書について垂れ流す予定です。
まあ、感覚としては個人的な備忘録にあたる何らかの思考の軌跡です。
では。
