搜狗做题网在matlab中看瑞利分布
来源：学生作业帮 编辑：搜狗做题网作业帮 分类：数学作业 时间：2024/07/23 10:03:25
在matlab中看瑞利分布
有此程序
画瑞利分布的两个函数：
其一：
N=10000;
x=randn(1,N);
y=randn(1,N);
r=sqrt(0.5*(x.^2+y.^2));
step=0.1;
range=0:step:3;
h=hist(r,range);
fr_approx=h/(step*sum(h));
fr=(range/0.5).*exp(-range.^2);
plot(range,fr_approx,'ko',range,fr,'k');
grid;
其二：
a=randn(1,10000);
b=randn(1,10000);
c=abs(a+j*b);
step=0.1;
x=0:step:5.5;
h=hist(c,x);
pdfc=h/(step*sum(h));
plot(x,pdfc,'ko');
程序一有r=sqrt(0.5*(x.^2+y.^2));程序二有c=abs(a+j*b);由于模的包络服从瑞利分布，对比两程序，程序二显然是在求模，而程序一为什么多了乘以0.若去掉0.画出的点图与概率密度曲线是不重合的？
请主要解释下为什么求模时乘了0.才与密度函数重合
有此程序
画瑞利分布的两个函数：
其一：
N=10000;
x=randn(1,N);
y=randn(1,N);
r=sqrt(0.5*(x.^2+y.^2));
step=0.1;
range=0:step:3;
h=hist(r,range);
fr_approx=h/(step*sum(h));
fr=(range/0.5).*exp(-range.^2);
plot(range,fr_approx,'ko',range,fr,'k');
grid;
其二：
a=randn(1,10000);
b=randn(1,10000);
c=abs(a+j*b);
step=0.1;
x=0:step:5.5;
h=hist(c,x);
pdfc=h/(step*sum(h));
plot(x,pdfc,'ko');
程序一有r=sqrt(0.5*(x.^2+y.^2));程序二有c=abs(a+j*b);由于模的包络服从瑞利分布，对比两程序，程序二显然是在求模，而程序一为什么多了乘以0.若去掉0.画出的点图与概率密度曲线是不重合的？
请主要解释下为什么求模时乘了0.才与密度函数重合
% by dynamic
% 2009-8-4 9:49:08
%
Matlab中提供了直接的函数,生成瑞利分布随机数raylrnd
下面我提供一个直接的函数吧
function x = RelayDist(x0,sigma,n)
x = zeros(n,1);
for i=1:n
r = MixMOD(x0,10,1);
k = 0;
while r(10) == 0
k = k + 1;
r(10) = power(2,k);
r = MixMOD(r(10),2,1);
end
y = -2*log(r(10));
x(i) = sigma*sqrt(y);
x0 = x(i);
end
function r = MixMOD(x0,n,type)
format long;
M1 = power(2,31);
M2 = power(2,35);
a1 = 314159269;
a2 = power(5,15);
c1 = 453806245;
c2 = 1;
r = zeros(n,1);
x = zeros(n+1,1);
x(1) = x0;
if type == 1
for i=2:n+1
y = a1*x(i-1)+c1;
x(i) = mod(y,M1);
r(i-1) = x(i)/M1;
end
else
for i=2:n+1
y = a2*x(i-1)+c2;
x(i) = mod(y,M2);
r(i-1) = x(i)/M2;
end
end
format short;
% 2009-8-4 9:49:08
%
Matlab中提供了直接的函数,生成瑞利分布随机数raylrnd
下面我提供一个直接的函数吧
function x = RelayDist(x0,sigma,n)
x = zeros(n,1);
for i=1:n
r = MixMOD(x0,10,1);
k = 0;
while r(10) == 0
k = k + 1;
r(10) = power(2,k);
r = MixMOD(r(10),2,1);
end
y = -2*log(r(10));
x(i) = sigma*sqrt(y);
x0 = x(i);
end
function r = MixMOD(x0,n,type)
format long;
M1 = power(2,31);
M2 = power(2,35);
a1 = 314159269;
a2 = power(5,15);
c1 = 453806245;
c2 = 1;
r = zeros(n,1);
x = zeros(n+1,1);
x(1) = x0;
if type == 1
for i=2:n+1
y = a1*x(i-1)+c1;
x(i) = mod(y,M1);
r(i-1) = x(i)/M1;
end
else
for i=2:n+1
y = a2*x(i-1)+c2;
x(i) = mod(y,M2);
r(i-1) = x(i)/M2;
end
end
format short;