Stanのパラメータ制限値を複数指定する
参考
やりたかったこと
渡された複数枚の図形を指定の順に並べる課題の課題遂行時間を分析したかった。この課題ではどれだけ習熟しても，図形を並び替えるのに最低限の時間（以下，ベースタイム）が存在すると考えられる。
課題遂行時間のデータは上のグラフのようになっている。このデータについて，試行が進むと，課題遂行時間が群ごとに異なるベースタイムに急速に寄っていくモデルを考える。この時，ベースタイムというパラメータの性質上，各群でその推定値は最短の課題遂行時間を超えてはならない。そこで，各群で，ベースラインのパラメータ推定値の上限を最短の課題遂行時間に設定したい。
data_raw <- read.csv("data.csv")
Trials <- array(c(data_raw$Trials),c(120, 4))
Time <- array(c(data_raw$Time),c(120, 4))
Time_min <- c(min(subset(data_raw$Time, A==1&B==1, c(Time))),
min(subset(data_raw$Time, A==1&B==2, c(Time))),
min(subset(data_raw$Time, A==2&B==1, c(Time))),
min(subset(data_raw$Time, A==2&B==2, c(Time))))
data_Time <- list(N=120,
Group=4,
Trials=Trials,
Time=Time,
Time_min=Time_min)> str(data_Time)
List of 5
$ N : num 120
$ Group : num 4
$ Trials : int [1:120, 1:4] 1 2 3 4 5 6 1 2 3 4 ...
$ Time : num [1:120, 1:4] 416 284 156 113 152 ...
$ Time_min: num [1:4] 47.2 54.1 28.7 19.1 なのでまず，何も考えずに上のようにRで最小値のベクトルを含むデータを作り，
data{
int N;
int Group;
int<lower=1, upper=6> Trials[N, Group];
real<lower=0> Time[N, Group];
vector[Group] Time_min;
}
parameters{
vector<lower=0, upper=Time_min[Group]>[Group] base_time; //ベースタイム
vector<lower=0>[Group] a; //1試行目でベースタイムに加算される時間
vector<lower=0>[Group] sigma; //分散
}
model{
base_time[Group] ~ normal(0, 100);
a[Group] ~ normal(0, 100);
sigma[Group] ~ cauchy(0, 100);
for(gr in 1:Group){
for(n in 1:N){
Time[n, gr] ~ normal(base_time[gr] + a[gr]*exp(-(Trials[n, gr]-1)), sigma[gr]);
}
}
}上のモデルで分析した結果が以下である。
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
base_time[1] 14.71 0.06 3.99 3.99 12.84 15.98 17.78 18.94 4640 1.00
base_time[2] 16.28 0.04 2.71 8.80 15.17 17.11 18.22 18.99 3732 1.00
base_time[3] 14.18 0.06 4.32 2.56 12.06 15.43 17.49 18.99 4928 1.00
base_time[4] 17.16 0.02 1.78 12.34 16.40 17.69 18.49 19.01 6178 1.00
a[1] 637.18 0.69 53.70 537.89 600.39 637.21 673.03 743.27 6039 1.00
a[2] 444.43 0.44 37.39 370.48 419.57 444.68 468.36 519.92 7268 1.00
a[3] 672.12 0.54 43.35 589.89 642.94 671.05 701.23 756.93 6491 1.00
a[4] 294.82 0.17 14.11 268.06 285.10 294.75 304.26 322.73 6499 1.00
sigma[1] 247.61 0.19 16.00 217.99 236.69 247.11 257.53 281.49 7104 1.00
sigma[2] 175.69 0.15 11.54 154.98 167.43 174.86 183.33 199.97 5967 1.00
sigma[3] 202.85 0.17 13.14 179.32 193.65 202.00 211.25 230.80 6202 1.00
sigma[4] 68.46 0.05 4.46 60.35 65.39 68.17 71.39 77.78 7083 1.00
lp__ -2621.75 0.07 2.68 -2627.77 -2623.31 -2621.41 -2619.75 -2617.55 1550 1.01ベースタイムの推定値を見ると，どうもベクトルで指定した上限値（47.2，54.1，28.7，19.1）の内、最小の値（19.1）だけが単一の上限値として認識されたようである。
なんで？
解決方法
まず，Stanのユーザーガイドによると上でやったような指定は出来ないもんらしい。しかし同時に，transformed parameterブロックを用いて，手動で最大値（あるいは最小値）を設定することで対応可能であり，わざわざ直感的に最初に書くであろう動かないコードと，それに対応する実際に動くコードを示してくれていた。
data{
int N;
int Group;
int<lower=1, upper=6> Trials[N, Group];
real<lower=0> Time[N, Group];
vector[Group] Time_min;
}
parameters{
vector<lower=0, upper=1>[Group] base_time_raw;
vector<lower=0>[Group] a;
vector<lower=0>[Group] sigma;
}
transformed parameters{
vector<lower=0>[Group] base_time = Time_min .* base_time_raw;
}
model{
base_time_raw[Group] ~ normal(0, 1);
a[Group] ~ normal(0, 100);
sigma[Group] ~ cauchy(0, 100);
for(gr in 1:Group){
for(n in 1:N){
Time[n, gr] ~ normal(base_time[gr] + a[gr]*exp(-(Trials[n, gr]-1)), sigma[gr]);
}
}
}
ユーザーガイドを参考にして書き直したモデルが上のものである。分析し直すと以下のようになった。
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
base_time_raw[1] 0.90 0.00 0.10 0.63 0.85 0.93 0.97 1.00 5694 1
base_time_raw[2] 0.94 0.00 0.06 0.79 0.93 0.96 0.98 1.00 5676 1
base_time_raw[3] 0.81 0.00 0.18 0.34 0.72 0.86 0.94 0.99 5027 1
base_time_raw[4] 0.90 0.00 0.10 0.64 0.86 0.93 0.97 1.00 5955 1
a[1] 599.84 0.70 49.32 501.97 566.66 599.40 633.26 697.53 5021 1
a[2] 396.13 0.47 33.07 331.59 373.76 395.72 417.86 461.71 4934 1
a[3] 660.75 0.54 41.44 578.57 632.91 660.56 688.97 740.82 5810 1
a[4] 294.51 0.20 14.15 265.96 285.24 294.58 303.88 322.57 4940 1
sigma[1] 235.87 0.20 16.21 206.80 224.65 234.95 246.23 270.75 6264 1
sigma[2] 158.30 0.15 10.62 139.23 150.69 157.71 165.21 180.71 5310 1
sigma[3] 200.37 0.18 13.57 175.94 191.02 199.48 209.03 229.01 5443 1
sigma[4] 68.40 0.06 4.59 60.07 65.21 68.15 71.37 78.06 5684 1
base_time[1] 42.22 0.06 4.90 29.85 40.30 43.67 45.70 47.02 5694 1
base_time[2] 51.13 0.04 3.05 42.84 50.07 52.03 53.28 54.04 5676 1
base_time[3] 23.14 0.07 5.15 9.71 20.68 24.71 27.08 28.59 5027 1
base_time[4] 17.18 0.02 1.84 12.27 16.41 17.74 18.52 19.01 5955 1
lp__ -2615.60 0.06 2.56 -2621.35 -2617.16 -2615.27 -2613.78 -2611.50 1696 1ベースタイムの推定値を見るに，無事それぞれの群ごとに上限値が設定されたようである。
追記
やってはみたが，そもそもベースタイムは人によって異なるので，ベースラインのパラメータ推定値の上限を群毎の最短の課題遂行時間に制限するのは無理があることに気付いてしまった。群毎に正規分布で平均+分散という形で推定するべきだと思う。