这是一篇超级超级长的文章，今天我们将继续分析这个案例研究。
前文回顾：(PyStan)零售价格贝叶斯策略建模（上）
观察：
按类别划分的价格有相当多的变化，我们可以看到，例如，类别美容/香水/女性（指数为9），有大量样本（225），也有一个最严格的估计值范围。
尽管美容/护发/洗发水加护发素（指数16）的样品数量最少（仅一个），但也有一个最广泛的估计范围。
我们可以可视化和β的参数估计分布。
代码语言：javascript复制az.plot_trace(varying_intercept_fit, var_names = ['sigma_a', 'b']);
varying_intercept_fit['b'].mean()
装运系数的估计值约为-0.27，这可以解释为卖方支付的产品装运费约为买方支付的装运价格（exp（-0.27）=0.76）的0.76倍（在考虑类别后）。
代码语言：javascript复制plt.figure(figsize=(12, 6))
xvals = np.arange(2)
bp = varying_intercept_fit['a'].mean(axis=0) # mean a (intercept) by category
mp = varying_intercept_fit['b'].mean() # mean b (slope/shipping effect)
for bi in bp:
plt.plot(xvals, mp*xvals bi, 'bo-', alpha=0.4)
plt.xlim(-0.1,1.1)
plt.xticks([0, 1])
plt.title('Fitted relationships by category')
plt.xlabel("shipping")
plt.ylabel("log(price)");
观察：
从这个图中可以清楚地看出，我们为每个类别设置了相同的运输效果，但每个类别的价格水平不同。
有一个类别的拟合价格估计值非常低，有几个类别的拟合价格估计值相对较低。
有多个类别的价格估计相对较高。
大部分类别构成了一组类似的拟合。
我们可以看到，对于样本量较小的类别，类别级价格估计的部分集中是否比集中或未集中模型提供了更合理的估计。
Partial Pooling -坡度变化模型
我们还可以建立一个模型，允许类别根据运输安排（由买方支付或由卖方支付）而变化，从而影响价格。方程如下：
变坡模型：
代码语言：javascript复制varying_slope = """
data {
int<lower=0> J;
int<lower=0> N;
int<lower=1,upper=J> category[N];
vector[N] x;
vector[N] y;
}
parameters {
real a;
vector[J] b;
real mu_b;
real<lower=0,upper=100> sigma_b;
real<lower=0,upper=100> sigma_y;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- a x[i] * b[category[i]];
}
model {
sigma_b ~ uniform(0, 100);
b ~ normal (mu_b, sigma_b);
a ~ normal (0, 1);
sigma_y ~ uniform(0, 100);
y ~ normal(y_hat, sigma_y);
}
"""
拟合模型：
代码语言：javascript复制varying_slope_data = {'N': len(log_price),
'J': len(n_category),
'category': category 1, # Stan counts starting at 1
'x': shipping,
'y': log_price}
sm = pystan.StanModel(model_code=varying_slope)
varying_slope_fit = sm.sampling(data=varying_slope_data, iter=1000, chains=2)
按照前面的过程，我们将可视化20个类别。
代码语言：javascript复制b_sample = pd.DataFrame(varying_slope_fit['b'])
plt.figure(figsize=(20, 5))
g = sns.boxplot(data=b_sample.iloc[:,0:20], whis=np.inf, color="g")
# g.set_xticklabels(df.category_name.unique(), rotation=90) # label counties
g.set_title("Estimate of shipping effect, by category")
g;
观察：
乍一看，我们可能看不到这两个箱线图之间的任何区别。但是如果你看得更深入，你会发现在不同的坡度模型中，不同类别的中间估计值的变化比在不同的截距模型中的变化要小，尽管不确定性的范围仍然是最大的类别与最少的产品,并至少在最多的类别产品中会如此。
可视化拟合模型：
代码语言：javascript复制plt.figure(figsize=(10, 6))
xvals = np.arange(2)
b = varying_slope_fit['a'].mean()
m = varying_slope_fit['b'].mean(axis=0)
for mi in m:
plt.plot(xvals, mi*xvals b, 'bo-', alpha=0.4)
plt.xlim(-0.2, 1.2)
plt.xticks([0, 1])
plt.title("Fitted relationships by category")
plt.xlabel("shipping")
plt.ylabel("log(price)");
观察：
从这个图中可以清楚地看出，我们已经为每个类别安装了相同的价格水平，但在每个类别中具有不同的运输效果。
有两个类别具有非常小的运输效果，但大多数类别构成了一组大多数相似的拟合。
Partial Pooling -变坡度和截距
允许坡度和截距按类别变化的最一般方法。方程如下：
变坡截距模型：
代码语言：javascript复制varying_intercept_slope = """
data {
int<lower=0> N;
int<lower=0> J;
vector[N] y;
vector[N] x;
int category[N];
}
parameters {
real<lower=0> sigma;
real<lower=0> sigma_a;
real<lower=0> sigma_b;
vector[J] a;
vector[J] b;
real mu_a;
real mu_b;
}
model {
mu_a ~ normal(0, 100);
mu_b ~ normal(0, 100);
a ~ normal(mu_a, sigma_a);
b ~ normal(mu_b, sigma_b);
y ~ normal(a[category] b[category].*x, sigma);
}
"""
拟合模型：
代码语言：javascript复制varying_intercept_slope_data = {'N': len(log_price),
'J': len(n_category),
'category': category 1, # Stan counts starting at 1
'x': shipping,
'y': log_price}
sm = pystan.StanModel(model_code=varying_intercept_slope)
varying_intercept_slope_fit = sm.sampling(data=varying_intercept_slope_data,
iter=1000, chains=2)
拟合模型可视化:
代码语言：javascript复制plt.figure(figsize=(10, 6))
xvals = np.arange(2)
b = varying_intercept_slope_fit['a'].mean(axis=0)
m = varying_intercept_slope_fit['b'].mean(axis=0)
for bi,mi in zip(b,m):
plt.plot(xvals, mi*xvals bi, 'bo-', alpha=0.4)
plt.xlim(-0.1, 1.1);
plt.xticks([0, 1])
plt.title("fitted relationships by category")
plt.xlabel("shipping")
plt.ylabel("log(price)");
<p style="text-align: center;"><img src="http://imgcdn.atyun.com/2019/09/23.png" alt="(PyStan)零售价格贝叶斯策略建模（下）" width="933" height="585"></p>
虽然这些关系都非常相似，但我们可以看到，通过允许运输效果和价格变化，与不同的拦截模型相比，我们似乎能够捕获更多的自然变化。
语境效果
在某些情况下，在多个层次上使用预测因子可以揭示单个层次变量和组残差之间的相关性。我们可以通过将单个预测因子的平均值作为协变量包含在组截距模型中来解释这一点。
语境效应模型:
代码语言：javascript复制# Create new variable for mean of shipping across categories
xbar = df.groupby('category_name')['shipping'].mean().rename(category_lookup).values
x_mean = xbar[category]
contextual_effect = """
data {
int<lower=0> J;
int<lower=0> N;
int<lower=1,upper=J> category[N];
vector[N] x;
vector[N] x_mean;
vector[N] y;
}
parameters {
vector[J] a;
vector[3] b;
real mu_a;
real<lower=0,upper=100> sigma_a;
real<lower=0,upper=100> sigma_y;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- a[category[i]] x[i]*b[2] x_mean[i]*b[3];
}
model {
mu_a ~ normal(0, 1);
a ~ normal(mu_a, sigma_a);
b ~ normal(0, 1);
y ~ normal(y_hat, sigma_y);
}
"""
拟合模型：
代码语言：javascript复制contextual_effect_data = {'N': len(log_price),
'J': len(n_category),
'category': category 1, # Stan counts starting at 1
'x_mean': x_mean,
'x': shipping,
'y': log_price}
sm = pystan.StanModel(model_code=contextual_effect)
contextual_effect_fit = sm.sampling(data=contextual_effect_data,
iter=1000, chains=2)
预测
我们想预测“女式/运动服装/裤子、紧身裤、紧身裤”类的新产品，由卖家支付运费，我们只需要从模型中选取合适的截距。
代码语言：javascript复制category_lookup['Women/Athletic Apparel/Pants, Tights, Leggings']
预测模型：
代码语言：javascript复制contextual_pred = """
data {
int<lower=0> J;
int<lower=0> N;
int<lower=0,upper=J> wa;
real xbar_wa;
int<lower=1,upper=J> category[N];
vector[N] x;
vector[N] x_mean;
vector[N] y;
}
parameters {
vector[J] a;
vector[3] b;
real mu_a;
real<lower=0,upper=100> sigma_a;
real<lower=0,upper=100> sigma_y;
}
transformed parameters {
vector[N] y_hat;
real wa_mu;
for (i in 1:N)
y_hat[i] <- a[category[i]] x[i] * b[2] x_mean[i] * b[3];
wa_mu <- a[wa 1] b[2] xbar_wa * b[3];
}
model {
mu_a ~ normal(0, 1);
a ~ normal(mu_a, sigma_a);
b ~ normal(0, 1);
y ~ normal(y_hat, sigma_y);
}
generated quantities {
real y_wa;
y_wa <- normal_rng(wa_mu, sigma_y);
}
"""
做出预测：
代码语言：javascript复制contextual_pred_data = {'N': len(log_price),
'J': len(n_category),
'category': category 1,
'x_mean': x_mean,
'x': shipping,
'y': log_price,
'wa': 13,
'xbar_wa': xbar[13]}
sm = pystan.StanModel(model_code=contextual_pred)
contextual_pred_fit = sm.sampling(data=contextual_pred_data,
iter=1000, chains=2)
预测：
代码语言：javascript复制contextual_pred_fit.plot('y_wa');
观察：
本次fit抽样的平均值≈3，所以我们可以预计，在卖方支付运费时，“女装/运动服装/裤子、紧身衣、打底裤”类新产品的实测价格≈exp(3)≈20.09，虽然预测值的范围比较大。
Jupyter笔记本可以在Github上找到。好好享受周末余下的时光吧。
参考文献：
1.https://github.com/widdowquinn/Teaching-Stan-Hierarchical-Modelling?source=post_page—–fd0571ed778———————-
2.https://mc-stan.org/users/documentation/case-studies/radon.html?source=post_page—–fd0571ed778———————-
原文链接：
https://towardsdatascience.com/bayesian-strategy-for-modeling-retail-price-with-pystan-fd0571ed778
end
