如何用倒数法生成数
How to generate number by inversion method
我有一个计算范围内数字的算法,例如
- 计算 F(x),x 从 0 到 n
- 生成统一数z
If(F(x)<=z<F(x+1))
然后 k=x
其中,F(x)是由Binomial_distribution计算得到的累积分布函数。
例如,F(x)=0.12, 0.2, 0.5, ..1.
z=0.3 然后 return k=1
我是用matlab实现的。然而,它 returned 错误的结果。你能看到我的实施对我有帮助吗?如有错误,请给出解决方案
clear all;
number_generated = 500000;
p=0.3;
n=20
%% Compute F(x): x=0..20
sum_F=0;
for x=0:n
sum_F=sum_F+nchoosek(n,x)*p^x*(1-p)^(n-x);
F(x+1)=sum_F;% matlab index from 1
end
k_arr=[]; %% Store k
for num=1:number_generated
%% Generate z_i uniform in the interval (0,1)
z=rand();
%% Find k such that F(x)<=z<F(x+1)
for i=1:length(F)-1
if(F(i)<=z & z<F(i+1))
k=i-1; % matlab index from 1
break;
end
end
k_arr=[k_arr k]; %% Record k in array
end
您的代码基本上是正确的,但是您忘记了 0
和 F(1)
之间的第一个间隔。
对您的代码的其他评论:
总是预分配数组
在逻辑条件下使用&&
代替&
。
更新代码:
clear all;
close all;
number_generated = 500000;
p=0.3;
n=20;
% Preallocate F array.
F = NaN(n + 2, 1);
% Set first value to 0
F(1) = 0;
% Save binomial distribution for plotting.
bin = zeros(n + 1, 1);
% Compute F(x): x=0..20
sum_F=0;
for x=0:n
% Save binomial distribution for plotting.
bin(x + 1) = nchoosek(n,x)*p^x*(1-p)^(n-x);
sum_F=sum_F+bin(x + 1);
% This is now x+2 because F(1) is 0.
F(x+2)=sum_F;% matlab index from 1
end
% Preallocate k_arr array.
k_arr=NaN(number_generated, 1); %% Store k
for num=1:number_generated
% Generate z_i uniform in the interval (0,1)
z=rand();
% Find k such that F(x)<=z<F(x+1)
for i=1:length(F)-1
if(z >= F(i) && z < F(i+1))
k=i-1; % matlab index from 1
break;
end
end
k_arr(num) = k; % Record k in array
end
% Plot expected result
figure
stairs((0:n) - 0.5, bin);
xlim([-1 20]);
% Plot sampled result
figure
histogram(k_arr);
xlim([-1 20]);
预期分布:
样本分布:
我有一个计算范围内数字的算法,例如
- 计算 F(x),x 从 0 到 n
- 生成统一数z
If(F(x)<=z<F(x+1))
然后k=x
其中,F(x)是由Binomial_distribution计算得到的累积分布函数。
例如,F(x)=0.12, 0.2, 0.5, ..1.
z=0.3 然后 return k=1
我是用matlab实现的。然而,它 returned 错误的结果。你能看到我的实施对我有帮助吗?如有错误,请给出解决方案
clear all;
number_generated = 500000;
p=0.3;
n=20
%% Compute F(x): x=0..20
sum_F=0;
for x=0:n
sum_F=sum_F+nchoosek(n,x)*p^x*(1-p)^(n-x);
F(x+1)=sum_F;% matlab index from 1
end
k_arr=[]; %% Store k
for num=1:number_generated
%% Generate z_i uniform in the interval (0,1)
z=rand();
%% Find k such that F(x)<=z<F(x+1)
for i=1:length(F)-1
if(F(i)<=z & z<F(i+1))
k=i-1; % matlab index from 1
break;
end
end
k_arr=[k_arr k]; %% Record k in array
end
您的代码基本上是正确的,但是您忘记了 0
和 F(1)
之间的第一个间隔。
对您的代码的其他评论:
总是预分配数组
在逻辑条件下使用
&&
代替&
。
更新代码:
clear all;
close all;
number_generated = 500000;
p=0.3;
n=20;
% Preallocate F array.
F = NaN(n + 2, 1);
% Set first value to 0
F(1) = 0;
% Save binomial distribution for plotting.
bin = zeros(n + 1, 1);
% Compute F(x): x=0..20
sum_F=0;
for x=0:n
% Save binomial distribution for plotting.
bin(x + 1) = nchoosek(n,x)*p^x*(1-p)^(n-x);
sum_F=sum_F+bin(x + 1);
% This is now x+2 because F(1) is 0.
F(x+2)=sum_F;% matlab index from 1
end
% Preallocate k_arr array.
k_arr=NaN(number_generated, 1); %% Store k
for num=1:number_generated
% Generate z_i uniform in the interval (0,1)
z=rand();
% Find k such that F(x)<=z<F(x+1)
for i=1:length(F)-1
if(z >= F(i) && z < F(i+1))
k=i-1; % matlab index from 1
break;
end
end
k_arr(num) = k; % Record k in array
end
% Plot expected result
figure
stairs((0:n) - 0.5, bin);
xlim([-1 20]);
% Plot sampled result
figure
histogram(k_arr);
xlim([-1 20]);
预期分布:
样本分布: