Gibbs Gauss采样入门

2023-10-25 09:28:11 浏览数 (2)

Gibbs Gauss采样入门

引言

Gibbs采样是一种马尔可夫链蒙特卡洛（MCMC）方法，用于从一个高维的概率分布中采样。在多元统计学和机器学习领域广泛应用。本文将介绍Gibbs采样的概念和步骤，并通过一个简单的例子演示如何使用Gibbs采样来采样从高斯分布中。

Gibbs采样步骤

Gibbs采样适用于联合分布(conditional distribution)可解析求解的情况下。下面是Gibbs采样的基本步骤：

  1. 初始化模型参数：从联合分布中的任意初始点开始，得到一个初始的样本点。
  2. 对于每个参数，按照条件分布采样：通过固定其他参数的值，从条件分布中采样得到新的参数值。
  3. 重复步骤2，直到得到足够数量的样本点。
  4. 收集样本点：将步骤3中得到的样本点作为最终的采样结果。 下面将通过一个示例来展示如何使用Gibbs采样来从高斯分布中采样。

Gibbs采样示例：采样从高斯分布中

考虑一个二维的高斯分布： $$ P(x, y) = frac{1}{2pi sigma_{x}sigma_{y}sqrt{1-rho^2}} expleft(-frac{1}{2(1-rho^2)}left(frac{(x-mu_x)^2}{sigma_x^2} - frac{2rho(x-mu_x)(y-mu_y)}{sigma_xsigma_y} frac{(y-mu_y)^2}{sigma_y^2}right)right) $$ 其中，$x$和$y$是服从正态分布的随机变量，$mu_x, mu_y$分别是$x$和$y$的均值，$sigma_{x}, sigma_{y}$分别是$x$和$y$的标准差，$rho$是$x$和$y$的相关系数。 下面通过Gibbs采样来从上述二维高斯分布中采样。

步骤1：初始化参数

首先，我们初始化模型参数。选择一个任意的初始点 $(x_0, y_0)$。

代码语言：javascript复制
pythonCopy codeimport numpy as np
# 初始化参数
mu_x = 2
mu_y = 2
sigma_x = 1
sigma_y = 1
rho = 0.5
# 初始化初始点
x0 = 0
y0 = 0

步骤2：按条件分布采样

待采样的二维高斯分布是联合分布，我们需要分别采样$x$和$y$这两个变量。根据高斯分布的特性，我们可以得到给定某个变量的条件下另一个变量的条件分布。 对于$x$： $$ P(x|y) = Nleft(mu_x rhofrac{sigma_x}{sigma_y}(y - mu_y), (1 - rho^2)sigma_x^2right) $$ 对于$y$： $$ P(y|x) = Nleft(mu_y rhofrac{sigma_y}{sigma_x}(x - mu_x), (1 - rho^2)sigma_y^2right) $$ 根据上述条件分布，我们可以分别得到给定$y_0$的条件下，采样$x_1$的值，以及给定$x_1$的条件下，采样$y_1$的值。

代码语言：javascript复制
pythonCopy codedef sample_x_given_y(y, mu_x, mu_y, sigma_x, sigma_y, rho):
    """给定y的条件下，采样x"""
    mu = mu_x   rho * sigma_x / sigma_y * (y - mu_y)
    sigma = np.sqrt((1 - rho ** 2) * sigma_x ** 2)
    return np.random.normal(mu, sigma)
def sample_y_given_x(x, mu_x, mu_y, sigma_x, sigma_y, rho):
    """给定x的条件下，采样y"""
    mu = mu_y   rho * sigma_y / sigma_x * (x - mu_x)
    sigma = np.sqrt((1 - rho ** 2) * sigma_y ** 2)
    return np.random.normal(mu, sigma)

步骤3和4：重复采样和收集样本点

接下来，根据上述步骤2中得到的条件分布，进行重复采样，直到得到足够数量的样本点。

代码语言：javascript复制
pythonCopy code# 重复采样并收集样本点
num_samples = 1000
samples = []
# 初始化样本点
x = x0
y = y0
# 重复采样
for _ in range(num_samples):
    # 采样x
    x = sample_x_given_y(y, mu_x, mu_y, sigma_x, sigma_y, rho)
    # 采样y
    y = sample_y_given_x(x, mu_x, mu_y, sigma_x, sigma_y, rho)
    # 收集样本点
    samples.append((x, y))

在上述代码中，重复采样了1000个样本点，并将其存储在​​samples​​列表中。

结果可视化

最后，我们可以将采样结果绘制成散点图，用于可视化二维高斯分布的采样结果。

代码语言：javascript复制
pythonCopy codeimport matplotlib.pyplot as plt
# 提取样本点的x和y值
x_samples = [s[0] for s in samples]
y_samples = [s[1] for s in samples]
# 绘制散点图
plt.scatter(x_samples, y_samples, alpha=0.6)
plt.title("Gibbs Sampling - Sampling from Bivariate Gaussian Distribution")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

运行上述代码后，将会得到一个散点图，展示了从二维高斯分布中采样的结果。

结论

本文介绍了Gibbs采样的基本步骤，并通过一个简单的示例演示了如何使用Gibbs采样从二维高斯分布中进行采样。Gibbs采样是一种强大的采样方法，在统计学、机器学习等领域有广泛应用。通过了解和使用Gibbs采样，我们可以更好地理解和应用MCMC方法。

实际应用场景可以是一个学生成绩管理系统。下面是一个使用Python实现的示例代码：

代码语言：javascript复制
pythonCopy codeclass Student:
    def __init__(self, stu_id, name):
        self.stu_id = stu_id   # 学生学号
        self.name = name   # 学生姓名
        self.scores = {}   # 学生成绩字典，键为科目名，值为分数
    def add_score(self, subject, score):
        self.scores[subject] = score
    def get_average_score(self):
        if len(self.scores) == 0:
            return 0
        else:
            total_score = sum(self.scores.values())
            return total_score / len(self.scores)
class GradeManagement:
    def __init__(self):
        self.students = {}   # 学生字典，键为学生学号，值为学生对象
    def add_student(self, stu_id, name):
        student = Student(stu_id, name)
        self.students[stu_id] = student
    def add_score(self, stu_id, subject, score):
        if stu_id in self.students:
            student = self.students[stu_id]
            student.add_score(subject, score)
    def get_student_average_score(self, stu_id):
        if stu_id in self.students:
            student = self.students[stu_id]
            return student.get_average_score()
# 示例代码的使用
grade_management = GradeManagement()
# 添加学生
grade_management.add_student(1, "张三")
grade_management.add_student(2, "李四")
# 添加成绩
grade_management.add_score(1, "数学", 90)
grade_management.add_score(1, "英语", 85)
grade_management.add_score(2, "数学", 95)
grade_management.add_score(2, "英语", 92)
# 获取学生平均成绩
average_score_1 = grade_management.get_student_average_score(1)
average_score_2 = grade_management.get_student_average_score(2)
print("学生1的平均成绩：", average_score_1)
print("学生2的平均成绩：", average_score_2)

以上代码实现了一个简单的学生成绩管理系统，可以添加学生、添加成绩，并获取学生的平均成绩。

Gibbs采样是一种用于对多维概率分布进行采样的方法，它通过逐个变量更新的方式利用条件概率进行采样。虽然Gibbs采样在实际应用中有着广泛的应用，并且其算法简单易懂，但也存在一些缺点：

  1. 收敛速度慢：Gibbs采样通常需要进行多轮迭代才能得到准确的采样结果。尤其是当所采样的分布具有较高的维度时，收敛速度更加慢。
  2. 无法跳出局部极值：Gibbs采样是一种局部采样方法，每次迭代仅通过条件概率对一个变量进行更新。这种逐个变量更新的特性使得Gibbs采样容易陷入和停滞在局部极值点，无法快速获取全局最优解。
  3. 采样效率低：由于Gibbs采样每次仅更新一个变量，因此需要大量的迭代次数才能收敛。因而在实际应用中需要投入较多的计算资源和时间。 类似的方法包括Metropolis-Hastings采样和重要性采样。 Metropolis-Hastings采样是一种基于马尔可夫链的采样方法，通过引入接受-拒绝策略，可以在某些情况下跳出局部极值。相较于Gibbs采样，Metropolis-Hastings采样的优势是可以从一个概率分布中采样，而不仅限于条件概率分布。 重要性采样是一种通过引入权重的方式对概率分布进行采样的方法。它通过从一个已知的简单概率分布中采样，并利用被采样样本对目标分布进行加权，从而得到目标分布的采样。与Gibbs采样相比，重要性采样在一定程度上可以更有效地探索高概率区域，但也需要引入权重计算和进行权重调整，增加了计算复杂度。
第二期热点征文-编程语言

0 人点赞