这段代码通过Metropolis-Hastings算法生成了100,000个样本，这些样本服从一个混合高斯分布。它首先定义了一个概率密度函数，然后进行迭代，每次迭代中根据概率密度函数的比较决定是否接受新的样本。最后，通过三个图形展示采样轨迹、样本分布直方图以及目标概率密度函数，直观地展示了采样过程和结果。
n=100000; t=zeros(2,n);
for i=2:n
u=randn(2,1)+t(:,i-1);
if rand<=pdf(u)/pdf(t(:,i-1))
t(:,i)=u;
else
t(:,i)=t(:,i-1);
end
end
figure(1); clf; hold on
plot(t(1,1:500),t(2,1:500),'b-');
figure(2); clf; hold on
c=4; m=30; Q=zeros(m,m);
a=linspace(-c,c,m); b=-c+2*c/m*[1:m];
for k=1:n
x=find(ceil(t(1,k)-b)==0,1);
y=find(ceil(t(2,k)-b)==0,1);
Q(x,y)=Q(x,y)+1;
end
surf(a,a,Q'/(2*c/m)^2/n);
figure(3); clf; hold on
P=zeros(m,m);
for X=1:m, for Y=1:m
P(X,Y)=pdf([a(X);a(Y)]);
end, end
surf(a,a,P');
function p=pdf(x)
S1=[0.5 0; 0 4]; m1=x-[0; -2]; w1=0.7;
S2=[4 0; 0 1]; m2=x-[2; 0]; w2=0.3;
p=w1*sqrt(det(S1))/(2*pi)*exp(-m1'*S1*m1/2) ...
+w2*sqrt(det(S2))/(2*pi)*exp(-m2'*S2*m2/2);
直实概率密度：
生成样本直方图：
生成的采样轨迹：
以下是每一行代码的功能解释：
1. `n=100000; t=zeros(2,n);`: 定义了变量 `n` 并赋值为 100,000，同时创建了一个大小为 2x100,000 的零矩阵 `t`。
3. `for i=2:n`: 开始一个循环，从 2 到 `n`。
4. `u=randn(2,1)+t(:,i-1);`: 生成一个大小为 2x1 的随机向量 `u`，其中的元素是从标准正态分布中随机采样得到的，然后将其与 `t` 的上一列相加。
6. `if rand<=pdf(u)/pdf(t(:,i-1))`: 判断随机数 `rand` 是否小于等于 `u` 的概率密度函数值除以 `t` 的上一列的概率密度函数值。
7. `t(:,i)=u;`: 如果条件满足，将 `u` 赋值给 `t` 的当前列。
9. `else`: 如果条件不满足，执行下面的代码。
10. `t(:,i)=t(:,i-1);`: 将 `t` 的上一列赋值给 `t` 的当前列。
12. `figure(1); clf; hold on`: 创建一个名为 "figure(1)" 的图形窗口，并清空其中的内容，同时保持现有图形。
13. `plot(t(1,1:500),t(2,1:500),'b-');`: 在图形窗口中绘制 `t` 的前 500 列的第一行作为 x 坐标，前 500 列的第二行作为 y 坐标，使用蓝色实线进行连接。
15. `figure(2); clf; hold on`: 创建一个名为 "figure(2)" 的图形窗口，并清空其中的内容，同时保持现有图形。
16. `c=4; m=30; Q=zeros(m,m);`: 定义变量 `c` 并赋值为 4，定义变量 `m` 并赋值为 30，创建一个大小为 `m`x`m` 的零矩阵 `Q`。
17. `a=linspace(-c,c,m); b=-c+2*c/m*[1:m];`: 生成一个由 `m` 个等间距元素组成的从 `-c` 到 `c` 的行向量 `a`，生成一个由 `m` 个等间距元素组成的从 `-c` 到 `c` 的行向量 `b`。
18. `for k=1:n`: 开始一个循环，从 1 到 `n`。
19. `x=find(ceil(t(1,k)-b)==0,1);`: 找到向量 `t` 的第一行的第 `k` 列减去向量 `b` 向
上取整后等于 0 的第一个索引，并将结果赋值给 `x`。
20. `y=find(ceil(t(2,k)-b)==0,1);`: 找到向量 `t` 的第二行的第 `k` 列减去向量 `b` 向上取整后等于 0 的第一个索引，并将结果赋值给 `y`。
21. `Q(x,y)=Q(x,y)+1;`: 将 `Q` 的第 `x` 行、第 `y` 列的元素加 1。
23. `surf(a,a,Q'/(2*c/m)^2/n);`: 在图形窗口中创建一个 3D 表面图，其中 `a` 是 x 和 y 坐标，`Q'/(2*c/m)^2/n` 是 z 坐标。
25. `figure(3); clf; hold on`: 创建一个名为 "figure(3)" 的图形窗口，并清空其中的内容，同时保持现有图形。
26. `P=zeros(m,m);`: 创建一个大小为 `m`x`m` 的零矩阵 `P`。
27. `for X=1:m, for Y=1:m`: 开始两个嵌套循环，分别从 1 到 `m`。
28. `P(X,Y)=pdf([a(X);a(Y)]);`: 将概率密度函数 `pdf` 在 `[a(X);a(Y)]` 处的值赋值给 `P` 的第 `X` 行、第 `Y` 列的元素。
31. `surf(a,a,P');`: 在图形窗口中创建一个 3D 表面图，其中 `a` 是 x 和 y 坐标，`P'` 是 z 坐标。
33. `function p=pdf(x)`: 定义一个名为 `pdf` 的函数，输入参数为 `x`。
35. `S1=[0.5 0; 0 4]; m1=x-[0; -2]; w1=0.7;`: 定义变量 `S1` 并赋值为矩阵 `[0.5 0; 0 4]`，定义变量 `m1` 并赋值为 `x` 减去向量 `[0; -2]`，定义变量 `w1` 并赋值为 0.7。
37. `S2=[4 0; 0 1]; m2=x-[2; 0]; w2=0.3;`: 定义变量 `S2` 并赋值为矩阵 `[4 0; 0 1]`，定义变量 `m2` 并赋值为 `x` 减去向量 `[2; 0]`，定义变量 `w2` 并赋值为 0.3。
39. `p=w1*sqrt(det(S1))/(2*pi)*exp(-m1'*S1*m1/2) ...`: 计算概率密度函数 `p` 的值，使用了高斯分布的公式。
至此，这段代码通过循环生成了
一个二维数据集 `t`，并在图形窗口中绘制了相关的可视化结果。
这段代码展示了如何使用Metropolis-Hastings采样方法来生成服从特定概率分布的样本。
首先，定义了一个函数`pdf(x)`，用于计算给定输入`x`的概率密度函数值。这个概率密度函数是一个混合高斯分布，由两个高斯分布组成。
然后，通过一个循环进行Metropolis-Hastings采样。在每次迭代中，生成一个随机向量`u`，该向量是当前样本`t`的前一个样本加上一个随机扰动。然后计算接受率，即当前样本的概率密度函数值与新样本的概率密度函数值之比。如果接受率大于一个随机数（在0和1之间均匀分布），则接受新样本，否则保持当前样本不变。
在循环结束后，我们得到了一系列采样点`t`，它们服从目标概率分布。这些采样点被用于绘制两个图形。
第一个图形（figure(1)）绘制了前500个采样点的二维轨迹。这样可以观察到样本是如何在空间中移动的。
第二个图形（figure(2)）绘制了采样点在二维网格上的分布。首先创建了一个二维计数矩阵`Q`，记录每个网格单元中的样本点数量。然后遍历所有采样点，根据采样点的值确定其所属网格单元，并增加相应的计数。最后，将计数矩阵除以总采样点数和单元面积，得到每个网格单元的概率估计，并使用`surf`函数绘制成一个表面图。
第三个图形（figure(3)）绘制了目标概率分布的真实概率密度函数。首先创建了一个与二维网格相同大小的矩阵`P`，用于存储每个网格单元的概率密度函数值。然后遍历所有网格单元，使用`pdf`函数计算对应网格点的概率密度函数值，并存储在矩阵`P`中。最后，使用`surf`函数绘制成一个表面图。
通过比较第二个图形和第三个图形，我们可以看到Metropolis-Hastings采样方法生成的样本分布逐渐接近于目标概率分布的真实概率密度函数。
扫一扫
5753
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
