UPDATE: 2024-01-05 14:14:28.09333

はじめに

このノートは「StanとRでベイズ統計モデリング」の内容を写経することで、ベイズ統計への理解を深めていくために作成している。

基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、「StanとRでベイズ統計モデリング」を読み進めるための自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。

今回は第11章「離散値をとるパラメータを使う」の後半から写経していく。

11.3 ゼロ過剰ポアソン分布

ある飲食店への来店回数を取得したデータをもとに来店回数が多い人の特徴を調べたい。来店回数はポアソン分布を仮定できそうではあるが、0が多いことから、ゼロ過剰ポアソン分布を利用するほうが望ましい。

  • Sex: 0が男性、1が女性
  • Sake: 0が飲まない、1が飲む
library(dplyr)
library(rstan)
library(ggplot2)
library(GGally)

options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap11/input/data-ZIP.txt')
head(d)
##   Sex Sake Age Y
## 1   0    1  18 5
## 2   1    0  18 2
## 3   1    1  18 1
## 4   0    0  19 3
## 5   0    0  19 5
## 6   1    0  19 7

11.3.1 解析目的とデータの分布の確認

解析の目的は「リピーターになりそうな人を知りたい」とする。また「説明変数でリピーターになるかをどれほど予測できるかを知りたい、影響度合いも知りたい」とする。

  • 1行1列目のグラフから、女性が少なそう
  • 2行2列目のグラフから、お酒を飲む人が少なそう
  • 3行3列目のグラフから、年齢に特徴はなさそう
  • 4行4列目のグラフから、来店回数が0が突出して多く、0を除くと正規分布のような形状
  • 4行1列目のグラフから、女性だと来店回数が少なそう
  • 4行3列目のグラフから、年齢が高いと来店回数が多くなりそう
  • 2行4列目のグラフから、お酒を飲むかつ来店回数が0の人が少ないため順位相関係数が少し高い
d$Sex <- as.factor(d$Sex)
d$Sake <- as.factor(d$Sake)

N_col <- ncol(d)
ggp <- ggpairs(d, upper = 'blank', diag = 'blank', lower = 'blank')

# 対角成分のヒストグラムを作成
for (i in 1:N_col) {
  x <- d[,i]
  p <- ggplot(data.frame(x), aes(x = x)) +
    theme_bw(base_size = 14) +
    theme(axis.text.x = element_text(angle = 40, vjust = 1, hjust = 1))
  if (class(x) == 'factor') {
    p <- p + geom_bar(color = 'grey5', fill = 'grey80')
  } else {
    bw <- ifelse(colnames(d)[i] == 'Y', 1, (max(x)-min(x))/10)
    p <- p + geom_histogram(binwidth = bw, color = 'grey5', fill = 'grey80') +
      geom_line(aes(y = after_stat(count)*bw), stat = 'density')
  }
  p <- p + geom_label(data = data.frame(x = -Inf, y = Inf, label = colnames(d)[i]), aes(x = x, y = y, label = label), hjust = 0, vjust = 1)
  ggp <- putPlot(ggp, p, i, i)
}

# 上三角成分の相関係数を作成
zcolat <- seq(-1, 1, length = 81)
zcolre <- c(zcolat[1:40]+1, rev(zcolat[41:81]))

for (i in 1:(N_col-1)) {
  for (j in (i+1):N_col) {
    x <- as.numeric(d[,i])
    y <- as.numeric(d[,j])
    r <- cor(x, y, method = 'spearman', use = 'pairwise.complete.obs')
    zcol <- lattice::level.colors(r, at = zcolat, col.regions = grey(zcolre))
    textcol <- ifelse(abs(r) < 0.4, 'grey20', 'white')
    ell <- ellipse::ellipse(r, level = 0.95, type = 'l', npoints = 50, scale = c(.2, .2), centre = c(.5, .5))
    p <- ggplot(data.frame(ell), aes(x = x, y = y)) + theme_bw() + theme(
      plot.background = element_blank(),
      panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
      panel.border = element_blank(), axis.ticks = element_blank()) +
      geom_polygon(fill = zcol, color = zcol) +
      geom_text(data = NULL, x = .5, y = .5, label = 100*round(r, 2), size = 6, col = textcol)
    ggp <- putPlot(ggp, p, i, j)
  }
}

# 下三角成分の箱ヒゲ図を作成
for (j in 1:(N_col-1)) {
  for (i in (j+1):N_col) {
    x <- d[,j]
    y <- d[,i]
    if (class(x) == 'factor' && class(y) == 'factor') {
      p <- ggplot(reshape2::melt(table(x,y)), aes(x = x, y = y)) +
        theme_bw(base_size = 14) +
        geom_point(aes(size = value), color = 'grey80') +
        geom_text(aes(label = value)) +
        scale_size_area(max_size = 8) +
        scale_x_continuous(breaks = 0:1, limits = c(-0.5,1.5)) + scale_y_continuous(breaks = 0:1, limits = c(-0.5,1.5))
    } else {
      p <- ggplot(data.frame(x, y, Y = as.numeric(d$Y)), aes(x = x, y = y, color = Y)) +
        theme_bw(base_size = 14) +
        theme(axis.text.x = element_text(angle = 40, vjust = 1, hjust = 1))
      if (class(x) == 'factor') {
        p <- p + geom_boxplot(alpha = 3/6, outlier.size = 0, fill = 'white') +
          geom_point(position = position_jitter(w = 0.4, h = 0), size = 1)
      } else {
        p <- p + geom_point(size = 1)
      }
      p <- p + scale_color_gradient(low = 'grey65', high = 'grey5')
    }
    ggp <- putPlot(ggp, p, i, j)
  }
}
print(ggp, left = 0.3, bottom = 0.3)

11.3.2 メカニズムの想像

来店回数\(Y\)の0が突出している点が特徴的なので、このメカニズムを考える。1人が来店するメカニズムとして、「とにかく1回来店する」場合と「気に入って複数回来店する」場合のパターンがある。そして、とにかく1回来店するのは確率\(q\)のベルヌイ分布に従うと考え、リピーターの来店回数は平均\(\lambda\)のポアソン分布に従うと考える。

つまり、確率\(q\)のコインを投げて、裏が出れば来店回数は0、表が出れば平均\(\lambda\)のポアソン分布に従って生成される。そして、\(q\)はロジスティック回帰、\(\lambda\)はポアソン回帰でパラメタの推定を分ける。コイン投げは離散パラメタとなるので、周辺化消去が必要になる。

\(y=0\)のときは「ベルヌイ分布で0が出る確率」と「ベルヌイ分布で1が出る確率」と「ポアソン分布で0が出る確率」の和となる。\(y=0\)のときはベルヌイ分布の可能性が高いが、\(y=1\)のときはどちらの分布か判断し難くく、\(y\ge2\)のときはポアソン分布に従うと考える。そのため、\(Bernoulli(0|q) + Bernoulli(1|q) × Poisson(y=0|\lambda)\)とする。

\[ \begin{eqnarray} ZIP(y|q,\lambda) = \begin{cases} Bernoulli(0|q) + Bernoulli(1|q) × Poisson(y=0|\lambda) \ if \ y = 0 \\ Bernoulli(1|q) × Poisson(y|\lambda) \ if \ y \ge 1 \\ \end{cases} \end{eqnarray} \]

11.3.3 モデル式の記述

\(\boldsymbol{ X }\)はN×4の行列で、\(\overrightarrow{ b}\)は長さ4のベクトル。

モデル11-7

\[ \begin{eqnarray} q[n] &=& inv\_logit((\boldsymbol{ X } \overrightarrow{ b_{1}})[n]) \\ \lambda[n] &=& (\boldsymbol{ X } \overrightarrow{ b_{2}})[n] \\ Y[n] &\sim& ZIP(q[n], \lambda[n]) \end{eqnarray} \]

11.3.4 Stanで実装

Stanのモデルは下記の通り。

functions {
  real ZIP_lpmf(int Y, real q, real lambda) {
    if (Y == 0) {
      return log_sum_exp(
        bernoulli_lpmf(0 | q),
        bernoulli_lpmf(1 | q) + poisson_log_lpmf(0 | lambda)
      );
    } else {
      return bernoulli_lpmf(1 | q) + poisson_log_lpmf(Y | lambda);
    }
  }
}

data {
  int N;
  int D;
  int<lower=0> Y[N];
  matrix[N,D] X;
}

parameters {
  vector[D] b[2];
}

transformed parameters {
  vector[N] q_x;
  vector[N] q;
  vector[N] lambda;

  q_x = X*b[1];
  lambda = X*b[2];
  for (n in 1:N)
    q[n] = inv_logit(q_x[n]);
}

model {
  for (n in 1:N)
    Y[n] ~ ZIP(q[n], lambda[n]);
}

11.3.5 推定結果の解釈

データを用意する。

d$Age <- d$Age/10
X <- cbind(1, d[,-ncol(d)])
data <- list(N = nrow(d), D = ncol(X), Y = d$Y, X = X)
data
## $N
## [1] 200
## 
## $D
## [1] 4
## 
## $Y
##   [1]  5  2  1  3  5  7  5  0  8  8  0  7  6  0  0  1  4  0  0  3  0  0  4  8  7
##  [26]  3  7  0  2  4  0  9  0  7  0  0  9  4  6  1 13  5  7  1  6  2 10  0  2  1
##  [51]  6  9  0  9  8  3 11  4  0  6  3  1  4  3  0  0  5  0  8  7  2  4  0  7  6
##  [76]  5  0  0  0  4  2 11 10  5  0  0  0  0  0  2  7 11  0  3  6  4  1  0  0  5
## [101]  0  0  6  0  5  9  0  0 10  0  2  0  6 10  6  8  5  8  9  3  6  5  9  0  0
## [126] 11  3  0  5  6  6  7  9  2  0 12  0  0  5  4  0  0  7  0  9  6  0  4  8  0
## [151]  0  4  0  8  6  3  8 11  0  7  0 15  0  0  0  9  0 12  0  9  7  9  0 15  9
## [176]  6  3  0  0  0 12  2  9  5 17  6  0  0 10 15  0  0  5  4 13  6  0  0 12  6
## 
## $X
##     1 Sex Sake Age
## 1   1   0    1 1.8
## 2   1   1    0 1.8
## 3   1   1    1 1.8
## 4   1   0    0 1.9
## 5   1   0    0 1.9
## 6   1   1    0 1.9
## 7   1   0    1 2.0
## 8   1   0    0 2.0
## 9   1   0    0 2.1
## 10  1   0    0 2.1
## 11  1   0    0 2.1
## 12  1   0    1 2.1
## 13  1   0    0 2.1
## 14  1   1    0 2.2
## 15  1   0    0 2.2
## 16  1   1    0 2.2
## 17  1   0    1 2.2
## 18  1   0    0 2.2
## 19  1   1    0 2.3
## 20  1   1    0 2.3
## 21  1   0    0 2.3
## 22  1   1    0 2.3
## 23  1   1    0 2.3
## 24  1   1    0 2.3
## 25  1   0    1 2.3
## 26  1   0    0 2.3
## 27  1   0    0 2.3
## 28  1   0    0 2.4
## 29  1   1    0 2.4
## 30  1   1    1 2.4
## 31  1   1    1 2.5
## 32  1   0    0 2.5
## 33  1   0    0 2.6
## 34  1   0    1 2.6
## 35  1   0    0 2.6
## 36  1   0    0 2.6
## 37  1   0    0 2.6
## 38  1   1    0 2.6
## 39  1   0    0 2.6
## 40  1   1    0 2.6
## 41  1   0    0 2.6
## 42  1   1    0 2.7
## 43  1   0    1 2.7
## 44  1   1    0 2.7
## 45  1   0    1 2.7
## 46  1   1    1 2.7
## 47  1   0    1 2.8
## 48  1   1    0 2.8
## 49  1   1    0 2.8
## 50  1   1    0 2.8
## 51  1   0    0 2.8
## 52  1   0    0 2.8
## 53  1   0    0 2.8
## 54  1   0    0 2.8
## 55  1   0    0 2.9
## 56  1   1    1 2.9
## 57  1   0    1 2.9
## 58  1   0    1 2.9
## 59  1   0    0 3.0
## 60  1   1    0 3.0
## 61  1   1    1 3.0
## 62  1   1    0 3.0
## 63  1   1    0 3.0
## 64  1   1    0 3.0
## 65  1   0    0 3.1
## 66  1   0    0 3.1
## 67  1   0    0 3.1
## 68  1   0    0 3.1
## 69  1   0    0 3.2
## 70  1   0    1 3.2
## 71  1   1    0 3.2
## 72  1   0    1 3.2
## 73  1   0    0 3.2
## 74  1   0    1 3.2
## 75  1   1    0 3.2
## 76  1   1    0 3.2
## 77  1   1    0 3.2
## 78  1   0    0 3.3
## 79  1   0    0 3.3
## 80  1   1    1 3.3
## 81  1   1    0 3.3
## 82  1   0    1 3.3
## 83  1   0    0 3.3
## 84  1   0    0 3.3
## 85  1   0    1 3.3
## 86  1   0    0 3.3
## 87  1   0    0 3.4
## 88  1   0    0 3.4
## 89  1   0    0 3.4
## 90  1   1    1 3.4
## 91  1   0    1 3.4
## 92  1   0    0 3.5
## 93  1   0    0 3.5
## 94  1   1    1 3.5
## 95  1   1    0 3.5
## 96  1   1    1 3.5
## 97  1   1    0 3.5
## 98  1   0    0 3.6
## 99  1   0    0 3.6
## 100 1   1    1 3.6
## 101 1   1    0 3.6
## 102 1   0    0 3.6
## 103 1   1    0 3.6
## 104 1   1    0 3.7
## 105 1   1    0 3.7
## 106 1   0    0 3.7
## 107 1   0    0 3.8
## 108 1   0    0 3.8
## 109 1   0    1 3.8
## 110 1   0    0 3.8
## 111 1   0    0 3.8
## 112 1   0    0 3.8
## 113 1   1    1 3.9
## 114 1   0    0 3.9
## 115 1   0    1 3.9
## 116 1   1    0 3.9
## 117 1   1    0 3.9
## 118 1   0    1 4.0
## 119 1   0    0 4.0
## 120 1   1    1 4.0
## 121 1   1    0 4.0
## 122 1   0    0 4.0
## 123 1   0    0 4.1
## 124 1   1    0 4.1
## 125 1   0    0 4.1
## 126 1   0    0 4.1
## 127 1   1    1 4.1
## 128 1   0    0 4.1
## 129 1   0    1 4.1
## 130 1   1    0 4.2
## 131 1   1    1 4.2
## 132 1   0    0 4.2
## 133 1   1    0 4.2
## 134 1   1    1 4.2
## 135 1   0    0 4.2
## 136 1   0    1 4.3
## 137 1   1    0 4.3
## 138 1   1    1 4.3
## 139 1   1    0 4.3
## 140 1   1    0 4.3
## 141 1   0    1 4.3
## 142 1   0    0 4.3
## 143 1   0    1 4.4
## 144 1   0    0 4.4
## 145 1   0    1 4.4
## 146 1   0    1 4.4
## 147 1   1    0 4.4
## 148 1   0    1 4.4
## 149 1   0    1 4.5
## 150 1   1    0 4.5
## 151 1   0    0 4.6
## 152 1   1    0 4.6
## 153 1   0    0 4.6
## 154 1   0    1 4.6
## 155 1   1    0 4.6
## 156 1   1    0 4.6
## 157 1   1    0 4.7
## 158 1   0    0 4.7
## 159 1   0    0 4.7
## 160 1   1    0 4.7
## 161 1   0    0 4.8
## 162 1   0    1 4.8
## 163 1   1    1 4.8
## 164 1   0    0 4.8
## 165 1   0    0 4.8
## 166 1   1    0 4.8
## 167 1   0    0 4.9
## 168 1   0    1 4.9
## 169 1   0    0 5.0
## 170 1   0    1 5.0
## 171 1   0    0 5.0
## 172 1   0    1 5.0
## 173 1   1    0 5.1
## 174 1   0    0 5.1
## 175 1   0    1 5.1
## 176 1   1    1 5.1
## 177 1   1    0 5.1
## 178 1   0    0 5.1
## 179 1   0    0 5.1
## 180 1   1    0 5.1
## 181 1   0    1 5.2
## 182 1   1    1 5.2
## 183 1   0    1 5.2
## 184 1   1    0 5.2
## 185 1   0    1 5.2
## 186 1   1    0 5.2
## 187 1   0    0 5.3
## 188 1   0    0 5.3
## 189 1   0    1 5.3
## 190 1   0    0 5.4
## 191 1   1    0 5.4
## 192 1   0    0 5.4
## 193 1   1    1 5.4
## 194 1   1    0 5.4
## 195 1   0    0 5.4
## 196 1   0    1 5.5
## 197 1   0    0 5.5
## 198 1   0    0 5.5
## 199 1   0    1 5.5
## 200 1   1    0 5.5

ここでは、stan_model()関数で最初にコンパイルしておいてから、

model117 <- stan_model('note_ahirubayes15-117.stan')

sampling()関数でサンプリングする。

fit <- sampling(object = model117, data = data, seed = 1989)

推定結果は下記の通り。b[1,]は来店確率に関わる回帰係数で、b[2,]はリピーターが何回来店するかに関わる回帰係数。b[,1−4]は切片、SexSakeAgeに対応する。Ageは10で割っているので、10倍する必要がある。1回来る来店確率が高いのは、「女性(b[1,2]=1.6)でお酒を飲んでいる人(b[1,3]=3.4)」である。リピーターになる確率が高いのは、「男性(b[2,2]=-0.7)でお酒を飲んでいない人(b[2,3]=-0.2)で、年齢が高い人(b[2,4]=0.2)」である。

print(fit, prob = c(0.025, 0.5, 0.975), pars = c('b', 'q', 'lambda'), digits_summary = 1)
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##             mean se_mean  sd 2.5%  50% 97.5% n_eff Rhat
## b[1,1]      -4.0       0 1.3 -6.6 -3.9  -1.7  2089    1
## b[1,2]       1.6       0 0.4  0.8  1.6   2.5  2918    1
## b[1,3]       3.4       0 0.9  2.0  3.2   5.4  2121    1
## b[1,4]      -0.4       0 0.2 -0.7 -0.4   0.0  3029    1
## b[2,1]       2.4       0 0.2  2.0  2.4   2.7  2065    1
## b[2,2]      -0.7       0 0.1 -0.9 -0.7  -0.6  2728    1
## b[2,3]      -0.2       0 0.1 -0.3 -0.2   0.0  2408    1
## b[2,4]       0.2       0 0.0  0.1  0.2   0.3  2618    1
## q[1]         1.0       0 0.0  0.9  1.0   1.0  3471    1
## q[2]         0.9       0 0.1  0.7  0.9   1.0  2949    1
## q[3]         1.0       0 0.0  1.0  1.0   1.0  2623    1
## q[4]         0.6       0 0.1  0.4  0.6   0.7  3198    1
## q[5]         0.6       0 0.1  0.4  0.6   0.7  3198    1
## q[6]         0.9       0 0.1  0.7  0.9   0.9  2964    1
## q[7]         1.0       0 0.0  0.9  1.0   1.0  3495    1
## q[8]         0.6       0 0.1  0.4  0.6   0.7  3213    1
## q[9]         0.5       0 0.1  0.4  0.5   0.7  3229    1
## q[10]        0.5       0 0.1  0.4  0.5   0.7  3229    1
## q[11]        0.5       0 0.1  0.4  0.5   0.7  3229    1
## q[12]        1.0       0 0.0  0.9  1.0   1.0  3505    1
## q[13]        0.5       0 0.1  0.4  0.5   0.7  3229    1
## q[14]        0.8       0 0.1  0.7  0.8   0.9  3011    1
## q[15]        0.5       0 0.1  0.4  0.5   0.7  3246    1
## q[16]        0.8       0 0.1  0.7  0.8   0.9  3011    1
## q[17]        1.0       0 0.0  0.9  1.0   1.0  3516    1
## q[18]        0.5       0 0.1  0.4  0.5   0.7  3246    1
## q[19]        0.8       0 0.1  0.7  0.8   0.9  3027    1
## q[20]        0.8       0 0.1  0.7  0.8   0.9  3027    1
## q[21]        0.5       0 0.1  0.4  0.5   0.7  3250    1
## q[22]        0.8       0 0.1  0.7  0.8   0.9  3027    1
## q[23]        0.8       0 0.1  0.7  0.8   0.9  3027    1
## q[24]        0.8       0 0.1  0.7  0.8   0.9  3027    1
## q[25]        1.0       0 0.0  0.9  1.0   1.0  3525    1
## q[26]        0.5       0 0.1  0.4  0.5   0.7  3250    1
## q[27]        0.5       0 0.1  0.4  0.5   0.7  3250    1
## q[28]        0.5       0 0.1  0.4  0.5   0.7  3254    1
## q[29]        0.8       0 0.1  0.7  0.8   0.9  3044    1
## q[30]        1.0       0 0.0  1.0  1.0   1.0  2597    1
## q[31]        1.0       0 0.0  1.0  1.0   1.0  2611    1
## q[32]        0.5       0 0.1  0.4  0.5   0.6  3259    1
## q[33]        0.5       0 0.1  0.4  0.5   0.6  3264    1
## q[34]        1.0       0 0.0  0.9  1.0   1.0  3547    1
## q[35]        0.5       0 0.1  0.4  0.5   0.6  3264    1
## q[36]        0.5       0 0.1  0.4  0.5   0.6  3264    1
## q[37]        0.5       0 0.1  0.4  0.5   0.6  3264    1
## q[38]        0.8       0 0.1  0.7  0.8   0.9  3079    1
## q[39]        0.5       0 0.1  0.4  0.5   0.6  3264    1
## q[40]        0.8       0 0.1  0.7  0.8   0.9  3079    1
## q[41]        0.5       0 0.1  0.4  0.5   0.6  3264    1
## q[42]        0.8       0 0.1  0.7  0.8   0.9  3097    1
## q[43]        1.0       0 0.0  0.9  1.0   1.0  3552    1
## q[44]        0.8       0 0.1  0.7  0.8   0.9  3097    1
## q[45]        1.0       0 0.0  0.9  1.0   1.0  3552    1
## q[46]        1.0       0 0.0  1.0  1.0   1.0  2657    1
## q[47]        1.0       0 0.0  0.9  1.0   1.0  3556    1
## q[48]        0.8       0 0.1  0.7  0.8   0.9  3115    1
## q[49]        0.8       0 0.1  0.7  0.8   0.9  3115    1
## q[50]        0.8       0 0.1  0.7  0.8   0.9  3115    1
## q[51]        0.5       0 0.1  0.4  0.5   0.6  3276    1
## q[52]        0.5       0 0.1  0.4  0.5   0.6  3276    1
## q[53]        0.5       0 0.1  0.4  0.5   0.6  3276    1
## q[54]        0.5       0 0.1  0.4  0.5   0.6  3276    1
## q[55]        0.5       0 0.1  0.4  0.5   0.6  3283    1
## q[56]        1.0       0 0.0  1.0  1.0   1.0  2654    1
## q[57]        1.0       0 0.0  0.9  1.0   1.0  3559    1
## q[58]        1.0       0 0.0  0.9  1.0   1.0  3559    1
## q[59]        0.5       0 0.1  0.3  0.5   0.6  3289    1
## q[60]        0.8       0 0.1  0.7  0.8   0.9  3154    1
## q[61]        1.0       0 0.0  1.0  1.0   1.0  2647    1
## q[62]        0.8       0 0.1  0.7  0.8   0.9  3154    1
## q[63]        0.8       0 0.1  0.7  0.8   0.9  3154    1
## q[64]        0.8       0 0.1  0.7  0.8   0.9  3154    1
## q[65]        0.5       0 0.1  0.3  0.5   0.6  3297    1
## q[66]        0.5       0 0.1  0.3  0.5   0.6  3297    1
## q[67]        0.5       0 0.1  0.3  0.5   0.6  3297    1
## q[68]        0.5       0 0.1  0.3  0.5   0.6  3297    1
## q[69]        0.4       0 0.1  0.3  0.4   0.6  3304    1
## q[70]        0.9       0 0.0  0.9  1.0   1.0  3561    1
## q[71]        0.8       0 0.1  0.7  0.8   0.9  3193    1
## q[72]        0.9       0 0.0  0.9  1.0   1.0  3561    1
## q[73]        0.4       0 0.1  0.3  0.4   0.6  3304    1
## q[74]        0.9       0 0.0  0.9  1.0   1.0  3561    1
## q[75]        0.8       0 0.1  0.7  0.8   0.9  3193    1
## q[76]        0.8       0 0.1  0.7  0.8   0.9  3193    1
## q[77]        0.8       0 0.1  0.7  0.8   0.9  3193    1
## q[78]        0.4       0 0.1  0.3  0.4   0.5  3311    1
## q[79]        0.4       0 0.1  0.3  0.4   0.5  3311    1
## q[80]        1.0       0 0.0  1.0  1.0   1.0  2617    1
## q[81]        0.8       0 0.1  0.7  0.8   0.9  3213    1
## q[82]        0.9       0 0.0  0.9  1.0   1.0  3559    1
## q[83]        0.4       0 0.1  0.3  0.4   0.5  3311    1
## q[84]        0.4       0 0.1  0.3  0.4   0.5  3311    1
## q[85]        0.9       0 0.0  0.9  1.0   1.0  3559    1
## q[86]        0.4       0 0.1  0.3  0.4   0.5  3311    1
## q[87]        0.4       0 0.1  0.3  0.4   0.5  3318    1
## q[88]        0.4       0 0.1  0.3  0.4   0.5  3318    1
## q[89]        0.4       0 0.1  0.3  0.4   0.5  3318    1
## q[90]        1.0       0 0.0  1.0  1.0   1.0  2605    1
## q[91]        0.9       0 0.0  0.9  0.9   1.0  3556    1
## q[92]        0.4       0 0.1  0.3  0.4   0.5  3325    1
## q[93]        0.4       0 0.1  0.3  0.4   0.5  3325    1
## q[94]        1.0       0 0.0  1.0  1.0   1.0  2591    1
## q[95]        0.8       0 0.1  0.6  0.8   0.9  3251    1
## q[96]        1.0       0 0.0  1.0  1.0   1.0  2591    1
## q[97]        0.8       0 0.1  0.6  0.8   0.9  3251    1
## q[98]        0.4       0 0.1  0.3  0.4   0.5  3331    1
## q[99]        0.4       0 0.1  0.3  0.4   0.5  3331    1
## q[100]       1.0       0 0.0  1.0  1.0   1.0  2582    1
## q[101]       0.8       0 0.1  0.6  0.8   0.9  3269    1
## q[102]       0.4       0 0.1  0.3  0.4   0.5  3331    1
## q[103]       0.8       0 0.1  0.6  0.8   0.9  3269    1
## q[104]       0.8       0 0.1  0.6  0.8   0.9  3286    1
## q[105]       0.8       0 0.1  0.6  0.8   0.9  3286    1
## q[106]       0.4       0 0.1  0.3  0.4   0.5  3336    1
## q[107]       0.4       0 0.1  0.3  0.4   0.5  3340    1
## q[108]       0.4       0 0.1  0.3  0.4   0.5  3340    1
## q[109]       0.9       0 0.0  0.8  0.9   1.0  3478    1
## q[110]       0.4       0 0.1  0.3  0.4   0.5  3340    1
## q[111]       0.4       0 0.1  0.3  0.4   0.5  3340    1
## q[112]       0.4       0 0.1  0.3  0.4   0.5  3340    1
## q[113]       1.0       0 0.0  0.9  1.0   1.0  2572    1
## q[114]       0.4       0 0.1  0.3  0.4   0.5  3343    1
## q[115]       0.9       0 0.0  0.8  0.9   1.0  3451    1
## q[116]       0.7       0 0.1  0.6  0.7   0.9  3308    1
## q[117]       0.7       0 0.1  0.6  0.7   0.9  3308    1
## q[118]       0.9       0 0.0  0.8  0.9   1.0  3422    1
## q[119]       0.4       0 0.1  0.3  0.4   0.5  3346    1
## q[120]       1.0       0 0.0  0.9  1.0   1.0  2570    1
## q[121]       0.7       0 0.1  0.6  0.7   0.9  3320    1
## q[122]       0.4       0 0.1  0.3  0.4   0.5  3346    1
## q[123]       0.4       0 0.1  0.3  0.4   0.5  3347    1
## q[124]       0.7       0 0.1  0.6  0.7   0.9  3340    1
## q[125]       0.4       0 0.1  0.3  0.4   0.5  3347    1
## q[126]       0.4       0 0.1  0.3  0.4   0.5  3347    1
## q[127]       1.0       0 0.0  0.9  1.0   1.0  2569    1
## q[128]       0.4       0 0.1  0.3  0.4   0.5  3347    1
## q[129]       0.9       0 0.0  0.8  0.9   1.0  3390    1
## q[130]       0.7       0 0.1  0.6  0.7   0.8  3359    1
## q[131]       1.0       0 0.0  0.9  1.0   1.0  2567    1
## q[132]       0.4       0 0.1  0.3  0.4   0.5  3348    1
## q[133]       0.7       0 0.1  0.6  0.7   0.8  3359    1
## q[134]       1.0       0 0.0  0.9  1.0   1.0  2567    1
## q[135]       0.4       0 0.1  0.3  0.4   0.5  3348    1
## q[136]       0.9       0 0.0  0.8  0.9   1.0  3324    1
## q[137]       0.7       0 0.1  0.6  0.7   0.8  3377    1
## q[138]       1.0       0 0.0  0.9  1.0   1.0  2565    1
## q[139]       0.7       0 0.1  0.6  0.7   0.8  3377    1
## q[140]       0.7       0 0.1  0.6  0.7   0.8  3377    1
## q[141]       0.9       0 0.0  0.8  0.9   1.0  3324    1
## q[142]       0.4       0 0.1  0.2  0.3   0.5  3348    1
## q[143]       0.9       0 0.0  0.8  0.9   1.0  3288    1
## q[144]       0.3       0 0.1  0.2  0.3   0.5  3348    1
## q[145]       0.9       0 0.0  0.8  0.9   1.0  3288    1
## q[146]       0.9       0 0.0  0.8  0.9   1.0  3288    1
## q[147]       0.7       0 0.1  0.6  0.7   0.8  3394    1
## q[148]       0.9       0 0.0  0.8  0.9   1.0  3288    1
## q[149]       0.9       0 0.0  0.8  0.9   1.0  3252    1
## q[150]       0.7       0 0.1  0.6  0.7   0.8  3409    1
## q[151]       0.3       0 0.1  0.2  0.3   0.5  3309    1
## q[152]       0.7       0 0.1  0.5  0.7   0.8  3424    1
## q[153]       0.3       0 0.1  0.2  0.3   0.5  3309    1
## q[154]       0.9       0 0.1  0.8  0.9   1.0  3216    1
## q[155]       0.7       0 0.1  0.5  0.7   0.8  3424    1
## q[156]       0.7       0 0.1  0.5  0.7   0.8  3424    1
## q[157]       0.7       0 0.1  0.5  0.7   0.8  3437    1
## q[158]       0.3       0 0.1  0.2  0.3   0.5  3290    1
## q[159]       0.3       0 0.1  0.2  0.3   0.5  3290    1
## q[160]       0.7       0 0.1  0.5  0.7   0.8  3437    1
## q[161]       0.3       0 0.1  0.2  0.3   0.5  3272    1
## q[162]       0.9       0 0.1  0.8  0.9   1.0  3173    1
## q[163]       1.0       0 0.0  0.9  1.0   1.0  2546    1
## q[164]       0.3       0 0.1  0.2  0.3   0.5  3272    1
## q[165]       0.3       0 0.1  0.2  0.3   0.5  3272    1
## q[166]       0.7       0 0.1  0.5  0.7   0.8  3443    1
## q[167]       0.3       0 0.1  0.2  0.3   0.5  3259    1
## q[168]       0.9       0 0.1  0.8  0.9   1.0  3151    1
## q[169]       0.3       0 0.1  0.2  0.3   0.5  3251    1
## q[170]       0.9       0 0.1  0.8  0.9   1.0  3128    1
## q[171]       0.3       0 0.1  0.2  0.3   0.5  3251    1
## q[172]       0.9       0 0.1  0.8  0.9   1.0  3128    1
## q[173]       0.7       0 0.1  0.5  0.7   0.8  3449    1
## q[174]       0.3       0 0.1  0.2  0.3   0.4  3241    1
## q[175]       0.9       0 0.1  0.8  0.9   1.0  3105    1
## q[176]       1.0       0 0.0  0.9  1.0   1.0  2516    1
## q[177]       0.7       0 0.1  0.5  0.7   0.8  3449    1
## q[178]       0.3       0 0.1  0.2  0.3   0.4  3241    1
## q[179]       0.3       0 0.1  0.2  0.3   0.4  3241    1
## q[180]       0.7       0 0.1  0.5  0.7   0.8  3449    1
## q[181]       0.9       0 0.1  0.8  0.9   1.0  3083    1
## q[182]       1.0       0 0.0  0.9  1.0   1.0  2507    1
## q[183]       0.9       0 0.1  0.8  0.9   1.0  3083    1
## q[184]       0.6       0 0.1  0.5  0.6   0.8  3449    1
## q[185]       0.9       0 0.1  0.8  0.9   1.0  3083    1
## q[186]       0.6       0 0.1  0.5  0.6   0.8  3449    1
## q[187]       0.3       0 0.1  0.1  0.3   0.4  3221    1
## q[188]       0.3       0 0.1  0.1  0.3   0.4  3221    1
## q[189]       0.9       0 0.1  0.7  0.9   1.0  3060    1
## q[190]       0.3       0 0.1  0.1  0.3   0.4  3210    1
## q[191]       0.6       0 0.1  0.4  0.6   0.8  3443    1
## q[192]       0.3       0 0.1  0.1  0.3   0.4  3210    1
## q[193]       1.0       0 0.0  0.9  1.0   1.0  2489    1
## q[194]       0.6       0 0.1  0.4  0.6   0.8  3443    1
## q[195]       0.3       0 0.1  0.1  0.3   0.4  3210    1
## q[196]       0.9       0 0.1  0.7  0.9   1.0  3031    1
## q[197]       0.3       0 0.1  0.1  0.3   0.4  3199    1
## q[198]       0.3       0 0.1  0.1  0.3   0.4  3199    1
## q[199]       0.9       0 0.1  0.7  0.9   1.0  3031    1
## q[200]       0.6       0 0.1  0.4  0.6   0.8  3439    1
## lambda[1]    1.6       0 0.1  1.5  1.6   1.8  2895    1
## lambda[2]    1.1       0 0.1  0.9  1.1   1.2  3487    1
## lambda[3]    0.9       0 0.1  0.7  0.9   1.1  2818    1
## lambda[4]    1.8       0 0.1  1.7  1.8   2.0  2362    1
## lambda[5]    1.8       0 0.1  1.7  1.8   2.0  2362    1
## lambda[6]    1.1       0 0.1  0.9  1.1   1.3  3525    1
## lambda[7]    1.7       0 0.1  1.5  1.7   1.8  2942    1
## lambda[8]    1.8       0 0.1  1.7  1.8   2.0  2362    1
## lambda[9]    1.9       0 0.1  1.7  1.9   2.0  2363    1
## lambda[10]   1.9       0 0.1  1.7  1.9   2.0  2363    1
## lambda[11]   1.9       0 0.1  1.7  1.9   2.0  2363    1
## lambda[12]   1.7       0 0.1  1.5  1.7   1.9  2969    1
## lambda[13]   1.9       0 0.1  1.7  1.9   2.0  2363    1
## lambda[14]   1.1       0 0.1  1.0  1.1   1.3  3646    1
## lambda[15]   1.9       0 0.1  1.7  1.9   2.0  2364    1
## lambda[16]   1.1       0 0.1  1.0  1.1   1.3  3646    1
## lambda[17]   1.7       0 0.1  1.6  1.7   1.9  2999    1
## lambda[18]   1.9       0 0.1  1.7  1.9   2.0  2364    1
## lambda[19]   1.2       0 0.1  1.0  1.2   1.3  3687    1
## lambda[20]   1.2       0 0.1  1.0  1.2   1.3  3687    1
## lambda[21]   1.9       0 0.1  1.8  1.9   2.0  2367    1
## lambda[22]   1.2       0 0.1  1.0  1.2   1.3  3687    1
## lambda[23]   1.2       0 0.1  1.0  1.2   1.3  3687    1
## lambda[24]   1.2       0 0.1  1.0  1.2   1.3  3687    1
## lambda[25]   1.7       0 0.1  1.6  1.7   1.9  3031    1
## lambda[26]   1.9       0 0.1  1.8  1.9   2.0  2367    1
## lambda[27]   1.9       0 0.1  1.8  1.9   2.0  2367    1
## lambda[28]   1.9       0 0.1  1.8  1.9   2.1  2372    1
## lambda[29]   1.2       0 0.1  1.0  1.2   1.3  3729    1
## lambda[30]   1.0       0 0.1  0.8  1.0   1.2  2699    1
## lambda[31]   1.0       0 0.1  0.8  1.0   1.2  2675    1
## lambda[32]   1.9       0 0.1  1.8  1.9   2.1  2378    1
## lambda[33]   2.0       0 0.1  1.8  2.0   2.1  2386    1
## lambda[34]   1.8       0 0.1  1.7  1.8   1.9  3151    1
## lambda[35]   2.0       0 0.1  1.8  2.0   2.1  2386    1
## lambda[36]   2.0       0 0.1  1.8  2.0   2.1  2386    1
## lambda[37]   2.0       0 0.1  1.8  2.0   2.1  2386    1
## lambda[38]   1.2       0 0.1  1.1  1.2   1.4  3808    1
## lambda[39]   2.0       0 0.1  1.8  2.0   2.1  2386    1
## lambda[40]   1.2       0 0.1  1.1  1.2   1.4  3808    1
## lambda[41]   2.0       0 0.1  1.8  2.0   2.1  2386    1
## lambda[42]   1.2       0 0.1  1.1  1.2   1.4  3844    1
## lambda[43]   1.8       0 0.1  1.7  1.8   2.0  3199    1
## lambda[44]   1.2       0 0.1  1.1  1.2   1.4  3844    1
## lambda[45]   1.8       0 0.1  1.7  1.8   2.0  3199    1
## lambda[46]   1.1       0 0.1  0.9  1.1   1.3  2627    1
## lambda[47]   1.8       0 0.1  1.7  1.8   2.0  3251    1
## lambda[48]   1.3       0 0.1  1.1  1.3   1.4  3876    1
## lambda[49]   1.3       0 0.1  1.1  1.3   1.4  3876    1
## lambda[50]   1.3       0 0.1  1.1  1.3   1.4  3876    1
## lambda[51]   2.0       0 0.1  1.9  2.0   2.1  2409    1
## lambda[52]   2.0       0 0.1  1.9  2.0   2.1  2409    1
## lambda[53]   2.0       0 0.1  1.9  2.0   2.1  2409    1
## lambda[54]   2.0       0 0.1  1.9  2.0   2.1  2409    1
## lambda[55]   2.0       0 0.1  1.9  2.0   2.1  2424    1
## lambda[56]   1.1       0 0.1  0.9  1.1   1.3  2576    1
## lambda[57]   1.9       0 0.1  1.7  1.9   2.0  3308    1
## lambda[58]   1.9       0 0.1  1.7  1.9   2.0  3308    1
## lambda[59]   2.0       0 0.1  1.9  2.0   2.2  2432    1
## lambda[60]   1.3       0 0.1  1.2  1.3   1.4  3923    1
## lambda[61]   1.1       0 0.1  1.0  1.1   1.3  2550    1
## lambda[62]   1.3       0 0.1  1.2  1.3   1.4  3923    1
## lambda[63]   1.3       0 0.1  1.2  1.3   1.4  3923    1
## lambda[64]   1.3       0 0.1  1.2  1.3   1.4  3923    1
## lambda[65]   2.1       0 0.1  2.0  2.1   2.2  2439    1
## lambda[66]   2.1       0 0.1  2.0  2.1   2.2  2439    1
## lambda[67]   2.1       0 0.1  2.0  2.1   2.2  2439    1
## lambda[68]   2.1       0 0.1  2.0  2.1   2.2  2439    1
## lambda[69]   2.1       0 0.1  2.0  2.1   2.2  2447    1
## lambda[70]   1.9       0 0.1  1.8  1.9   2.0  3441    1
## lambda[71]   1.3       0 0.1  1.2  1.3   1.5  3940    1
## lambda[72]   1.9       0 0.1  1.8  1.9   2.0  3441    1
## lambda[73]   2.1       0 0.1  2.0  2.1   2.2  2447    1
## lambda[74]   1.9       0 0.1  1.8  1.9   2.0  3441    1
## lambda[75]   1.3       0 0.1  1.2  1.3   1.5  3940    1
## lambda[76]   1.3       0 0.1  1.2  1.3   1.5  3940    1
## lambda[77]   1.3       0 0.1  1.2  1.3   1.5  3940    1
## lambda[78]   2.1       0 0.1  2.0  2.1   2.2  2455    1
## lambda[79]   2.1       0 0.1  2.0  2.1   2.2  2455    1
## lambda[80]   1.2       0 0.1  1.0  1.2   1.4  2474    1
## lambda[81]   1.4       0 0.1  1.2  1.4   1.5  3935    1
## lambda[82]   1.9       0 0.1  1.8  1.9   2.1  3483    1
## lambda[83]   2.1       0 0.1  2.0  2.1   2.2  2455    1
## lambda[84]   2.1       0 0.1  2.0  2.1   2.2  2455    1
## lambda[85]   1.9       0 0.1  1.8  1.9   2.1  3483    1
## lambda[86]   2.1       0 0.1  2.0  2.1   2.2  2455    1
## lambda[87]   2.1       0 0.1  2.0  2.1   2.2  2467    1
## lambda[88]   2.1       0 0.1  2.0  2.1   2.2  2467    1
## lambda[89]   2.1       0 0.1  2.0  2.1   2.2  2467    1
## lambda[90]   1.2       0 0.1  1.0  1.2   1.4  2450    1
## lambda[91]   2.0       0 0.1  1.8  2.0   2.1  3551    1
## lambda[92]   2.1       0 0.1  2.0  2.1   2.3  2483    1
## lambda[93]   2.1       0 0.1  2.0  2.1   2.3  2483    1
## lambda[94]   1.2       0 0.1  1.1  1.2   1.4  2426    1
## lambda[95]   1.4       0 0.1  1.3  1.4   1.5  3893    1
## lambda[96]   1.2       0 0.1  1.1  1.2   1.4  2426    1
## lambda[97]   1.4       0 0.1  1.3  1.4   1.5  3893    1
## lambda[98]   2.2       0 0.1  2.1  2.2   2.3  2501    1
## lambda[99]   2.2       0 0.1  2.1  2.2   2.3  2501    1
## lambda[100]  1.3       0 0.1  1.1  1.3   1.4  2403    1
## lambda[101]  1.4       0 0.1  1.3  1.4   1.5  3857    1
## lambda[102]  2.2       0 0.1  2.1  2.2   2.3  2501    1
## lambda[103]  1.4       0 0.1  1.3  1.4   1.5  3857    1
## lambda[104]  1.4       0 0.1  1.3  1.4   1.6  3813    1
## lambda[105]  1.4       0 0.1  1.3  1.4   1.6  3813    1
## lambda[106]  2.2       0 0.1  2.1  2.2   2.3  2523    1
## lambda[107]  2.2       0 0.1  2.1  2.2   2.3  2554    1
## lambda[108]  2.2       0 0.1  2.1  2.2   2.3  2554    1
## lambda[109]  2.0       0 0.1  1.9  2.0   2.1  3871    1
## lambda[110]  2.2       0 0.1  2.1  2.2   2.3  2554    1
## lambda[111]  2.2       0 0.1  2.1  2.2   2.3  2554    1
## lambda[112]  2.2       0 0.1  2.1  2.2   2.3  2554    1
## lambda[113]  1.3       0 0.1  1.1  1.3   1.5  2343    1
## lambda[114]  2.2       0 0.1  2.1  2.2   2.3  2589    1
## lambda[115]  2.1       0 0.1  1.9  2.1   2.2  3920    1
## lambda[116]  1.5       0 0.1  1.3  1.5   1.6  3701    1
## lambda[117]  1.5       0 0.1  1.3  1.5   1.6  3701    1
## lambda[118]  2.1       0 0.1  2.0  2.1   2.2  3960    1
## lambda[119]  2.2       0 0.1  2.1  2.2   2.4  2628    1
## lambda[120]  1.3       0 0.1  1.2  1.3   1.5  2325    1
## lambda[121]  1.5       0 0.1  1.4  1.5   1.6  3637    1
## lambda[122]  2.2       0 0.1  2.1  2.2   2.4  2628    1
## lambda[123]  2.3       0 0.1  2.1  2.3   2.4  2669    1
## lambda[124]  1.5       0 0.1  1.4  1.5   1.6  3570    1
## lambda[125]  2.3       0 0.1  2.1  2.3   2.4  2669    1
## lambda[126]  2.3       0 0.1  2.1  2.3   2.4  2669    1
## lambda[127]  1.4       0 0.1  1.2  1.4   1.5  2310    1
## lambda[128]  2.3       0 0.1  2.1  2.3   2.4  2669    1
## lambda[129]  2.1       0 0.1  2.0  2.1   2.2  3987    1
## lambda[130]  1.5       0 0.1  1.4  1.5   1.7  3502    1
## lambda[131]  1.4       0 0.1  1.2  1.4   1.5  2297    1
## lambda[132]  2.3       0 0.1  2.2  2.3   2.4  2709    1
## lambda[133]  1.5       0 0.1  1.4  1.5   1.7  3502    1
## lambda[134]  1.4       0 0.1  1.2  1.4   1.5  2297    1
## lambda[135]  2.3       0 0.1  2.2  2.3   2.4  2709    1
## lambda[136]  2.1       0 0.1  2.0  2.1   2.2  3969    1
## lambda[137]  1.6       0 0.1  1.4  1.6   1.7  3434    1
## lambda[138]  1.4       0 0.1  1.2  1.4   1.6  2285    1
## lambda[139]  1.6       0 0.1  1.4  1.6   1.7  3434    1
## lambda[140]  1.6       0 0.1  1.4  1.6   1.7  3434    1
## lambda[141]  2.1       0 0.1  2.0  2.1   2.2  3969    1
## lambda[142]  2.3       0 0.1  2.2  2.3   2.4  2751    1
## lambda[143]  2.2       0 0.1  2.0  2.2   2.3  3937    1
## lambda[144]  2.3       0 0.1  2.2  2.3   2.4  2792    1
## lambda[145]  2.2       0 0.1  2.0  2.2   2.3  3937    1
## lambda[146]  2.2       0 0.1  2.0  2.2   2.3  3937    1
## lambda[147]  1.6       0 0.1  1.4  1.6   1.7  3367    1
## lambda[148]  2.2       0 0.1  2.0  2.2   2.3  3937    1
## lambda[149]  2.2       0 0.1  2.1  2.2   2.3  3893    1
## lambda[150]  1.6       0 0.1  1.5  1.6   1.7  3303    1
## lambda[151]  2.4       0 0.1  2.2  2.4   2.5  2881    1
## lambda[152]  1.6       0 0.1  1.5  1.6   1.8  3241    1
## lambda[153]  2.4       0 0.1  2.2  2.4   2.5  2881    1
## lambda[154]  2.2       0 0.1  2.1  2.2   2.3  3840    1
## lambda[155]  1.6       0 0.1  1.5  1.6   1.8  3241    1
## lambda[156]  1.6       0 0.1  1.5  1.6   1.8  3241    1
## lambda[157]  1.6       0 0.1  1.5  1.6   1.8  3183    1
## lambda[158]  2.4       0 0.1  2.3  2.4   2.5  2894    1
## lambda[159]  2.4       0 0.1  2.3  2.4   2.5  2894    1
## lambda[160]  1.6       0 0.1  1.5  1.6   1.8  3183    1
## lambda[161]  2.4       0 0.1  2.3  2.4   2.5  2904    1
## lambda[162]  2.2       0 0.1  2.1  2.2   2.4  3718    1
## lambda[163]  1.5       0 0.1  1.3  1.5   1.7  2260    1
## lambda[164]  2.4       0 0.1  2.3  2.4   2.5  2904    1
## lambda[165]  2.4       0 0.1  2.3  2.4   2.5  2904    1
## lambda[166]  1.7       0 0.1  1.5  1.7   1.8  3128    1
## lambda[167]  2.4       0 0.1  2.3  2.4   2.6  2911    1
## lambda[168]  2.3       0 0.1  2.1  2.3   2.4  3654    1
## lambda[169]  2.4       0 0.1  2.3  2.4   2.6  2916    1
## lambda[170]  2.3       0 0.1  2.2  2.3   2.4  3591    1
## lambda[171]  2.4       0 0.1  2.3  2.4   2.6  2916    1
## lambda[172]  2.3       0 0.1  2.2  2.3   2.4  3591    1
## lambda[173]  1.7       0 0.1  1.6  1.7   1.9  2988    1
## lambda[174]  2.5       0 0.1  2.3  2.5   2.6  2919    1
## lambda[175]  2.3       0 0.1  2.2  2.3   2.4  3529    1
## lambda[176]  1.6       0 0.1  1.4  1.6   1.7  2264    1
## lambda[177]  1.7       0 0.1  1.6  1.7   1.9  2988    1
## lambda[178]  2.5       0 0.1  2.3  2.5   2.6  2919    1
## lambda[179]  2.5       0 0.1  2.3  2.5   2.6  2919    1
## lambda[180]  1.7       0 0.1  1.6  1.7   1.9  2988    1
## lambda[181]  2.3       0 0.1  2.2  2.3   2.4  3469    1
## lambda[182]  1.6       0 0.1  1.4  1.6   1.7  2256    1
## lambda[183]  2.3       0 0.1  2.2  2.3   2.4  3469    1
## lambda[184]  1.7       0 0.1  1.6  1.7   1.9  2948    1
## lambda[185]  2.3       0 0.1  2.2  2.3   2.4  3469    1
## lambda[186]  1.7       0 0.1  1.6  1.7   1.9  2948    1
## lambda[187]  2.5       0 0.1  2.3  2.5   2.7  2921    1
## lambda[188]  2.5       0 0.1  2.3  2.5   2.7  2921    1
## lambda[189]  2.3       0 0.1  2.2  2.3   2.5  3413    1
## lambda[190]  2.5       0 0.1  2.4  2.5   2.7  2920    1
## lambda[191]  1.8       0 0.1  1.6  1.8   1.9  2880    1
## lambda[192]  2.5       0 0.1  2.4  2.5   2.7  2920    1
## lambda[193]  1.6       0 0.1  1.4  1.6   1.8  2243    1
## lambda[194]  1.8       0 0.1  1.6  1.8   1.9  2880    1
## lambda[195]  2.5       0 0.1  2.4  2.5   2.7  2920    1
## lambda[196]  2.4       0 0.1  2.2  2.4   2.5  3310    1
## lambda[197]  2.5       0 0.1  2.4  2.5   2.7  2918    1
## lambda[198]  2.5       0 0.1  2.4  2.5   2.7  2918    1
## lambda[199]  2.4       0 0.1  2.2  2.4   2.5  3310    1
## lambda[200]  1.8       0 0.1  1.6  1.8   2.0  2850    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jan  5 14:14:40 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

\(q, \lambda\)の順位相関係数の中央値は負であるため、来店確率が高いことと、リピーターで多数来店することは、負の関係がある。

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)
r <- sapply(1:N_mcmc,
            function(i) cor(ms$lambda[i,], ms$q[i,], method = 'spearman'))
quantile(r, prob = c(0.025, 0.25, 0.5, 0.75, 0.975))
##       2.5%        25%        50%        75%      97.5% 
## -0.8081299 -0.6974237 -0.6496459 -0.6050023 -0.4827052