机器学习系列：随机采样方法之接受-拒绝采样

接受-拒绝采样或舍选法Acceptance-Rejection Sampling是一种通用的蒙特卡洛采样方法由‌约翰·冯·诺依曼John von Neumann‌于‌20 世纪 40 年代中期‌约 1947-1951 年间在曼哈顿计划期间正式提出并系统化。‌‌该方法用于从难以直接采样的复杂目标分布中生成样本。其核心思想是利用一个易于采样的‌建议分布‌来包裹目标分布并通过概率判断接受或拒绝样本从而间接获得服从目标分布的样本集。一、接受-拒绝采样原理假设我们有一个目标概率密度函数为了对其进行采样我们可以选择一个已知的、易于采样的建议分布proposal distribution比如均匀分布或高斯分布。然后我们需要找到一个合适的缩放常数, 使得能够完全覆盖目标函数即满足对所有都成立。采样的具体步骤如下从建议分布中随机采样一个点计算接受概率生成一个均匀分布的随机数如果则接受这个样本否则拒绝这个样本重新从步骤1开始这个过程可以形象地理解为我们在曲线下方随机投点如果点落在曲线下方则接受这个样本否则拒绝这个样本。通过这种方式我们最终得到的样本将服从目标分布。 具体过程是这样的先通过分布产生一个随机数如下图中x轴上的坐标点)然后用区间上的均匀分布产生一个随机数这等效于在区间上得到一个服从均匀分布的随机数(如图y轴上的点)。 如果落在区间中(即在曲线下方)则接受这个随机数如果落在上(即在曲线上方)则拒绝这个随机数。二、接受-拒绝采样的证明假设目标分布为而建议分布为且存在常数使得对所有得满足。 这里要证明得是, 也就是被接受随机数的分布与目标分布一致。证明之前先明确几个事件与关系1接受事件2均匀分布的概率密度函数3随机变量是独立的。证明首先由条件分布计算公式计算联合概率密度, 根据联合分布和条件分布的关系由于与是独立的所以而是,给定时, 事件Accept的概率密度即于是联合概率密度为于是得证。注证明过程参考了Warren的博客三、两个例子例1.以均匀分布作为建议分布为正态分布设计一个采样程序。function mainForRS1() %用接受-拒绝采样方法生成服从正态分布的随机数 clear all clc %绘制目标分布概率密度函数 x-5:0.1:5; mu0;sigma1.0; yNormpdf(x,mu,sigma); plot(x,y) hold on %采样 samplesSamplingNorm(20000,mu,sigma); %绘制样本直方图 histogram(samples,NumBins,150,Normalization,pdf,EdgeColor,b) end function samplesSamplingNorm(numberSample,mu,sigma) %以均匀分布作为建议分布 %用接受-拒绝采样方法生成服从正态分布的随机数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %输入: % numberSample: 样本数量 % mu: 均值 % sigma 标准差 %输出 % samples: 样本集 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %确定分布的范围[a,b] bmu5*sigma; % 均匀分布的上限 amu-5*sigma; % 均匀分布的下限 % xa:0.1:b; % ynormpdf(x,mu,sigma); % plot(x,y) % % 样本集 samples zeros(1,numberSample); %存放样本 %一个合适的M M(b-a)/(sqrt(2*pi)*sigma); accepted0; for i1:numberSample while 1 %从辅助分布采样 xa(b-a)*rand(); % 计算接受概率 alphaNormpdf(x,mu,sigma)/(M*uniformpdf(x,a,b)); %生成一个[0,1]上均匀分布随机数 urand(); if ualpha %接受这个样本 accepted accepted 1; samples(accepted)x; break; % 接受后退出循环继续下一个样本的生成 end end end end function pdfuniformpdf(x,a,b) %区间[a,b]上的均匀分布(建议分布) if xa xb pdf1/(b-a); else pdf0; end end function yNormpdf(x,mu,sigma) %正态分布密度函数(目标分布) yexp(-1/2*(x-mu).*(x-mu)/(sigma*sigma)); yy./(sqrt(2*pi)*sigma); end运行效果图样本直方图例 2. 要采样的目标分布密度函数为这里我们用一个正态分布作为建议分布, 其中选取缩放常数, 则蓝色和红色的曲线图形如下特别地本例中的正态分布随机数的生成用例1的函数SamplingNorm获得。程序如下function mainForRS2() clear all clc %绘制目标分布概率密度函数与建议分布曲线 x-4:0.1:6; y1q(x); mu1.2;sigma1.0; M3; y2M*Normpdf(x, mu, sigma); plot(x,y1,x,y2) hold on % 样本集 numberSample1e5; %样本数 samples zeros(1,numberSample); %存放样本 %一个合适的M M3.0; accepted0; for i1:numberSample while 1 %从辅助分布采样 xSamplingNorm(1,mu,sigma); %xnormrnd(mu, sigma); % 计算接受概率 alphaq(x)/(M*Normpdf(x, mu, sigma)); %生成一个[0,1]上均匀分布随机数 urand(); if ualpha %接受这个样本 accepted accepted 1; samples(accepted)x; break; % 接受后退出循环继续下一个样本的生成 end end end histogram(samples,NumBins,150,Normalization,pdf,EdgeColor,b) end function yq(x) %目标分布概率密度函数 y(0.3*exp(-(x-0.3).^2)0.7*exp(-(x-2).^2/0.3))/1.2113; end function yNormpdf(x, mu, sigma) %建议分布概率密度函数正态分布 y1/(sqrt(2*pi)*sigma)*exp(-(x-mu).^2/(2*sigma*sigma)); end运行效果图样本直方图为四、接受-拒绝采样的不足接受-拒绝采样的不足主要体现在以下几方面1当目标分布与提议分布形状差异大、常数M取值过大时大量样本会被直接拒绝有效样本生成速度极慢。2在高维分布中很难找到合适的提议分布来紧密包裹目标分布接受概率会呈指数级下降。复杂分布下的实用性大幅降低。3M值的确定对使用者的经验要求很高。