从 Matlab 中的任意离散概率密度函数生成随机样本
Generate random samples from arbitrary discrete probability density function in Matlab
我在 Matlab 中将任意概率密度函数离散为矩阵,这意味着对于每一对 x,y,概率存储在矩阵中:
A(x,y) = 概率
这是一个 100x100 矩阵,我希望能够从该矩阵中生成二维 (x,y) 的随机样本,并且如果可能的话,还能够计算均值和其他矩的PDF。我想这样做是因为在重采样之后,我想将样本拟合到近似的高斯混合模型。
我找遍了所有地方,但没有找到像这样具体的东西。我希望你能帮助我。
谢谢。
我不相信 matlab 具有生成具有任意分布的多元随机变量的内置功能。事实上,单变量随机数也是如此。但是虽然后者可以很容易地基于累积分布函数生成,但多元分布不存在 CDF,因此生成此类数字要混乱得多(主要问题是 2 个或更多变量具有相关性这一事实)。所以你的这部分问题远远超出了本站的范围。
由于只有一半的答案总比没有答案好,下面介绍了如何使用 matlab 以数字方式计算平均矩和更高矩:
%generate some dummy input
xv=linspace(-50,50,101);
yv=linspace(-30,30,100);
[x y]=meshgrid(xv,yv);
%define a discretized two-hump Gaussian distribution
A=floor(15*exp(-((x-10).^2+y.^2)/100)+15*exp(-((x+25).^2+y.^2)/100));
A=A/sum(A(:)); %normalized to sum to 1
%plot it if you like
%figure;
%surf(x,y,A)
%actual half-answer starts here
%get normalized pdf
weight=trapz(xv,trapz(yv,A));
A=A/weight; %A normalized to 1 according to trapz^2
%mean
mean_x=trapz(xv,trapz(yv,A.*x));
mean_y=trapz(xv,trapz(yv,A.*y));
因此,关键是您可以使用对 trapz
的两次连续调用对矩形网格执行二重积分。这允许您计算与网格具有相同形状的任何数量的积分,但缺点是矢量分量必须独立计算。如果你只想计算可以用 x
和 y
参数化的东西(它们自然与你的网格大小相同),那么你可以不用做任何额外的思考。
您还可以为集成定义一个函数:
function res=trapz2(xv,yv,A,arg)
if ~isscalar(arg) && any(size(arg)~=size(A))
error('Size of A and var must be the same!')
end
res=trapz(xv,trapz(yv,A.*arg));
end
这样你就可以计算像
这样的东西
weight=trapz2(xv,yv,A,1);
mean_x=trapz2(xv,yv,A,x);
注意:我在示例中使用 101x100 网格的原因是对 trapz
的双重调用应该以正确的顺序执行。如果您在调用中交换 xv
和 yv
,由于与 A
的定义不一致,您会得到错误的答案,但如果 A
是正方形,这将不明显.我建议在开发阶段避免对称数量。
如果你真的有一个由 A
定义的离散概率密度函数(与仅由 A
描述的连续概率密度函数相反),你可以 "cheat" 通过将您的二维问题转化为一维问题。
%define the possible values for the (x,y) pair
row_vals = [1:size(A,1)]'*ones(1,size(A,2)); %all x values
col_vals = ones(size(A,1),1)*[1:size(A,2)]; %all y values
%convert your 2D problem into a 1D problem
A = A(:);
row_vals = row_vals(:);
col_vals = col_vals(:);
%calculate your fake 1D CDF, assumes sum(A(:))==1
CDF = cumsum(A); %remember, first term out of of cumsum is not zero
%because of the operation we're doing below (interp1 followed by ceil)
%we need the CDF to start at zero
CDF = [0; CDF(:)];
%generate random values
N_vals = 1000; %give me 1000 values
rand_vals = rand(N_vals,1); %spans zero to one
%look into CDF to see which index the rand val corresponds to
out_val = interp1(CDF,[0:1/(length(CDF)-1):1],rand_vals); %spans zero to one
ind = ceil(out_val*length(A));
%using the inds, you can lookup each pair of values
xy_values = [row_vals(ind) col_vals(ind)];
希望对您有所帮助!
芯片
我在 Matlab 中将任意概率密度函数离散为矩阵,这意味着对于每一对 x,y,概率存储在矩阵中: A(x,y) = 概率
这是一个 100x100 矩阵,我希望能够从该矩阵中生成二维 (x,y) 的随机样本,并且如果可能的话,还能够计算均值和其他矩的PDF。我想这样做是因为在重采样之后,我想将样本拟合到近似的高斯混合模型。
我找遍了所有地方,但没有找到像这样具体的东西。我希望你能帮助我。
谢谢。
我不相信 matlab 具有生成具有任意分布的多元随机变量的内置功能。事实上,单变量随机数也是如此。但是虽然后者可以很容易地基于累积分布函数生成,但多元分布不存在 CDF,因此生成此类数字要混乱得多(主要问题是 2 个或更多变量具有相关性这一事实)。所以你的这部分问题远远超出了本站的范围。
由于只有一半的答案总比没有答案好,下面介绍了如何使用 matlab 以数字方式计算平均矩和更高矩:
%generate some dummy input
xv=linspace(-50,50,101);
yv=linspace(-30,30,100);
[x y]=meshgrid(xv,yv);
%define a discretized two-hump Gaussian distribution
A=floor(15*exp(-((x-10).^2+y.^2)/100)+15*exp(-((x+25).^2+y.^2)/100));
A=A/sum(A(:)); %normalized to sum to 1
%plot it if you like
%figure;
%surf(x,y,A)
%actual half-answer starts here
%get normalized pdf
weight=trapz(xv,trapz(yv,A));
A=A/weight; %A normalized to 1 according to trapz^2
%mean
mean_x=trapz(xv,trapz(yv,A.*x));
mean_y=trapz(xv,trapz(yv,A.*y));
因此,关键是您可以使用对 trapz
的两次连续调用对矩形网格执行二重积分。这允许您计算与网格具有相同形状的任何数量的积分,但缺点是矢量分量必须独立计算。如果你只想计算可以用 x
和 y
参数化的东西(它们自然与你的网格大小相同),那么你可以不用做任何额外的思考。
您还可以为集成定义一个函数:
function res=trapz2(xv,yv,A,arg)
if ~isscalar(arg) && any(size(arg)~=size(A))
error('Size of A and var must be the same!')
end
res=trapz(xv,trapz(yv,A.*arg));
end
这样你就可以计算像
这样的东西weight=trapz2(xv,yv,A,1);
mean_x=trapz2(xv,yv,A,x);
注意:我在示例中使用 101x100 网格的原因是对 trapz
的双重调用应该以正确的顺序执行。如果您在调用中交换 xv
和 yv
,由于与 A
的定义不一致,您会得到错误的答案,但如果 A
是正方形,这将不明显.我建议在开发阶段避免对称数量。
如果你真的有一个由 A
定义的离散概率密度函数(与仅由 A
描述的连续概率密度函数相反),你可以 "cheat" 通过将您的二维问题转化为一维问题。
%define the possible values for the (x,y) pair
row_vals = [1:size(A,1)]'*ones(1,size(A,2)); %all x values
col_vals = ones(size(A,1),1)*[1:size(A,2)]; %all y values
%convert your 2D problem into a 1D problem
A = A(:);
row_vals = row_vals(:);
col_vals = col_vals(:);
%calculate your fake 1D CDF, assumes sum(A(:))==1
CDF = cumsum(A); %remember, first term out of of cumsum is not zero
%because of the operation we're doing below (interp1 followed by ceil)
%we need the CDF to start at zero
CDF = [0; CDF(:)];
%generate random values
N_vals = 1000; %give me 1000 values
rand_vals = rand(N_vals,1); %spans zero to one
%look into CDF to see which index the rand val corresponds to
out_val = interp1(CDF,[0:1/(length(CDF)-1):1],rand_vals); %spans zero to one
ind = ceil(out_val*length(A));
%using the inds, you can lookup each pair of values
xy_values = [row_vals(ind) col_vals(ind)];
希望对您有所帮助!
芯片