Twitter で行った4択アンケートの結果を本気で分析する（打ち切りと丸めによって生成されたヒストグラムのモデリング）

経緯

実家のお風呂の温度が一人暮らしと比べて3℃も高いので長湯できない
— しょラー (@shora_kujira16) 2021年12月28日

お風呂の温度
— しょラー (@shora_kujira16) 2021年12月28日

事前の想像では最頻値は 40 or 41 になるのではないかと考えていたため Twitter のアンケート機能で利用可能な4択を上記のように設定したのですが、予想が外れて「42℃以上」が最頻値となってしまいました。これでは結局何℃にしている人が一番多いのかよくわかりません。

アンケートの回答の総数が 28 だったため、割合から逆算して以下のようなデータが得られました。

温度	回答数
39℃以下	2
40℃	5
41℃	10
42℃以上	11

分析で知りたいこと

一般のご家庭の平均的なお風呂の設定温度を知りたい
お風呂の温度を43℃に設定しているご家庭がどれくらいあるのか知りたい

愚直な分析で困ること

打ち切り

39℃以下を39、42℃以上を42と扱っても、それっぽい平均は計算できます。しかし「42℃以上」と回答した人は11人もいるので43℃や44℃に設定している人も何人かいそうです。よって真の平均はそれより少し高くなりそうです。

丸め

各家庭が本当に好むお風呂の温度は本来実数でなければ表現できないはずですが、アンケートの都合により回答を整数で入力させました。これにより、たとえば41.6℃が最適な人でも「42℃以上」と回答することになります。この丸めによって平均値がずれることがありそうです。

仮定

各家庭が本当に好むお風呂の温度が実数で回答できて打ち切りがない場合にアンケートの結果が正規分布に従うと仮定します。

モデル

お風呂の温度の好みが平均 $\mu$ 、分散 $\sigma^2$ の正規分布から生成され、打ち切りと丸めによって 39.5, 40.5, 41.5 を区切りとして「39℃以下」「40℃」「41℃」「42℃以上」に振り分けられるとします。このときそれぞれの項目に割り振られる確率 $\theta$ は累積分布関数 (normal_cdf) を使って表すことができます。

$\theta$ が分かっていて回答の総数が $N$ のとき、所望のアンケート結果 $y$ が得られる確率は多項分布 $\mathrm{Multinomial}(y|\theta)$ として表すことができます。今回はこの考え方に従って Stan の multinomial 関数を使いましたが、多項係数の部分は尤度の計算においては推定したい未知パラメータの MCMC に対して不変のため、（私の理解が間違っていなければ）必ずしも使う必要はないはずです。

Stan のモデルコードは以下のとおりです。

data {
  int<lower=0> Y[4];
}

parameters {
  real mu;
  real<lower=0> sigma;
}

transformed parameters {
  vector<lower=0>[4] theta;
  theta[1] = normal_cdf(39.5 | mu, sigma);
  theta[2] = normal_cdf(40.5 | mu, sigma) - normal_cdf(39.5 | mu, sigma);
  theta[3] = normal_cdf(41.5 | mu, sigma) - normal_cdf(40.5 | mu, sigma);
  theta[4] = 1 - normal_cdf(41.5 | mu, sigma);
}

model {
  Y ~ multinomial(theta);
}

generated quantities {
  int<lower=0> Y_pred[4];
  Y_pred = multinomial_rng(theta, sum(Y));
}