２種類の広告介入の効果を分解してみる

13877 ワード
状態空間モデル 時系列解析 Stan 時系列解析 テキストリンク

１つだとつまらないので、2つのtrendで

１つはt=30から(trend)、もう一つは80から(trend2)。

介入時の情報は広告の目的ならわかるので、trendの一回微分（速度）に対してexp(b*is_start)と入れた。
介入時点は微分不可だが、可能としても別にいい気もしたけど収束しないので・・・

観測変数のsigmaの制限と、事前分布は外したかったけどこれも収束しないのでやむなく。

library(tidyverse)

T=100

z1=c(rep(0,30),seq(0,20,length.out = 70))+2*rnorm(T) #広告1
z2=c(rep(0,80),seq(0,100,length.out = 20))+3*rnorm(T) #広告2

y=z1+z2+3*rnorm(T)

ts.plot(y)

z1_begin=rep(0,T)
z1_begin[30]=1
z2_begin=rep(0,T)
z2_begin[80]=1

df1=data.frame(y=y,T=T,z1,z2)

'''
![a7dc6424-ccb1-45d2-927d-30df92f3e46e.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38159/7fdf03fb-9282-34f3-e56e-1a4576be792a.png)

stanコードは以下
```r

scode="
data{
int N;
real y[N];
int arn;
int z1_begin[N];
int z2_begin[N];
}

parameters {
real<upper=15> sigma;
real trend[N];
real<upper=2> s_trend;

real trend2[N];
real<upper=10> s_trend2;

real b1;
real b2;

#real ar[N];
#real c_ar[arn];
#real s_ar;

}

model{
real ssum;
s_trend ~ normal(1,1);
#s_ar ~ normal(1, 1);
trend[1:30]~normal(0,2);
trend2[1:80]~normal(0,2);

for(t in 3:N){
trend[t]~normal(trend[t-1] -exp(b1*z1_begin[t-1])*(trend[t-1]-trend[t-2]),s_trend);
trend2[t]~normal(trend2[t-1] -exp(b2*z2_begin[t-1])*(trend2[t-1]-trend2[t-2]),s_trend2);


}


for(t in 1:N){
    y[t]~normal(trend[t]+trend2[t],sigma);
    #y[t]~normal(trend[t]+ar[t],sigma);
}

}"



library(rstan)
S=1
y=y
N=length(y)
time=1:N
sdat=list(y=y,N=N,z1_begin=z1_begin,z2_begin=z2_begin,
          arn=1)

t2i=rep(0,N)
t2i[80:100]=seq(0,100,length.out = 21)

fit=stan(model_code = scode,data=sdat,
         init=function() { list(#trend=rep(mean(y),N),
                                #trend2=t2i,

                                ar=rep(0,N)
         )},warmup=1000,iter=2000,chains=1)

trend=get_posterior_mean(fit,par="trend")
trend2=get_posterior_mean(fit,par="trend2")

sigma=get_posterior_mean(fit,par="sigma")
s_trend=get_posterior_mean(fit,par="s_trend")
ar=get_posterior_mean(fit,par="ar")

b=get_posterior_mean(fit,par="b1")
b2=get_posterior_mean(fit,par="b2")

print(sigma)

plot.ts=function(df_plt){
  df_plt=df_plt %>% gather(names(df_plt)[names(df_plt)!="time"],key=variable,value=value)
  df_plt  %>% 
    ggplot(aes(x=time,y=value,colour=variable))+geom_line()
}

df_plt=data_frame(y,trend=trend,trend2=trend2,time)

print(plot.ts(df_plt))

上手くいったら、介入時刻の情報無しでやってみます