matlab如何用吉布斯抽样，或者MCMC抽样，条件抽样等方法，去实现简单的二维分布抽样

假设二维分布是类似图1这种二元指数多项式，参数我是随意给的，啥都行，其图形大致为图二。x范围是0-20，y范围是0-2π
式子满足概率密度函数的特性，其x从0-20，y从0-2π的积分为1，因为函数比较复杂，累积分布无法直接积分得出式子
要如何用吉布斯抽样，或者MCMC抽样，或者条件抽样等效率高点的其他方法，去实现二维分布抽样，只要是效率高点的，可以对这个二维分布进行抽样的方法都可以
我之前用过接受拒绝抽样去抽样，但是接受率只有10%左右，效率太低了，得考虑其他方法

img


img

创建测试数据
作为第一步，我们创建一些测试数据，用于拟合我们的模型。让我们假设预测变量和响应变量之间存在线性关系，因此我们采用线性模型并添加一些噪声。

我将x值平衡在零附近以“去相关”斜率和截距。

trueA <- 5
trueB <- 0
trueSd <- 10
sampleSize <- 31
# 创建独立的x值
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
# 根据ax + b + N（0，sd）创建因变量
y <-  trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)
plot(x,y, main="Test Data")

img


定义统计模型
下一步是指定统计模型。我们已经知道数据是用x和y之间的线性关系y = a * x + b和带有标准差sd的正常误差模型N（0，sd）创建的，所以让我们使用相同的模型进行拟合，看看如果我们可以检索我们的原始参数值。

从模型中导出似然函数
为了估计贝叶斯分析中的参数，我们需要导出我们想要拟合的模型的似然函数。可能性是我们期望观察到的数据以我们所看到的模型的参数为条件发生的概率（密度）。因此，假设我们的线性模型y = b + a * x + N（0，sd）将参数（a，b，sd）作为输入，我们必须返回在此模型下获得上述测试数据的概率（这听起来更复杂，正如你在代码中看到的那样，我们只是计算预测y = b + a * x和观察到的y之间的差异，然后我们必须查找概率密度（使用dnorm）发生这种偏差

likelihood <- function(param){
    a = param[1]
    b = param[2]
    sd = param[3]
    pred = a*x + b
     sumll = sum(singlelikelihoods)
     (sumll)  
}
 slopevalues <- function(x){return(likelihood(c(x, trueB, trueSd)))}

img

斜率参数的对数似然曲线

作为说明，代码的最后几行绘制了斜率参数a的一系列参数值的似然性。 

为什么我们使用对数
您可能已经注意到我返回似然函数中概率的对数，这也是我对所有数据点的概率求和的原因（乘积的对数等于对数之和）。我们为什么要做这个？ 因为很多小概率乘以的可能性很快就会变得非常小（比如10 ^ -34）。在某些阶段，计算机程序正在进入数字舍入问题。 

定义先验
作为第二步，与贝叶斯统计中一样，我们必须为每个参数指定先验分布。为了方便起见，我对所有三个参数使用了均匀分布和正态分布。 无信息先验通常是1 / sigma的比例（如果你想了解原因，请看这里）。

#先验分布
prior <- function(param){
    a = param[1]
    b = param[2]
    sd = param[3]
    aprior =  (a, min=0, max=10, log = T)
    bprior = dnorm(b, sd = 5, log = T)
 }

后验
先验和可能性的乘积是MCMC将要处理的实际数量。这个函数被称为后验 。同样，在这里我们使用总和，因为我们使用对数。

posterior <- function(param){
   return ( (param) + prior(param))
}

MCMC
现在，实际上是Metropolis-Hastings算法。该算法最常见的应用之一（如本例所示）是从贝叶斯统计中的后验密度中提取样本。然而，原则上，该算法可用于从任何可积函数中进行采样。因此，该算法的目的是在参数空间中跳转，但是以某种方式使得在某一点上的概率与我们采样的函数成比例（这通常称为目标函数）。在我们的例子中，这是上面定义的后验。

这是通过

  • 从随机参数值开始
  • 根据称为提议函数的某个概率密度，选择接近旧值的新参数值
  • 以概率p（新）/ p（旧）跳到这个新点，其中p是目标函数，p> 1表示跳跃

当我们运行这个算法时，它访问的参数的分布会收敛到目标分布p。那么，让我们在R中得到 ：

########Metropolis算法# ################

proposalfunction <- function(param){
    return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
      for (i in 1:iterations){
         if (runif(1) < probab){
            chain[i+1,] = proposal
        }else{
            chain[i+1,] = chain[i,]
        }
    }
    return(chain)
}
 chain = run_metropolis_MCMC(startvalue, 10000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))

使用后验的对数可能在开始时有点混乱，特别是当您查看计算接受概率的行时（probab = exp（后验（建议） - 后验（链[i，]）） ）。要理解我们为什么这样做，请注意p1 / p2 = exp [log（p1）-log（p2）]。

算法的第一步可能受初始值的偏差，因此通常被丢弃用于进一步分析 。要看的一个有趣的输出是接受率： 接受标准拒绝提案的频率是多少？接受率可以受提议函数的影响：通常，提案越接近，接受率越大。然而，非常高的接受率通常是无益的：这意味着算法“停留”在同一点 。可以证明，20％到30％的接受率对于典型应用来说是最佳的

###概要#######################
par(mfrow = c(2,3))
hist( [-(1:burnIn),1],nclass=30, , main="Posterior of a", xlab="True value = red line" )
abline(v = mean(chain[-(1:burnIn),1]))
#进行比较：
summary(lm(y~x))

生成的图应该类似于下图。您可以看到我们或多或少地检索了用于创建数据的原始参数，并且您还看到我们获得了一个围绕最高后验值的特定区域，这些区域也有一些数据支持，这是贝叶斯相当于置信度的间隔。

img

图：上排显示斜率（a）的后验估计，截距（b）和误差的标准偏差（sd）。下一行显示马尔可夫链参数值。

如有帮助，请采纳，十分感谢！

和这个是不是一样的 https://blog.csdn.net/zb1165048017/article/details/51778694

实际上在你的上一个问题中已经有人回答了一个合理答案😅，这是用的高斯过程的抽样，就是Gaussian的抽样。高斯过程的抽样具体可参考https://blog.csdn.net/weixin_44618906/article/details/113809458

Tristan Ursell
2D Random Number Generator for a Given Discrete Distribution
March 2012
[x0,y0]=pinky(Xin,Yin,dist_in,varargin);
'Xin' is a vector specifying the equally spaced values along the x-axis.
'Yin' is a vector specifying the equally spaced values along the y-axis.
'dist_in' (dist_in > 0) is a matrix with dimensions length(Yin) x length(Xin), whose values specify a 2D discrete probability distribution. The distribution does not need to be normalized.
 
'res' (res > 1) is a multiplicative increase in the resolution of chosen random number, using cubic spline interpolation of the values in 'dist_in'. Using the 'res' option can significantly slow down the code, due to the computational costs of interpolation, but allows one to generate more continuous values from the distribution.
 
[x0,y0] is the output doublet of random numbers consistent with dist_in.
 
The 'gendist' function required by this script, is included in this m-file.
 
Example with an anisotropic Gaussian at native resolution:
 
Xin=-10:0.1:10;
Yin=-5:0.1:5;
 
Xmat = ones(length(Yin),1)*Xin;
Ymat = Yin'*ones(1,length(Xin));
 
D1=1;
D2=3;
Dist=exp(-Ymat.^2/(2*D1^2)-Xmat.^2/(2*D2^2));
 
N=10000;
 
vals=zeros(2,N);
for i=1:N
[vals(1,i),vals(2,i)]=pinky(Xin,Yin,Dist);
end
 
figure;
hold on
imagesc(Xin,Yin,Dist)
colormap(gray)
plot(vals(1,:),vals(2,:),'r.')
xlabel('Xin')
ylabel('Yin')
axis equal tight
box on
 
Example with multiple peaks at 10X resolution:
 
Xin=-10:0.1:10;
Yin=-5:0.1:5;
 
Xmat = ones(length(Yin),1)*Xin;
Ymat = Yin'*ones(1,length(Xin));
 
D1=0.5;
D2=1;
Dist=exp(-(Ymat-3).^2/(4*D2^2)-(Xmat+6).^2/(2*D1^2))+exp(-(Ymat+2).^2/(2*D1^2)-(Xmat+1).^2/(2*D2^2))+exp(-(Ymat-1).^2/(2*D1^2)-(Xmat-2).^2/(2*D2^2));
 
N=10000;
res=10;
vals=zeros(2,N);
for i=1:N
[vals(1,i),vals(2,i)]=pinky(Xin,Yin,Dist,res);
end
 
figure;
hold on
imagesc(Xin,Yin,Dist)
colormap(gray)
plot(vals(1,:),vals(2,:),'r.')
xlabel('Xin')
ylabel('Yin')
axis equal tight
box on

可以参考一下

拿参加数模的经验和上次帮一个人解决的经验来说…也是问的二维分布的问题。。
实现二维分布需要以下几步
第一步，基于实验获得一条数据，要确定是否服从分布，将其转换为直方图，
第二步，从绘制的直方图中可以看出，该数据近似正态分布。从原始数据计算出的正态分布的关键值为：平均值，方差和标准差等等，
第三步，f =@（x）normpdf（x，均值，标准差），总概率密度f *等于拟合频数，然后拟合图。
第四步，如果要绘制频率概率图，则该方法与此类似，只需将频率除以总数即可得到频率，而正态分布f（x）为概率。
望采纳谢谢啦

感觉有点复杂