2016-03-06 211 views
3

我有根据确定的概率密度函数(截短拉普拉斯),以生成随机值:生成赋予了PDF随机值

enter image description here

Sigma_Phi = 1; 
B = 1/(1-exp(-sqrt(2)*pi/Sigma_Phi)); 
Phi = -60:0.001:60; 

for iter=1:length(Phi) 
    if Phi(iter) < pi && Phi(iter)>=0 
     P(iter) = ((B*exp(-abs(sqrt(2)*Phi(iter)/Sigma_Phi)))/(sqrt(2)*Sigma_Phi)); 
    elseif Phi(iter) >= -pi && Phi(iter)<=0 
     P(iter) = ((B*exp(-abs(sqrt(2)*Phi(iter)/Sigma_Phi)))/(sqrt(2)*Sigma_Phi)); 
    else 
     P(iter)=0; 
    end 
end 

然后,我有生成PDF以下的随机值,但我不知道该怎么做。

回答

3

下面是一个简单的,一般的解决方案使用inverse transform sampling

% for reproducibility 
rng(333) 

Sigma_Phi = 1; 
B = 1/(1-exp(-sqrt(2)*pi/Sigma_Phi)); 
Phi = -10:0.001:10; 

P = nan(length(Phi),1); 
for iter=1:length(Phi) 
    if Phi(iter) < pi && Phi(iter)>=0 
     P(iter) = ((B*exp(-abs(sqrt(2)*Phi(iter)/Sigma_Phi)))/(sqrt(2)*Sigma_Phi)); 
    elseif Phi(iter) >= -pi && Phi(iter)<=0 
     P(iter) = ((B*exp(-abs(sqrt(2)*Phi(iter)/Sigma_Phi)))/(sqrt(2)*Sigma_Phi)); 
    else 
     P(iter)=0; 
    end 
end 

% create cdf 
cdf   = cumtrapz(Phi, P); 

% keep only the unique values: needed for interpolation 
idx_mid = (Phi < pi) & (Phi >= -pi); 
cdf = cdf(idx_mid); 
Phi = Phi(idx_mid); 
P = P(idx_mid); 

% number of required random draws 
n = 1e4; 
% generate uniformly distributed random numbers from [0,1] 
r = rand(n,1); 

% generate random numbers from the desired pdf; inverse transform sampling 
laplrnd = interp1(cdf, Phi, r); 

% Verfication plot 
[f,x] = hist(laplrnd,100); 

bar(x,f/trapz(x,f)) 
hold on 
plot(Phi, P, 'red', 'Linewidth', 1.2) 
legend('histogram from random values', 'analytical pdf') 

enter image description here