a*bの値が与えられたときにaとbをベイズで推定してみる

基礎代謝から体重の推移のモデリングをしていたときに未知数が多く、どのくらい識別できるのか不安になったので、実験として。

データは単純にa,bを正規分布から生成

library(rstan)

N=1000
a=rnorm(N,mean=100,sd=10)
b=rnorm(N,mean=2,sd=0.2)

lst=list(N=N,y=a*b)
fit=stan(model_code = scode,data=lst,
         warmup=300,iter=1000,chain=1)

data {
  int N;
  real y[N];
}
parameters {
  real<lower=0> a;
  real b;
  real sigma;
}

model {
  sigma~cauchy(0.01,0.01);
  //a~normal(140,30); 
  //b~normal(2.5,0.5);
  a~normal(130,30);
  b~normal(1.5,0.5);
  
  y~normal(a*b,sigma);
}

print(fit)

事前分布を色々かえると、a,b同じ方向にずらしたときは比較的うまくいき、上下反対方向にずらすと的外れになっていくっぽい。ただし一応正しい値のほうにズレていく。

なぜ識別できるのか考えてみると、y,bが与えられたときのaの事後分布見ると
$P(a|y,b)∝P(y,b|a)*P(a)$
となるが、
$P(y,b|a)$ の部分は

yが平均付近でbが（真の）平均付近ならaも平均付近の確率が高く
yが平均付近でbが高ければaは低い確率が高く
それにaの事前分布をかけた値が高い所からサンプリングされる確率が高くなる。
どっちにしろ真の確率分布*事前分布が高い所からサンプリングされる確率が高くなり、収束していくのもわかる気がする。