【机器学习】Tensorflow概率编程：线性混合模型

　　线性模型是我们最常见到的、最理想的数学模型，基本的线性模型是数据科学入门的基本案例。然而现实生活中的线性问题，很大几率不适用于基本的线性模型，需要使用线性混合模型来描述。Tensorflow edward提供对这类问题的解决方案。
　　大部分数据科学教材头几章都会讲解线性回归等知识。这些方法基于简单的线性模型，有一些强假设，比如：只有目标特征会影响样本；目标特征彼此独立等。然而现实问题并不一定符合这样理想化的假设。譬如说，某类网站点击量会随着一年不同季节的变化而出现一定差异；调查某类人群的心理状况，单个人调查多次，可能出现不同的结果。显然，这些样本由一些主要的、明显独立的主要特征来决定，但又受到一些其他特征的干扰，这些特征可能彼此都不独立，也可能无法量化（比如人员ID）。这些可能需要使用更加复杂的模型来描述。
　　线性混合模型要解决的问题，是计算主要特征（固定特征\效应）与样本的关系，同时考虑一些次要特征（随机特征\效应）的影响。如果我们只计算固定特征，那么表明我们默认结果只和固定特征、目标样本有关，所做的计算完全去迎合我们的样本，这是一个非常强的假设，能否适用于其他样本其他情况，这个不予以考虑；如果我们补充考虑一些其他次要特征带来的影响，那么整个模型弹性度就上升，能适用于更广的范围，同时这也使得固定特征的解空间的范围变小，即解更加明确。那么观测上述作用，我们注意到这这些作用与某个概念的作用非常相似——正则化。
　　我们来看看线性混合模型的结构。线性模式的通用形式是$Y=WX+B$，为了对应Edward官网教程，也可写成$Y=X\beta$，统一了斜率项和截距项。而线性混合模型同时考虑了固定特征和随机特征。设随机特征为$Z$，那么写成$Y=WX+B_f+Z\eta+B_r$，后两项是随机特征；亦可写成$Y=X\beta+Z\eta,\eta\sim N(0,\sigma^2)$，对$\eta$的限定是为了控制随机项对结果的影响。一些教材中往往写成$Y=X\beta+Z\eta+\epsilon$，其中$\epsilon$指的是误差项。随机特征的考虑是比较复杂的，有相关性算法用于判别；工程上则比较考验工程师的经验，用经验判断往往有很好的效果。
那么Tensorflow Edward怎么实现线性混合模型呢？我们来看看官网给出的例子。
　　首先我们需要一块数据，这块数据叫InstEval，下载地址为https://github.com/blei-lab/edward/blob/master/examples/data/insteval.csv。这是一个关于教师教学与学生成绩的关系列表。数据里面$Y$对应的项就是y，即成绩等级/分数，s对应学生ID，dept对应部门ID，d对应教师ID，service表示是否进行教学，是一个二值项。我们将service设为固定特征。ID只是作为样本的序号来使用，与我们明显是不能作为特征输入的；我们知道不同的教师、部门等会产生影响，但是不明白不同部门、教师等具体有什么样的影响与关系，这些作为随机特征。为了考虑不同教师、部门的影响，ID必须经过数据预处理，变为category数据才能使用，即

data = pd.read_csv('../examples/data/insteval.csv')
data['dcodes'] = data['d'].astype('category').cat.codes
data['deptcodes'] = data['dept'].astype('category').cat.codes
data['s'] = data['s'] - 1

这些category项目并不带有真正的特征；我们没有证据认为教师和部门之间具有强的相关性，默认其彼此独立的；我们没有真正意义上的随机特征，因此随机项$Z\eta+B_r$前面的$Z\eta$就可以不要了，只保留截距项$B_r$，并且令随机项~$N(0,\sigma^2)$。由于信息并不丰富，只能做出上面非常保守的假设。最终我们的线性混合模型表述为：

$$Y=Service\beta_{service}+B_f+B_s+B_{dept}+B_{d}$$
根据这个模型，我们就可以利用Tensorflow Edward来建模求解，求解方法是KL变分推断。代码如下：

n_s = 2972  # 学生总数
n_d = 1128  # 老师总数
n_dept = 14  # 部门总数
# 数据输入
s_ph = tf.placeholder(tf.int32, [None])#学生编号category
d_ph = tf.placeholder(tf.int32, [None])#教师编号category
dept_ph = tf.placeholder(tf.int32, [None])#部门编号category
service_ph = tf.placeholder(tf.float32, [None])#二值项，固定特征

#固定特征参数项
mu = tf.Variable(tf.random_normal([]))#Bf
service = tf.Variable(tf.random_normal([]))#beta
#随机特征截距的参数
sigma_s = tf.sqrt(tf.exp(tf.Variable(tf.random_normal([]))))#学生Bs的方差
sigma_d = tf.sqrt(tf.exp(tf.Variable(tf.random_normal([]))))#教师Bd的方差
sigma_dept = tf.sqrt(tf.exp(tf.Variable(tf.random_normal([]))))#部门Bdept方差

# 随机特征截距
eta_s = Normal(loc=tf.zeros(n_s), scale=sigma_s * tf.ones(n_s))
eta_d = Normal(loc=tf.zeros(n_d), scale=sigma_d * tf.ones(n_d))
eta_dept = Normal(loc=tf.zeros(n_dept), scale=sigma_dept * tf.ones(n_dept))
#随机特征项+固定特征项
yhat = tf.gather(eta_s, s_ph) + \
    tf.gather(eta_d, d_ph) + \
    tf.gather(eta_dept, dept_ph) + \
    mu + service * service_ph#这里tf.gather实际作用是样本采样
y = Normal(loc=yhat, scale=tf.ones(n_obs))

这里各个随机特征的截距项是一个服从$N(0,\sigma^2)$的随机变量。然后使用KL变分推断方法进行计算：

q_eta_s = Normal(
    loc=tf.Variable(tf.random_normal([n_s])),
    scale=tf.nn.softplus(tf.Variable(tf.random_normal([n_s]))))
q_eta_d = Normal(
    loc=tf.Variable(tf.random_normal([n_d])),
    scale=tf.nn.softplus(tf.Variable(tf.random_normal([n_d]))))
q_eta_dept = Normal(
    loc=tf.Variable(tf.random_normal([n_dept])),
    scale=tf.nn.softplus(tf.Variable(tf.random_normal([n_dept]))))

latent_vars = {
    eta_s: q_eta_s,
    eta_d: q_eta_d,
    eta_dept: q_eta_dept}
data = {
    y: y_train,
    s_ph: s_train,
    d_ph: d_train,
    dept_ph: dept_train,
    service_ph: service_train}
inference = ed.KLqp(latent_vars, data)

训练数据都是采样好了的，tensorflow对用户自定义采样支持的很好。这里包含另外一个假设，那就是我们认为特征与观测样本的关系可以完全因子分解，用数个多维正态分布的组合能良好的拟合映射关系。
　　为了评估训练的好坏，我们需要计算当前模型的预测结果与实际结果的差值。Tensorflow Edward的语法还不算太完善。比如Edward目前没有原生的用于保存模型、参数的API，需要直接使用Tensorflow的API；Edward单独提供了copy函数，这个API是用于构建预测功能的Edward模块的，请注意python的=号表达式一般是传递引用，因此copy构建的预测模块可以随着训练的进行也同步优化。例如我们要使用训练好的模型来预测新数据：

yhat_test = ed.copy(yhat, {
    eta_s: q_eta_s.mean(),
    eta_d: q_eta_d.mean(),
    eta_dept: q_eta_dept.mean()})
yhat_test.eval(feed_dict={.......})#输入你的数据

一开始我们的参数都是随机的，差值分布非常随意:

训练完毕后，差值为0的部分非常集中，差值分布呈现类似正态分布的特点：

　　至此，大功告成，考虑进随机项之后，结果依然非常稳定。官网给出了与R语言的结果对比，差异非常小，Edward可以替代部分R的功能。
引用：
http://edwardlib.org/tutorials/linear-mixed-effects-models
Gelman, A., & Hill, J. L. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.