廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

【rstan】グループ化時間&右切断ワイブル分布+ポアソンノイズで需要予測モデル

なにかの製品の月ごとなり週ごとなりの出荷数のデータがあるとします。
このデータを「製品が発売されてから消費者が購入に至るまでの待ち時間」を計測したものだと捉えなおすと、ワイブル分布を仮定して分析するのもさほど不自然ではないように思えます。

そう考えると製品の出荷数をどこかの時点までで切り取って観測していることになり、これは右切断データです。

月ごとなり週ごとなりで集計しているのでグループ化時間データでもあります。

時間 t_i (i=1,\ldots,T)での出荷数を y_i とすると尤度は以下のようになります。

\displaystyle \prod_{i=2}^{T} \left(\frac{F(t_{i})-F(t_{i-1})}{F(t_T)}\right)^{y_{i-1}}

R でシミュレーションしてみます。

set.seed(1)
x <- sort(rweibull(1000,2,100))
tim <-seq(0,150,by=10)
len <-length(tim)
y <-sapply(tim[-1],function(t)sum(t-10<=x & x<t))
plot(tim[-1],y,type="b",ylim=c(0,max(y)),xlab="time",ylab="sales")

f:id:abrahamcow:20171013001411p:plain

なんとなくこんな感じの時系列データを見たことがあるような気がします。

上述の尤度関数で最尤推定してみます。

ll <- function(par){
  m <- exp(par[1])
  eta <- exp(par[2])
sum(y*(log(pweibull(tim[-1],m,eta)/pweibull(tim[len],m,eta)-
             pweibull(tim[-len],m,eta)/pweibull(tim[len],m,eta))))
}
opt1 <-optim(c(1,4),ll,control = list(fnscale=-1))
probfun <-function(par){
  m <- par[1]
  eta <- par[2]
  (pweibull(tim[-1],m,eta)-pweibull(tim[-len],m,eta))/
    pweibull(tim[len],m,eta)
}
plot(tim[-1],y,type="b",ylim=c(0,max(y)),xlab="time",ylab="sales")
lines(tim[-1],sum(y)*probfun(exp(opt1$par)),col="royalblue",lwd=2)

f:id:abrahamcow:20171013001544p:plain

うまく推定できました。

製品の総需要の合計は以下のように求まります。

> sum(y)/pweibull(tim[len],exp(opt1$par[1]),exp(opt1$par[2]))
[1] 1022.3

さて、実際のデータ z_t は製品の本来の需要 y_t に加えて、売れそうだからいっぱい仕入れとこうとか返品しようとか、そういったノイズがのっていることが想定されます。

とりあえずポアソン分布でノイズを加えてみます。

\displaystyle z_{t} \sim {\rm Poisson}(y_t)

y2 <-rpois(len-1,y)

z_t から y_t およびワイブル分布のパラメータを推定する問題を考えます。

y_t が離散変数だと微分できなくて推定が大変になるので、正規分布で近似することにします。

\displaystyle p_i = \left(\frac{F(t_{i})-F(t_{i-1})}{F(t_T)}\right)

として、

\displaystyle y_{t} \sim {\rm Normal}\left(p_t \sum y_t,p_t(1-p_t) \sum y_t\right)

と近似しました。

本当は共分散も考えたほうがいいのかな、と思ったけどまあいいことにします。

Stan のコードは以下のようになりました。

data{
  int<lower=1> Tim;
  int<lower=0> Y[Tim-1];
  real<lower=0> time[Tim];
}
parameters{
  real<lower=0> mu[Tim-1];
  real<lower=0> m;
  real<lower=0> eta;
}
transformed parameters{
  real summu;
  vector<lower=0, upper=1>[Tim-1] prob;
  summu=sum(mu);
  prob[1] = exp(weibull_lcdf(time[2]|m,eta)-weibull_lcdf(time[Tim]|m,eta));
  for(i in 3:Tim){
    prob[i-1] = exp(weibull_lcdf(time[i]|m,eta)-weibull_lcdf(time[Tim]|m,eta))-exp(weibull_lcdf(time[i-1]|m,eta)-weibull_lcdf(time[Tim]|m,eta));
  }
}
model{
  Y ~ poisson(mu);
  for(i in 1:(Tim-1)){
    mu[i] ~ normal(prob[i]*summu,sqrt(prob[i]*(1-prob[i])*summu)); 
  }
}
generated quantities{
  real<lower=0> Ym;
  Ym = summu/exp(weibull_lcdf(time[Tim]|m,eta));
}

いい感じに推定できました。

f:id:abrahamcow:20171013003305p:plain

R のコードは以下です。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
weibull_noiz <-stan_model("~/Documents/weibull_noiz.stan")
dat <- list(Tim=length(tim),Y=y2,time=tim)
fit1 <- sampling(weibull_noiz,dat)

traceplot(fit1,pars=c("m","eta","Ym"))
muhat <-get_posterior_mean(fit1,pars="mu")[,"mean-all chains"]
plot(tim[-1],y2,type="b",ylim=c(0,max(y2)),xlab="time",ylab="sales")
lines(tim[-1],muhat,col="royalblue",lwd=2)

f:id:abrahamcow:20171013003432p:plain

以上です。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)