跟着开源项目学因果推断——CausalImpact 贝叶斯结构时间序列模型（二十一）

df=pd.read_csv("air_visit_data.csv",parse_dates=["visit_date"])

#元のデータは欠損があるため日時のベースを作成
basedate=pd.DataFrame(list(pd.date_range('2016-05-01', '2016-08-20')))
basedate.columns=["visit_date"]

#店舗を指定、日時を限定
store1=df[df["air_store_id"]=="air_e55abd740f93ecc4"]
store1=store1[(store1["visit_date"]>="2016-05") & (store1["visit_date"]<="2016-08-20")]
store2=df[df["air_store_id"]=="air_1653a6c513865af3"]
store2=store2[(store2["visit_date"]>="2016-05") & (store2["visit_date"]<="2016-08-20")]

#データを結合
store1.rename(columns={"visitors":"visitors_store1"},inplace=True)
store2.rename(columns={"visitors":"visitors_store2"},inplace=True)
data=pd.merge(basedate,store1[["visit_date","visitors_store1"]],on="visit_date",how="left")
data=pd.merge(data,store2[["visit_date","visitors_store2"]],on="visit_date",how="left")
data=data.interpolate() #欠損値を線形補間
data=data.set_index("visit_date")

from causalimpact import CausalImpact
pre_period = ['2016-05-01', '2016-08-10']
post_period = ['2016-08-11', '2016-08-20']
ci = CausalImpact(data, pre_period, post_period,nseasons=[{'period': 7}],trend=True)
print(ci.summary())

ci.trained_model.summary()

pre_period = ['2016-05-01', '2016-08-10']
post_period = ['2016-08-11', '2016-08-20']
ci = CausalImpact(data.iloc[:,0], pre_period, post_period,nseasons=[{'period': 7}])
ci.plot(figsize=(12, 6))

# 创建数据
set.seed(1)
x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 100)
y <- 1.2 * x1 + rnorm(100)
y[71:100] <- y[71:100] + 10
data <- cbind(y, x1)


pre.period <- c(1, 70)
post.period <- c(71, 100)
# 画图
impact <- CausalImpact(data, pre.period, post.period)
plot(impact)

%matplotlib inline

import sys
import os


sys.path.append(os.path.abspath('../'))


import numpy as np
import pandas as pd
import tensorflow as tf
import tensorflow_probability as tfp
import matplotlib.pyplot as plt
from causalimpact import CausalImpact


tfd = tfp.distributions
plt.rcParams['figure.figsize'] = [15, 10]

# 生成数据
# This is an example presented in Google's R code.
# Uses TensorFlow Probability to simulate a random walk process.
observed_stddev, observed_initial = (tf.convert_to_tensor(value=1, dtype=tf.float32),
                                     tf.convert_to_tensor(value=0., dtype=tf.float32))
level_scale_prior = tfd.LogNormal(loc=tf.math.log(0.05 * observed_stddev), scale=1, name='level_scale_prior')
initial_state_prior = tfd.MultivariateNormalDiag(loc=observed_initial[..., tf.newaxis], scale_diag=(tf.abs(observed_initial) + observed_stddev)[..., tf.newaxis], name='initial_level_prior')
ll_ssm = tfp.sts.LocalLevelStateSpaceModel(100, initial_state_prior=initial_state_prior, level_scale=level_scale_prior.sample())
ll_ssm_sample = np.squeeze(ll_ssm.sample().numpy())

x0 = 100 * np.random.rand(100)
x1 = 90 * np.random.rand(100)
y = 1.2 * x0 + 0.9 * x1 + ll_ssm_sample
y[70:] += 10
data = pd.DataFrame({'x0': x0, 'x1': x1, 'y': y}, columns=['y', 'x0', 'x1'])

data.plot()
plt.axvline(69, linestyle='--', color='k')
plt.legend();

pre_period = [0, 69]
post_period = [70, 99]
ci = CausalImpact(data, pre_period, post_period)

ci.plot()

ci.plot(panels=['original', 'pointwise'], figsize=(15, 12))

print(ci.summary())

Posterior Inference {Causal Impact}
                          Average            Cumulative
Actual                    106.01             3180.15
Prediction (s.d.)         96.16 (3.6)        2884.67 (107.86)
95% CI                    [89.16, 103.25]    [2674.67, 3097.47]

Absolute effect (s.d.)    9.85 (3.6)         295.49 (107.86)
95% CI                    [2.76, 16.85]      [82.69, 505.48]

Relative effect (s.d.)    10.24% (3.74%)     10.24% (3.74%)
95% CI                    [2.87%, 17.52%]    [2.87%, 17.52%]

Posterior tail-area probability p: 0.0
Posterior prob. of a causal effect: 99.6%

For more details run the command: print(impact.summary('report'))

print(ci.summary(output='report'))

Analysis report {CausalImpact}


During the post-intervention period, the response variable had
an average value of approx. 106.01. By contrast, in the absence of an
intervention, we would have expected an average response of 96.16.
The 95% interval of this counterfactual prediction is [89.16, 103.25].
Subtracting this prediction from the observed response yields
an estimate of the causal effect the intervention had on the
response variable. This effect is 9.85 with a 95% interval of
[2.76, 16.85]. For a discussion of the significance of this effect,
see below.

ci.inferences.head()

print('Mean value of noise observation std: ', ci.model_samples[0].numpy().mean())
print('Mean value of level std: ', ci.model_samples[1].numpy().mean())
print('Mean value of linear regression x0, x1: ', ci.model_samples[2].numpy().mean(axis=0))

ci = CausalImpact(data, pre_period, post_period, model_args={'prior_level_sd':0.1})
ci.plot(figsize=(15, 12))

observed_stddev, observed_initial = (5., 0.)
num_seasons = 7

drift_scale_prior = tfd.LogNormal(loc=tf.math.log(.01 * observed_stddev), scale=1.)
initial_effect_prior = tfd.Normal(loc=observed_initial, scale=tf.abs(observed_initial) + observed_stddev)
initial_state_prior = tfd.MultivariateNormalDiag(loc=tf.stack([initial_effect_prior.mean()] * num_seasons, axis=-1),
                                                 scale_diag=tf.stack([initial_effect_prior.stddev()] * num_seasons, axis=-1))

s_ssm = tfp.sts.SeasonalStataeSpaceModel(
    num_timesteps=100,
    num_seasons=num_seasons,
    num_steps_per_season=1,
    drift_scale=drift_scale_prior.sample(),
    initial_state_prior=initial_state_prior
)

plt.plot(s_ssm.sample());

season_data = data.copy()
season_data['y'] += np.squeeze(s_ssm.sample().numpy())

ci = CausalImpact(season_data, pre_period, post_period, model_args={'nseasons': 7})

print(ci.summary())

ci.plot()

Posterior Inference {Causal Impact}
                          Average            Cumulative
Actual                    103.22             3096.61
Prediction (s.d.)         91.52 (3.48)       2745.52 (104.37)
95% CI                    [84.63, 98.26]     [2538.75, 2947.89]

Absolute effect (s.d.)    11.7 (3.48)        351.09 (104.37)
95% CI                    [4.96, 18.6]       [148.72, 557.85]

Relative effect (s.d.)    12.79% (3.8%)      12.79% (3.8%)
95% CI                    [5.42%, 20.32%]    [5.42%, 20.32%]

Posterior tail-area probability p: 0.0
Posterior prob. of a causal effect: 99.9%

For more details run the command: print(impact.summary('report'))