Introduction

第16局の総則6.10「溶出試験法」に装置、測定手法、「判定法」が規定。
判定法には判定法１と判定法２があり、前者は3極協調されている。後者は日本特有の規定。判定法２について、Stanを用いて推定し、局方の設定を理解してみる。
JAGSでもなんでもいいんだが、generated quantitiesがとてつもなく便利なのでStan利用。

溶出試験の判定法

判定法1
- 判定基準表を基に判定
- 表の水準S1、S2を満たさない場合はS3まで実施
判定法2
- 試料6個について試験し、個々の試料の溶出率全てが規格に規定する値の時
- 規格する値から外れた試料が1または2個の時、新たに試料6個を試験
  - 12個中、10個以上の試料の個々の溶出率が規定する値の時、適合とする

表：判定法１

	試験個数	判定基準
S1	6	個々の試料からの溶出が$Q+5%$以上
S2	6	12個(s1+s2)の試料の平均溶出率$\ge Q$
		且つ$Q - 15\%$のモノがない
S3	12	24個(S1+S2+S3)の試料の平均溶出率$\ge Q$
		$Q-15\% $未満のモノが2個以下
		$Q-25\% $未満のモノがない

推定用コード

頻度論では①に母数を点と見做すこと、②に大数の法則前提、（多くの場合、母数をサンプルから推定するとして）、③標準偏差の推定が初等以上の知識がいる…など利用が難しい。④信頼性区間を用いる区間推定が実は困難。。。。
など、高度な数学を習っていない薬剤師の私では計算困難。だがベイズを使えばとても簡単に推定できるので、その解釈から局方の設定原理を理解しようという雑文。

ドメイン知識を平均に対して用いる場合は下記になるだろう。事前分布のsigmaには無条件分布を設定。

平均母数muに一様分布やStanお任せ設定にすれば無条件情報になるだろう。

Domain.NormalDistModel

// DOMAIN INFORMATION -3*sd0 < mu0 < +3*sd0
data {
  int<lower=0> N;
  vector[N]    Y;
  real LowerPoint;
  real mu0;
  real sd0;
}
parameters {
  real mu;
  real sigma;
}
transformed parameters {
}
model {
// mu can be changed to uniform distribution(Non-informative prior)
     mu     ~ normal(100, sd0);
     sigma  ~ gamma(1.0E-3, 1.0E-3);

// posterior distribution vector
  Y ~ normal(mu, sigma);
}
generated quantities{
  real Q_sample;
  real Q_normal;
  real RR;
  Q_sample = normal_cdf(LowerPoint, mu, sigma);
  Q_normal = step(LowerPoint - normal_rng(mu0, sd0));
  RR = normal_cdf(LowerPoint, mu, sigma)/ normal_cdf(LowerPoint, mu0, sd0);;

規格は100% ± 5と仮定し、その中の最低の測定値に$3\sigma$以下のモノが局方仮定のうち、一つだけ見つかった時、規格とくらべてどんな差があるかを確認する。

set.seed(3141592)
s0 <- rnorm(12,100,5)
s1 <- s0[1:6]
s1[s1==min(s1)] <- 84.9
s2 <- c(s1:s0[7:12])

仮想データを設定し、下記のように推定する。規格外が１つか2つ発覚した場合、測定法２の２までで12サイズの乱数が必用。なんでいっぺんに作成しておき最初の6個を用いた。

model.stan <- stan_model(model_code = Domain.NormalDistModel)

fit6 <- sampling(object = model.stan, pars = parameters,
      data = list(N = length(s1), Y = s1, mu0=100, sd0=5, LowerPoint = 85),
      warmup = 10000, iter = 60000,
      chain  = 4, seed = 123, thin =1,
      control = list(adapt_delta = 0.99),
      show_messages = F, open_progress =F)

4 chains, each with iter=60000; warmup=10000; thin=2; 
post-warmup draws per chain=25000, total post-warmup draws=100000.

           mean se_mean    sd   2.5%    25%    50%    75%  97.5%  n_eff Rhat
mu        98.30    0.01  2.65  93.21  96.58  98.22  99.93 103.80  63023    1
sigma      7.74    0.02  2.96   4.18   5.77   7.08   8.91  15.27  35492    1
Q_sample   0.05    0.00  0.06   0.00   0.01   0.03   0.08   0.23  43099    1
Q_normal   0.00    0.00  0.04   0.00   0.00   0.00   0.00   0.00 100000    1
RR        39.86    0.22 45.94   0.37   7.33  22.86  55.94 169.42  43099    1
lp__     -17.12    0.01  1.28 -20.57 -17.63 -16.73 -16.19 -15.84  33857    1

理論的に$3\sigma$を下回る確率はpnorm(-3)と計算できるが、Q_normalとして図示したとき便利にしてみた。

事後平均（EAP)推定で39.9倍、中央値で22.8倍ですか。高いね。

背景に複雑な数理が（たぶん）潜んでいるんだろうが、有意性検定のアナロジーの元で中央値も平均も
20倍以上だから、有意水準5%有意になるのだろう。

宿題は、
- generated quantitiesパートでリスク差も計算（いれとけ俺）
- 6個中1の規格外があれば全部で12個で評価する部分のシミュレーション
- 下限値の値をいろいろ変えてみたら結果はどうなるか？的な
- 図示化
- わすれちゃいけない判定法１のモデル（やっとけ俺）

うんうんなるほど。勉強になりました。

局方（JP16)における溶出試験の理解

Introduction

溶出試験の判定法

推定用コード