本文仅用于阐明《统计计算》课程StatComp22034
R包中自制函数的背景知识与相关理论推导。
Individualized Regression For Clustering
模型概述
记 $y_i\in\mathbb R$ 为第 $i$ 个个体的响应变量(response)观测值,$x_i\in\mathbb R$ 为第 $i$ 个个体异质性的协变量(covariate)观测值,$z_i\in\mathbb R$ 为第 $i$ 个个体的群体同质性的协变量观测值。考虑一般的回归模型
其中 $\beta_i\in\mathbb R$ ($i=1,\dots,N$) 就衡量了不同个体的异质性,而 $\alpha\in\mathbb R$ 则是群体共享参数,$\varepsilon_i$ 是满足零均值,同方差的随机误差。
记 $\pmb\beta=(\beta_1,\dots,\beta_N)’$,$\pmb\gamma=(\gamma_1,\dots,\gamma_K)’$。定义
为了实现对个体的聚类,考虑目标函数为
其中 $\pmb y=(y_1,\dots,y_N)’$,$\pmb X=\text{diag}(x_1,\dots,x_N)$ 是 $N\times N$ 矩阵,$\pmb z=(z_1,\dots,z_N)’$,$\lambda>0$ 是惩罚系数。
参数优化
进一步, 可以使用ADMM算法解决这一问题. 改写优化问题为
则该问题的增广拉格朗日函数为
其中 $\rho>0$ 是事先设定的固定参数,$\pmb\nu\in\mathbb R^N$是拉格朗日乘子向量.
ADMM算法的好处在于把这个函数分成好几个容易求解的部分进行计算。具体来说,
- 设定初值 $\pmb\beta^{(0)}$,$\alpha^{(0)}$,再设定初值 $\pmb\delta^{(0)}$,$\pmb\gamma^{(0)}$ 和 $\pmb\nu^{(0)}$。
-
对第 $r+1$($r=0,1,\dots$)次迭代,按下述式子更新诸参数的值:
$\begin{array}{rl} \left(\begin{array}{c} \pmb\beta^{(r+1)}\\ \alpha^{(r+1)} \end{array}\right)=&\!\!\!\!\underset{\pmb\beta,\alpha}{\arg\min} ~L_\rho(\pmb\beta,\pmb\alpha,\pmb\delta^{(r)},\pmb\gamma^{(r)};\pmb\nu^{(r)}),\\ \left(\begin{array}{c} \pmb\delta^{(r+1)}\\ \pmb\gamma^{(r+1)} \end{array}\right)=&\!\!\!\!\underset{\pmb\delta,\pmb\gamma}{\arg\min} ~L_\rho(\pmb\beta^{(r+1)},\alpha^{(r+1)},\pmb\delta,\pmb\gamma;\pmb\nu^{(r)}),\\ \pmb\nu^{(r+1)}=&\!\!\!\!\pmb\nu^{(r)}+\rho(\pmb\beta^{(r+1)}-\pmb\delta^{(r+1)}). \end{array}$ -
重复迭代直至算法收敛或迭代次数过多。
$\dfrac{1}{N}\|\pmb\beta^{(r+1)}-\pmb\beta^{(r)}\|_2+\|\alpha^{(r+1)}-\alpha^{(r)}\|_2+\dfrac{1}{K}\|\pmb\gamma^{(r+1)}-\pmb\gamma^{(r)}\|_2<\varepsilon_1$ 或 $\|\pmb r^{(r+1)}-\pmb r^{(r)}\|_2<\varepsilon_2$, 其中 $\pmb r^{(r)}=\pmb\beta^{(r)}-\pmb\delta^{(r)}$。
$\pmb\beta$ 和 $\alpha$ 的更新
舍去与最小化函数无关的项,即
即求解方程
该方程有显式解
其中 $\widetilde{\pmb X}=(\pmb X,\pmb z)$ 是 $N\times(N+1)$ 矩阵,$\pmb I_N$ 是 $N\times N$ 单位矩阵。
$\pmb\delta$ 和 $\pmb\gamma$ 的更新
舍去与最小化函数无关的项,即
在给定 $\pmb\gamma$ 的条件下,注意到两部分对 $\pmb\delta$ 的计算是分离的,因此计算可知
注意这里若有多个 $k_0$ 满足条件,则取使得不等号左边最小的那个 $k_0$。
而当 $\pmb\delta$ 完成更新后,基于最小化 $D(\pmb\delta,\pmb\gamma)$ 可知,$\pmb\gamma^{(r+1)}$ 为 $\pmb\delta^{(r+1)}$ 的 $K$-mediods 结果。
EM Algorithm for Repeated Experiments
模型概述
记 $y_{ij}$ 为第 $i$ $(i=1,\dots,n)$ 个个体在第 $j$ $(j=1,\dots,t)$ 次重复实验中的反应时间,构建贝叶斯层次模型如下:
其中 $s_i=1$ 当且仅当第 $i$ 个个体患有精神分裂症,这些个体的平均反应时间比正常个体有 $\gamma>0$ 的偏移;$z_{ij}=1$ 当且仅当第 $i$ 个个体在第 $j$ 次重复实验中病症发作,发病率为 $\lambda$,从而平均以 $\beta>0$ 的弛豫时间才作出反应;$\alpha_i$ 是每个个体的平均反应时间,$\sigma_1^2$,$\sigma_2^2$ 分别是正常个体和病症发作个体反应时间的方差,$\tau_1^2$,$\tau_2^2$ 分别是正常个体和患病个体平均反应时间的方差.
记 $\pmb\theta=(\pmb\alpha’,\beta,\gamma,\mu,\lambda,\sigma_1^2,\sigma_2^2,\tau_1^2,\tau_2^2)$ 为模型的未知参数,其中 $\pmb\alpha=(\alpha_1,\dots,\alpha_n)’$,而 $z_{ij}$ 为隐变量。于是模型在完全数据下的似然函数为
其中 $f(\cdot\mid\mu,\sigma^2)$ 是正态分布 $N(\mu,\sigma^2)$ 的密度函数.
EM算法
考虑到模型中存在缺失数据 $\pmb z$,因此使用E-M算法如下。首先计算模型在完全数据下的对数似然函数
E步
注意到 $z_{ij}\mid\lambda,s_i\sim\text{Bernoulli}(\lambda s_i)$,因此
其中
M步
这样E步就完成了,接下来通过坐标下降法来依次更新 $\pmb\theta$ 的分量,通过最大化 $Q(\pmb\theta\mid\pmb\theta^{(r)})$ 可以得到相应参数的更新公式如下:
- 更新 $\lambda$:
$\displaystyle\lambda^{(r+1)}=\sum\limits_{i=1}^n\sum\limits_{j=1}^t\xi_{ij}^{(r)}\bigg/\sum\limits_{i=1}^n\sum\limits_{j=1}^ts_i.$ - 更新 $\pmb\alpha$:
$\alpha_i^{(r+1)}=\dfrac{\sum\limits_{j=1}^t\Big(\frac{(1-\xi_{ij}^{(r)})y_{ij}}{(\sigma_1^{(r)})^2}+\frac{\xi_{ij}^{(r)}(y_{ij}-\beta^{(r)})}{(\sigma_2^{(r)})^2}\Big)+\frac{(1-s_i)\mu^{(r)}}{(\tau_1^{(r)})^2}+\frac{s_i(\mu^{(r)}+\gamma^{(r)})}{(\tau_2^{(r)})^2}}{\sum\limits_{j=1}^t\Big(\frac{1-\xi_{ij}^{(r)}}{(\sigma_1^{(r)})^2}+\frac{\xi_{ij}^{(r)}}{(\sigma_2^{(r)})^2}\Big)+\frac{1-s_i}{(\tau_1^{(r)})^2}+\frac{s_i}{(\tau_2^{(r)})^2}},\quad i=1,\dots,n.$ - 更新 $\beta,\sigma_1^2,\sigma_2^2$:
$\displaystyle\beta^{(r+1)}=\sum\limits_{i=1}^n\sum\limits_{j=1}^t\xi_{ij}^{(r)}(y_{ij}-\alpha_i^{(r+1)})\bigg/\sum\limits_{i=1}^n\sum\limits_{j=1}^t\xi_{ij}^{(r)},$ $\displaystyle(\sigma_1^{(r+1)})^2=\sum\limits_{i=1}^n\sum\limits_{j=1}^t(1-\xi_{ij}^{(r)})(y_{ij}-\alpha_i^{(r+1)})^2\bigg/\sum\limits_{i=1}^n\sum\limits_{j=1}^t(1-\xi_{ij}^{(r)}),$ $\displaystyle(\sigma_2^{(r+1)})^2=\sum\limits_{i=1}^n\sum\limits_{j=1}^t\xi_{ij}^{(r)}(y_{ij}-\alpha_i^{(r+1)}-\beta^{(r+1)})^2\bigg/\sum\limits_{i=1}^n\sum\limits_{j=1}^t\xi_{ij}^{(r)}.$ - 更新 $\mu,\gamma,\tau_1^2,\tau_2^2$:
$\mu^{(r+1)}=\dfrac{\sum\limits_{i=1}^n(1-s_i)\alpha_i^{(r+1)}}{\sum\limits_{i=1}^n(1-s_i)},\quad\gamma^{(r+1)}=\dfrac{\sum\limits_{i=1}^ns_i\alpha_i^{(r+1)}}{\sum\limits_{i=1}^ns_i}-\mu^{(r+1)},$ $(\tau_1^{(r+1)})^2=\dfrac{\sum\limits_{i=1}^n(1-s_i)(\alpha_i^{(r+1)}-\mu^{(r+1)})^2}{\sum\limits_{i=1}^n(1-s_i)},\quad(\tau_2^{(r+1)})^2=\dfrac{\sum\limits_{i=1}^ns_i(\alpha_i^{(r+1)}-\mu^{(r+1)}-\gamma^{(r+1)})^2}{\sum\limits_{i=1}^ns_i}.$
声明
本文内容未经许可,严禁转载与挪用!