有人建议谈论shrinkage的时候，不可以缺少这个名场面，也许可以拿来当你的Graphical Abstract。 嗯，好主意，这里就直接当了标题。

Erik van Zwet和合作者最近一段时间(2020年末)在arXiv上连发表三篇关于shrinkage的三篇文章(shrinkage trilogy)。其中还有一篇还是和Andrew Gelman合作的。这篇推送将用贝叶斯的经典例子--抛硬币（哦不，应该用一个高大上的名称A/B testing）讲清楚什么是shrinkage。

什么是shrinkage？

Shrinkage is where extreme values in a sample are “shrunk” towards a central value, like the sample mean.

嗯，一句话搞定，下面的都是废话，可以滑到文末点在看了

A/B testing

AB Test 实验一般有 2 个目的：判断哪个更好：例如，有 2 个 UI 设计，究竟是 A 更好一些，还是 B 更好一些，我们需要实验判定 计算收益：例如，最近新上线了一个直播功能，那么直播功能究竟给平台带了来多少额外的 DAU，多少额外的使用时长，多少直播以外的视频观看时长等

https://www.zhihu.com/question/20045543

抛硬币太无聊了，换个有意思的例子。比如用自己当前的头像在某社交网络上搭讪了200个小姐姐，有20个加了你；用彭于晏的头像搭讪了100个小姐姐，有10个加了你，成功率都是0.1，但是今后应该使用哪个头像继续呢？国外某数据网站就使用了9000多个不同人的头像做了这样的测试，数据如下。其中AB是搭讪总次数，H代表成功次数，average是成功率。心中默念三遍，以免后面迷路。

我们可以看一下搭讪成功率最高的是谁：

df %>% arrange(desc(average)) %>% head(5) %>% paged_table()

再看一下搭讪成功率最低的是谁：

df %>% arrange(average) %>% head(5) %>% paged_table()

What? 直觉告诉你这不科学，学过的知识也告诉样本没有代表性。

选取超过500次的实验，做个Histogram看看：

df %>%
filter(AB >= 500) %>%
ggplot(aes(average)) +
geom_histogram(binwidth = .005)

Make Sense，那接下来就用这个数据吧(df_filtered)建模，当先验吧。

df_filtered <- df %>%filter(AB >= 500)

这样类似抛硬币的事情，一般用的是Beta分布建模

它有两个参数，这里使用MLE拟合找到这两个参数的值，得到α0为101.7319，β0为289.046。

library(stats4)
ll <- function(alpha, beta) {x <- df_filtered$Htotal <- df_filtered$AB-sum(VGAM::dbetabinom.ab(x, total, alpha, beta, log = TRUE))
}
# maximum likelihood estimation
m <- mle(ll, start = list(alpha = 1, beta = 10), method = "L-BFGS-B",lower = c(0.0001, .1))
ab <- coef(m)
alpha0 <- ab[1]
beta0 <- ab[2]> alpha0alpha
101.7319
> beta0beta
289.046

画个图看一下这个Beta分布

df_filtered  %>%ggplot() +geom_histogram(aes(average, y = ..density..), binwidth = .005) +stat_function(fun = function(x) dbeta(x, alpha0, beta0), color = "red",size = 1) +xlab("average")

使用得到的这个Beta分布作为先验分布，对每一个人的成功率计算其经验贝叶斯估计(empirical Bayes estimate)。相当于：我们已知成功率服从这样的一个分布Beta(101.4, 287.3)，基于先验信息，对于给定数据，计算他们的后验期望是多少。比如一个人的搭讪总次数是1000次，成功了300次，那么成功率的经验贝叶斯估计(期望E)就是:

为什么可以这么算

因为Beta分布的后验分布这么算:

将后验的α和β值代入计算Beta分布期望值的公式就得能得出以上结果。

对所有的人（包括搭讪次数少于500的人），计算他们成功率的贝叶斯估计：

df_eb <- df %>%mutate(eb_estimate = (H + alpha0) / (AB + alpha0 + beta0))

看一下根据贝叶斯估计，谁的成功率最高：

df_eb %>% arrange(desc(eb_estimate)) %>% head(5) %>% paged_table()

谁的成功率最低：

df_eb %>% arrange(eb_estimate) %>% head(5) %>% paged_table()

显然之前我们看到的那几个只搭讪了几次，成功率为0或者1的人并不在列。Rogers Hornsby有8000多的搭讪次数，获得了0.358的成功率，获得第一名当之无愧。使用先验信息+贝叶斯估计好像疗效不错~~

画个图比较一下直接计算的成功率和基于先验计算的成功率。

ggplot(df_eb, aes(average, eb_estimate, color = AB)) +geom_hline(yintercept = alpha0 / (alpha0 + beta0), color = "red", lty = 2) +geom_point() +geom_abline(color = "red") +scale_colour_gradient(trans = "log", breaks = 10 ^ (1:5)) +xlab("average") +ylab("Empirical Bayes average")

横轴(average)是无先验信息直接计算的成功率:H/AB

纵轴(Empirical Bayes average)是基于先验信息计算的贝叶斯估计: (H+α0)/(AB+α0+β0)

颜色的深浅(AB)代表了搭讪次数的多少。

水平红色虚线是先验的期望值α0/(α0+β0)=0.261，它代表的是在没有观测到数据时，基于先验信息我们的期望是多少。

红色斜线是对角线，落在这条线上的点表示两种计算方式得出的结果相同。

先看最两侧，搭讪次数10次以下，成功率离谱(0或者1)的几个点。由于搭讪次数较少，使用先验和贝叶斯进行计算的时候，仿佛先验在说，你们的搭讪次数太少了，没什么代表性，就暂且认为你们的成功率和我差不多(0.261)吧。这些点受到到先验的影响更多，距离对角线距离是更远的。

再看对角线附近的点会发现，搭讪的次数越多，采样更具有代表性，在使用先验和贝叶斯计算的时候，仿佛先验在说，考虑到你们搭讪次数比较多，采样还有点代表性，那我就多相信你们一点。这些点会较少受到先验信息的影响，距离对角线越近。

看到这里应该明白什么是shrinkage了：moving all our estimates towards the average。搞明白之后记得还要保持神秘，跟师妹说这是emprical Bayesian shrinkage toward a Beta prior~

那么问题来了，什么样的人搭讪成功率最高？

谷歌查了一下他们，结果发现好像全部都是打棒球的，事有蹊跷！！！

这篇文章基于http://varianceexplained.org/r/empirical_bayes_baseball/

整理，代码都源于该博文，这是如何获取所用的数据：

library(Lahman)
df<- Batting %>%filter(AB > 0) %>%anti_join(Pitching, by = "playerID") %>%group_by(playerID) %>%summarize(H = sum(H), AB = sum(AB)) %>%mutate(average = H / AB)

大部分的Emprical Bayes的例子都是抛硬币，给人的印像就是统计学家都是疯狂抛硬币的人。David Robinson使用了Baseball的例子，贴近生活也更易理解。该博主还有很多好文带你入门Empirical Bayes，比如怎么算FDR, 什么是credible intervals，怎么做Emprical Bayesian hierachical modeling等等，推荐前往阅读。如果有童鞋感兴趣撰文整理这些内容，欢迎投稿或者后台联系我~

—END—