非参数贝叶斯

非参数贝叶斯方法(Nonparametric Bayesian Methods)是一类允许参数空间随数据增长而无限扩展的贝叶斯模型。与固定维度的参数模型不同,非参数模型具有适应性复杂度,能够根据数据自动调整模型的复杂度。1

核心工具包括:

  • Dirichlet过程(DP):无限混合模型的基础
  • Chinese Restaurant Process(CRP):聚类的非参数先验
  • Indian Buffet Process(IBP):稀疏特征发现的非参数先验

1. Dirichlet过程

1.1 定义与构造

Dirichlet过程是定义在概率测度空间上的随机过程,是Dirichlet分布向无限维空间的推广。

定义:设 为基分布(base distribution), 为浓度参数(concentration parameter)。Dirichlet过程 满足:对于任意有限可测划分 ,有

Sethuraman的三参数stick-breaking构造:Dirichlet过程可以表示为离散概率测度的混合:

其中权重 通过stick-breaking过程生成:

stick-breaking的名称源于”折棍子”的过程:取长度为1的棍子,每次折断长度为 的段,剩余部分继续折断。

import numpy as np
 
def stick_breaking(alpha, n_components):
    """Stick-breaking构造生成Dirichlet过程样本"""
    V = np.random.beta(1, alpha, n_components)
    # 确保和为1,通过最后一维取1-sum
    pi = V * np.concatenate([np.ones(n_components-1), [1]])
    pi[1:] = pi[1:] * np.cumprod(1 - V[:-1])
    return pi
 
# 生成样本
alpha = 1.0  # 浓度参数
n_samples = 1000
weights = stick_breaking(alpha, n_samples)

Polya Urn构造:另一种等价的构造方式。设颜色对应atom ,从 urn 中采样:

这意味着新样本以概率 从基分布采样,以概率 复制已有观测( 为该atom的计数)。

Dirichlet分布序列极限:设 ,边缘化 可得2

时,经验分布 弱收敛到 ,但 本身保持随机性。

1.2 DP的后验推断

Dirichlet过程的后验同样是Dirichlet过程。给定观测数据

后验预测分布:新观测 的预测分布为:

这正是 Polya Urn 的采样规则。设已有唯一值 对应计数 ,则:

当新值被采样时,其分布为基分布

1.3 共轭先验与后验

当似然函数与基分布共轭时,后验推断可以高效计算。常见的共轭对:

似然函数基分布 后验参数更新
BernoulliBetaBeta先验
CategoricalDirichletDirichlet先验
Gaussian (σ²已知)GaussianGaussian先验
PoissonGammaGamma先验
ExponentialGammaGamma先验

高斯共轭示例:设数据 已知,基分布 。后验为:

这使得 DP混合高斯模型(DP-GMM) 的Gibbs采样变得可行。


2. Chinese Restaurant Process

2.1 CRP与聚类的关系

Chinese Restaurant Process(CRP,中餐馆过程)是Dirichlet过程诱导的聚类先验,因其类比直观而广泛应用。

类比:想象进入一家有无限张桌子的中餐馆:

  • 第一位顾客随机选择一张桌子坐下
  • 位顾客选择:
    • 张已有桌子:概率 (当前顾客数)
    • 新桌子:概率 (浓度参数)

与DP的联系:设桌子对应atom ,顾客对应从分布 采样得到的观测。顾客坐在某桌子意味着该观测等于对应atom的值。

2.2 座位概率与后续分布

对于第 位顾客,座位概率为:

其中 是第 张桌子当前的顾客数, 是总顾客数。

座位概率的性质

  • rich-get-richer:顾客多的桌子更容易吸引新顾客(幂律特性)
  • 穷尽性:随着顾客数趋于无穷,聚类数增长率为
  • 一致性:CRP产生无限维的随机划分

2.3 预测分布公式

假设已有聚类配置 ,其中 标识第 个观测属于哪个聚类。预测第 个观测的聚类:

当引入基分布 时,完整预测分布为:

import numpy as np
from collections import Counter
 
def crp_predictive(n_k, alpha):
    """
    CRP预测分布
    n_k: dict,键为聚类ID,值为该聚类计数
    返回各聚类及新聚类的选择概率
    """
    n = sum(n_k.values())
    probs = {}
    
    # 已有聚类
    for k, count in n_k.items():
        probs[k] = count / (alpha + n)
    
    # 新聚类
    probs['new'] = alpha / (alpha + n)
    
    return probs
 
def sample_crp(alpha, n_samples):
    """从CRP采样"""
    clusters = Counter()
    assignments = []
    
    for i in range(n_samples):
        probs = crp_predictive(dict(clusters), alpha)
        labels = list(probs.keys())
        p = list(probs.values())
        
        choice = np.random.choice(labels, p=p)
        if choice == 'new':
            choice = max(clusters.keys(), default=0) + 1
        
        assignments.append(choice)
        clusters[choice] += 1
    
    return assignments

聚类数的期望:在CRP下,聚类数 的期望为3

其中 是第 个调和数, 是欧拉-马歇罗尼常数。


3. Indian Buffet Process

3.1 IBP定义与二项过程的关系

Indian Buffet Process(IBP,印度自助餐过程)是用于建模稀疏二值特征的非参数先验。

类比:进入一家有无限道菜的自助餐厅:

  • 顾客按编号顺序进入(
  • 位顾客从左到右浏览菜品:
    • 已有菜品:以概率 尝试( 为已尝试过的顾客数)
    • 新菜品:按Poisson(α/n)分布的数量新增
  • 顾客离开后,新菜品对后续顾客变为”已有菜品”

形式化定义:设 的二值矩阵,表示顾客与特征的关系。 表示第 个顾客具有第 个特征。

IBP分布定义为:

其中 是特征 的受欢迎程度, 是使用的特征总数。

与Beta-Bernoulli过程的联系:IBP等价于以下生成过程:

时, 的边际分布诱导出IBP。

3.2 特征矩阵表示

IBP产生的特征矩阵具有稀疏且右对齐的结构——特征按受欢迎程度从左到右排列。

矩阵的期望性质

  • 特征总数期望
  • 每行(顾客)特征数期望
  • 每列(特征)1的数量期望

典型矩阵形式):

特征  1 2 3 4 5 6 7 ...
顾客1  1 1 1 0 0 1 1 ...
顾客2  1 1 0 0 1 0 1 ...
顾客3  1 1 1 0 0 0 0 ...
...

呈现稀疏性(多数为0)且特征呈幂律分布(左边特征更常见)。

3.3 稀疏特征发现模型

IBP常作为稀疏特征发现的先验。典型应用是IBP线性因子模型

其中:

  • :稀疏二值特征
  • :特征对应的因子向量
  • :噪声

预测分布(顾客视角):已知前 位顾客的特征矩阵 ,第 位顾客的特征分布为:

import numpy as np
 
def ibp_sample(alpha, n_customers, max_features=100):
    """
    IBP采样
    alpha: 浓度参数
    n_customers: 顾客数量
    """
    Z = np.zeros((n_customers, max_features), dtype=int)
    m = np.zeros(max_features)  # 每列1的数量
    n_assigned = np.zeros(max_features, dtype=int)  # 每列被分配的顾客数
    
    K_observed = 0  # 当前观察到的特征数
    
    for n in range(1, n_customers + 1):
        # 处理已有特征
        for k in range(K_observed):
            prob = m[k] / n
            Z[n-1, k] = np.random.binomial(1, prob)
        
        # 新增特征(Poisson(α/n))
        n_new = np.random.poisson(alpha / n)
        for i in range(n_new):
            if K_observed + i < max_features:
                Z[n-1, K_observed + i] = 1
        
        K_observed += n_new
        
        # 更新计数
        m[:K_observed] += Z[n-1, :K_observed]
        n_assigned[:K_observed] += 1
    
    # 移除空列
    Z = Z[:, :K_observed]
    return Z
 
def ibp_posterior_predictive(Z, n, alpha):
    """
    IBP后验预测分布
    Z: 当前特征矩阵
    n: 当前顾客数(Z的行数)
    """
    m = Z.sum(axis=0)  # 每列的1的数量
    
    probs_existing = m / n
    prob_new = alpha / n
    
    return probs_existing, prob_new

应用场景

  • 协同过滤:发现用户的隐式偏好特征
  • 基因表达分析:识别功能基因模块
  • 文本主题建模:发现稀疏主题
  • 神经网络剪枝:识别重要连接模式

4. 应用案例

4.1 无限混合高斯模型(DP-GMM)

Dirichlet Process Gaussian Mixture Model(DP-GMM)将DP作为混合权重的先验,实现自动确定聚类数

模型定义

其中 是高斯-Wishart先验:

等效的CRP视角:聚类分配 服从CRP():

当创建新聚类时,从 采样其参数。

import numpy as np
from scipy.stats import wishart, norm
 
class DPGMM:
    """Dirichlet Process Gaussian Mixture Model"""
    
    def __init__(self, alpha, mu_0, kappa_0, nu_0, Lambda_0):
        self.alpha = alpha
        self.prior_params = (mu_0, kappa_0, nu_0, Lambda_0)
        self.clusters = {}  # 存储各聚类的参数
        self.assignments = []
    
    def sample_params(self):
        """从先验采样混合成分参数"""
        mu_0, kappa_0, nu_0, Lambda_0 = self.prior_params
        # Wishart采样
        Sigma = wishart.rvs(df=nu_0, scale=Lambda_0)
        # 高斯采样
        mu = norm.rvs(loc=mu_0, scale=Sigma / kappa_0)
        return {'mu': mu, 'Sigma': Sigma}
    
    def conditional_prior(self, k):
        """创建新聚类k的参数"""
        if k not in self.clusters:
            self.clusters[k] = self.sample_params()
        return self.clusters[k]
    
    def sample_assignments(self, X):
        """Gibbs采样聚类分配"""
        n = len(X)
        for i in range(n):
            # 移除当前观测
            old_k = self.assignments[i] if i < len(self.assignments) else None
            if old_k in self.clusters:
                self.clusters[old_k]['n'] -= 1
                if self.clusters[old_k]['n'] == 0:
                    del self.clusters[old_k]
            
            # 计算各聚类的后验概率
            log_probs = {}
            for k, cl in self.clusters.items():
                if cl['n'] > 0:
                    log_prob = np.log(cl['n']) + norm.logpdf(
                        X[i], cl['mu'], np.sqrt(cl['Sigma']))
                    log_probs[k] = log_prob
            
            # 新聚类
            params = self.sample_params()
            log_prob_new = np.log(self.alpha) + norm.logpdf(
                X[i], params['mu'], np.sqrt(params['Sigma']))
            log_probs['new'] = log_prob_new
            
            # 归一化并采样
            max_log = max(log_probs.values())
            probs = np.exp([lp - max_log for lp in log_probs.values()])
            probs /= probs.sum()
            
            choices = list(log_probs.keys())
            new_k = np.random.choice(choices, p=probs)
            
            if new_k == 'new':
                new_k = max(self.clusters.keys(), default=0) + 1
                self.clusters[new_k] = params
                self.clusters[new_k]['n'] = 0
            
            self.clusters[new_k]['n'] += 1
            
            if i < len(self.assignments):
                self.assignments[i] = new_k
            else:
                self.assignments.append(new_k)
        
        return self.assignments

DP-GMM的特点

  • 自动聚类数:无需预先指定聚类数
  • 避免过拟合:通过控制聚类数的先验
  • 后验推断:常用Gibbs采样或变分推断

4.2 层级Dirichlet过程(HDP)

Hierarchical Dirichlet Process(HDP,层级Dirichlet过程)用于共享聚类的多组数据

问题背景:假设有 个文档组,每个组内部有无限个聚类(主题),但所有组共享同一个聚类基座(word types)。

模型结构

其中 是全局基分布, 是组的特定分布。

等效构造:HDP可写为stick-breaking形式:

其中

HDP-LDA应用:HDP最著名的应用是HDP-LDA主题模型4

class HDP_LDA:
    """Hierarchical Dirichlet Process LDA"""
    
    def __init__(self, gamma, alpha, eta):
        """
        gamma: top-level concentration
        alpha: second-level concentration
        eta: base distribution for word-topic
        """
        self.gamma = gamma
        self.alpha = alpha
        self.eta = eta
        self.topics = {}  # 全局atom字典 phi_k
        self.doc_topics = []  # 每篇文档的topic分配
    
    def _sample_global_topics(self):
        """采样全局主题基座 G_0"""
        # 使用stick-breaking
        # 实际实现需要更复杂的MCMC
        pass
    
    def conditional_doc_topic(self, doc_idx, word_idx, word):
        """
        文档-主题的条件分布
        p(z = k | rest) ∝ (n_{jk}^{-i} + απ_k) * p(w | φ_k)
        """
        pass
    
    def fit(self, documents, n_iter=100):
        """HDP-LDA Gibbs采样"""
        # 初始化
        # 迭代采样
        for it in range(n_iter):
            for doc_idx, doc in enumerate(documents):
                for pos, word in enumerate(doc):
                    # 采样主题分配
                    # 更新计数
                    pass

HDP的特点

  • 跨组共享聚类:允许信息在组间共享
  • 组内独立:每个组可以有自己独特的聚类比例
  • 应用:多文档主题建模、多任务学习、迁移学习

4.3 非参数神经网络的连接(Neural Processes)

Neural Processes(NPs,神经过程)是将随机过程神经网络结合的框架,深度对应DP的思想。

Neural Process定义:设 为上下文集, 为目标输入。NP预测

确定性NP(DNP)

变分NP(VNP):引入潜在变量

这与DP混合模型的结构相似—— 对应聚类/主题变量。

Conditional NP(CNP)

import torch
import torch.nn as nn
 
class Encoder(nn.Module):
    """编码器:映射 (x, y) -> r"""
    def __init__(self, input_dim, hidden_dim, r_dim):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, r_dim)
        )
    
    def forward(self, x, y):
        return self.net(torch.cat([x, y], dim=-1))
 
class Aggregator(nn.Module):
    """聚合:集合 -> 表示"""
    def forward(self, r):
        return r.mean(dim=0)  # 设置维度
 
class Decoder(nn.Module):
    """解码器:生成条件分布"""
    def __init__(self, r_dim, x_dim, hidden_dim):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(r_dim + x_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU()
        )
        self.mu_net = nn.Linear(hidden_dim, 1)
        self.sigma_net = nn.Linear(hidden_dim, 1)
    
    def forward(self, r, x_target):
        h = self.net(torch.cat([r.expand(x_target.size(0), -1), x_target], dim=-1))
        return self.mu_net(h), torch.exp(self.sigma_net(h))
 
class ConditionalNP(nn.Module):
    """Conditional Neural Process"""
    def __init__(self, x_dim, y_dim, r_dim, hidden_dim):
        super().__init__()
        self.encoder = Encoder(x_dim + y_dim, hidden_dim, r_dim)
        self.aggregator = Aggregator()
        self.decoder = Decoder(r_dim, x_dim, hidden_dim)
    
    def forward(self, context_x, context_y, target_x):
        # 编码上下文
        r = self.encoder(context_x, context_y)
        r = self.aggregator(r)
        # 解码目标
        mu, sigma = self.decoder(r, target_x)
        return torch.distributions.Normal(mu, sigma)

与DP的理论联系

  • 表示不确定性:NP通过潜在变量 表示函数空间的随机性
  • 数据自适应复杂度:随着上下文点增多,预测不确定性减小(类似DP后验)
  • 不变性:置换不变聚合对应DP的exchangeability

Attentive NP (ANP):引入注意力机制改进聚合:

这使得模型能够选择性关注相关上下文点,类似于CRP的rich-get-richer特性。


总结

非参数贝叶斯方法提供了一套优雅的工具来处理复杂度自适应的建模问题:

方法核心思想主要应用
DP无限离散分布无限混合模型
CRP聚类的exchangeable分配聚类、分类
IBP稀疏二值特征的发现特征选择、因子分析
HDP层级共享的聚类结构多文档主题建模
NP随机函数的神经表示元学习、函数逼近

这些方法的核心特性:

  • 适应性:复杂度随数据增长
  • 一致性:exchangeability保证统计性质
  • 可解释性:稀疏性、聚类结构提供可解释的表示

参考

Footnotes

  1. Teh, Y. W. (2010). Dirichlet processes. In Encyclopedia of Machine Learning. Springer.

  2. Blackwell, D., & MacQueen, J. B. (1973). Ferguson distributions via Pólya’s urn scheme. The Annals of Statistics, 2(2), 353-355.

  3. Antoniak, C. E. (1974). Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. The Annals of Statistics, 2(6), 1152-1174.

  4. Teh, Y. W., Jordan, M. I., Beal, M. J., & Blei, D. M. (2006). Hierarchical Dirichlet processes. Journal of the American Statistical Association, 101(476), 1566-1581.