PMM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 泊松混合模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 泊松分布\n",
    "\n",
    "泊松分布（Poisson distribution）的概率分布函数为：\n",
    "\n",
    "$$P(X = k) = \\frac{\\lambda^{k}}{k !} e^{- \\lambda}, k = 0, 1, \\cdots \\tag {1}$$\n",
    "\n",
    "泊松分布的参数$\\lambda$表示单位时间（单位面积）内随机事件的平均发生次数。泊松分布适合于描述单位时间内随机事件发生的次数。\n",
    "\n",
    "泊松分布的期望$\\lambda$：，方差：$\\lambda$，特征函数：$\\phi(t) = \\exp(\\lambda (e^{i t} - 1))$\n",
    "\n",
    "### 泊松分布与二项分布\n",
    "\n",
    "当二项分布的$n$很大而$p$很小时，二项分布可近似为泊松分布，其中$\\lambda = np$。通常当$n \\geq 20, p \\leq 0.05$时，就可以用泊松分布近似。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 泊松混合模型\n",
    "\n",
    "泊松混合模型（Poisson mixture model）的概率密度函数为：\n",
    "\n",
    "$$p(y; \\theta) = \\sum_{k}^{K} \\alpha_{k} \\frac{\\lambda^{k}}{k !} e^{- \\lambda} \\tag {2}$$\n",
    "\n",
    "其中，$\\alpha_{k} \\geq 0$，$\\sum_{k}^{K} \\alpha_{k} = 1$。\n",
    "\n",
    "给定由$K$个Poisson分布生成的$N$条独立数据构成数据集$D = {y_{j} | j = 1, 2, \\cdots, N}$，估计泊松混合模型的参数$\\theta_{k} = (\\alpha_{k}, \\lambda_{k} | k = 1, 2, \\cdots, K)$。\n",
    "\n",
    "1. 隐变量\n",
    "\n",
    "观测数据$y_{j}$已知，反映观测数据$y_{j}$来自第$k$个分模型的数据未知，以隐变量$\\gamma_{k}$表示：\n",
    "\n",
    "$$\\gamma_{jk} = \\begin{cases}\n",
    "1, & \\text{第}j\\text{个观测来自第}k\\text{个分模型} \\\\\n",
    "0, & 否则\n",
    "\\end{cases}\n",
    "j = 1, 2, \\cdots, N; k = 1, 2, \\cdots, K \\tag {3}$$\n",
    "\n",
    "未观测数据$\\gamma_{jk}$为`0-1`随机变量，则完全数据为\n",
    "\n",
    "$$(y_{j}, \\gamma_{j1}, \\cdots, \\gamma_{jK}), j = 1, 2, \\cdots, N$$\n",
    "\n",
    "完全数据的似然函数：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "P(y, \\gamma; \\theta) & = \\prod_{j = 1}^{N} P(y_{j}, \\gamma_{j1}, \\cdots, \\gamma_{jK}; \\theta) \\\\\n",
    "& = \\prod_{j = 1}^{N} \\prod_{k = 1}^{K} [\\alpha_{k} P(y_{j}; \\theta_{k})]^{\\gamma_{jk}} \\\\\n",
    "& = \\prod_{k = 1}^{K} \\alpha_{k}^{n_{k}} \\prod_{j = 1}^{N} [P(y_{j}; \\theta_{k})]^{\\gamma_{jk}} \\\\\n",
    "& = \\prod_{k = 1}^{K} \\alpha_{k}^{n_{k}} \\prod_{j = 1}^{N} \\left[\n",
    "    \\frac{\\lambda_{k}^{y_{j}}}{y_{j} !} e^{- \\lambda_{k}}\n",
    "\\right]^{\\gamma_{jk}}\n",
    "\\end{aligned}$$\n",
    "\n",
    "其中，$n_{k} = \\sum_{j = 1}^{N} \\gamma_{jk}$，$N = \\sum_{k = 1}^{K} n_{k}$。完全数据的对数似然函数为\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\log P(y, \\gamma; \\theta) = \\sum_{k = 1}^{K} \\left( n_{k} \\log \\alpha_{k} + \\sum_{j = 1}^{N} \\gamma_{jk} \\left[\n",
    "    y_{j} \\log \\lambda_{k} - \\log y_{j}! - \\lambda_{k}\n",
    "\\right] \\right) \\\\\n",
    "\\end{aligned} \\tag {4}$$\n",
    "\n",
    "2. E步，$Q$函数\n",
    "\n",
    "$$\\begin{aligned}\n",
    "Q(\\theta, \\theta^{(i)})\n",
    "& = \\text{E}_{\\gamma} \\left[ \\log P(y, \\gamma; \\theta) | y; \\theta^{(i)} \\right] \\\\\n",
    "& = \\sum_{k = 1}^{K} \\left(\n",
    "    \\sum_{j = 1}^{N} \\text{E}[\\gamma_{jk} | y_{j}; \\theta^{(i)}] \\log \\alpha_{k} + \\sum_{j = 1}^{N} \\text{E}[\\gamma_{jk} | y_{j}; \\theta^{(i)}] \\left[\n",
    "    y_{j} \\log \\lambda_{k} - \\log y_{j}! - \\lambda_{k}\n",
    "\\right] \\right) \\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "令$\\hat{\\gamma}_{jk} = \\text{E}[\\gamma_{jk} | y_{j}; \\theta^{(i)}]$，表示当前模型参数下第$j$个观测数据来自第$k$个分模型的概率，称为分模型$k$对观测数据$y_{j}$的响应度。\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\hat{\\gamma}_{jk} & = \\text{E}[\\gamma_{jk} | y_{j}; \\theta^{(i)}] \\\\\n",
    "& = \\frac{\n",
    "    P(y_{j} | \\gamma_{jk} = 1; \\theta^{(i)}) P(\\gamma_{jk} = 1; \\theta^{(i)})\n",
    "}{\n",
    "    \\sum_{k = 1}^{K} P(y_{j} | \\gamma_{jk} = 1; \\theta^{(i)}) P(\\gamma_{jk} = 1; \\theta^{(i)})\n",
    "} \\\\\n",
    "& = \\frac{\n",
    "    \\alpha_{k}^{(i)} \\frac{\\left( \\lambda_{k}^{(i)} \\right)^{y_{j}}}{y_{j} !} e^{- \\lambda_{k}^{(i)}}\n",
    "}{\n",
    "    \\sum_{k = 1}^{K} \\alpha_{k}^{(i)} \\frac{ \\left( \\lambda_{k}^{(i)} \\right)^{y_{j}}}{y_{j} !} e^{- \\lambda_{k}^{(i)}}\n",
    "} \\\\\n",
    "\\end{aligned} \\tag {5}$$\n",
    "\n",
    "则$n_{k} = \\sum_{j = 1}^{N} E[\\hat{\\gamma}_{jk}]$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "Q(\\theta, \\theta^{(i)}) = \\sum_{k = 1}^{K} \\left(\n",
    "    \\sum_{j = 1}^{N} \\hat{\\gamma}_{jk} \\log \\alpha_{k} + \\sum_{j = 1}^{N} \\hat{\\gamma}_{jk} \\left[\n",
    "    y_{j} \\log \\lambda_{k} - \\log y_{j}! - \\lambda_{k}\n",
    "\\right] \\right) \\\\\n",
    "\\end{aligned} \\tag {6}$$\n",
    "\n",
    "3. M步\n",
    "\n",
    "求使函数$Q(\\theta, \\theta^{(i)})$取极大值时的参数$\\theta^{(i + 1)} = (\\hat{\\lambda}_{k}, \\hat{\\alpha}_{k} | k = 1, 2, \\cdots, K)$：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\theta^{(i + 1)} & = \\argmax_{\\theta} Q(\\theta, \\theta^{(i)}) \\\\\n",
    "& \\text{s.t.} \\sum_{k = 1}^{K} \\alpha_{k} = 1, \\ \\alpha_{k} \\geq 1\n",
    "\\end{aligned}$$\n",
    "\n",
    "通过拉格朗日乘数法可得：\n",
    "\n",
    "$$\\hat{\\lambda}_{k} = \\frac{\\sum_{j = 1}^{N} \\hat{\\gamma}_{jk} y_{j}}{\\sum_{j = 1}^{N} \\hat{\\gamma}_{jk}} \\tag {7}$$\n",
    "\n",
    "$$\\hat{\\alpha}_{k} = \\frac{n_k}{N} = \\frac{\\sum_{j = 1}^{N} \\hat{\\gamma}_{jk}}{N} \\tag {8}$$\n",
    "\n",
    "### 泊松混合模型EM参数估计\n",
    "\n",
    "输入：观测数据$y_{1}, y_{2}, \\cdots,y_{N}$，泊松混合模型；\n",
    "\n",
    "输出：泊松混合模型参数。\n",
    "\n",
    "（1）参数初始化\n",
    "\n",
    "（2）E步：依据当前模型参数，计算分模型$k$对观测数据$y_{j}$的响应度\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\hat{\\gamma}_{jk} = \\frac{\n",
    "    \\alpha_{k}^{(i)} \\frac{\\left( \\lambda_{k}^{(i)} \\right)^{y_{j}}}{y_{j} !} e^{- \\lambda_{k}^{(i)}}\n",
    "}{\n",
    "    \\sum_{k = 1}^{K} \\alpha_{k}^{(i)} \\frac{ \\left( \\lambda_{k}^{(i)} \\right)^{y_{j}}}{y_{j} !} e^{- \\lambda_{k}^{(i)}}\n",
    "} \\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "（3）M步：计算新一轮迭代的模型参数\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\hat{\\lambda}_{k} & = \\frac{\\sum_{j = 1}^{N} \\hat{\\gamma}_{jk} y_{j}}{\\sum_{j = 1}^{N} \\hat{\\gamma}_{jk}} \\\\\n",
    "\\hat{\\alpha}_{k} & = \\frac{n_k}{N} = \\frac{\\sum_{j = 1}^{N} \\hat{\\gamma}_{jk}}{N} \\\\\n",
    "k & = 1, 2, \\cdots, K\n",
    "\\end{aligned}$$\n",
    "\n",
    "（4）重复（2）、（3），直至收敛（对数似然函数值不再有明显的变化）\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}