随机过程基础(1)---一维随机变量标准分布函数、常见性质MATLAB实现

基础概念大部分在高中数学、概率论和离散数学等中已经涉及。
本文根据均匀分布随机数生成正态分布随机变量部分代码摘自《Random Signal Processing》，Shaila Dinkar Apte著，Taylor Francis Group出版。
正态分布

老办法实现：

sigma=1;
mean=0;
x=(-10:0.1:10);
a=(x-mean).*(x-mean);
y=(1/(sigma*sqrt(2*pi)))*exp(-(((x-mean).*(x-mean))/(2*sigma^2)));
plot(x,y,'r-','linewidth',2);xlabel('x');ylabel('pdf');title('Gaussian pdf');

用现成的函数normpdf()：

plot(x,normpdf(x,mean,sigma),'r-','linewidth',2);xlabel('x');ylabel('pdf');title('Gaussian pdf');

所得结果与上相同。

使用normcdf查看标准正态概率分布：

函数norminv(P,MU,SIGMA)，计算参数为MU SIGMA的正态概率分布函数对应概率为P处的反函数值。

误差函数erf(x)$$erf(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-u^2}du$$反函数为erfinv(x)；
此外还有补余误差函数Y=erfc(X) =1-erf(x)。与标准正态分布的关系是$$\varphi (x)=\frac{1}{2}erfc(-\frac{x}{\sqrt{2}})$$其反函数为erfcinv(x)。
瑞利分布
来源：两个独立且方差相等的正态分布随机变量的平方和的开方
概率密度 f X ( x ) = { 2 b 2 ( x − α ) e − ( x − α ) 2 / 2 b 2 , x > α 0 ,     x < α f_X(x)=\{\begin{array}{lr} {\frac{2}{b^2}(x-\alpha)e^{-(x-\alpha)^2/2b^2} ,x > \alpha}\\ 0,~~~x<\alpha \end{array} fX​(x)={b22​(x−α)e−(x−α)2/2b2,x>α0,   x<α​概率分布$$F_X(x)=P(X\le x)=1-e^{\frac{-(x-\alpha {)}^2}{2b^2}}x\ge \alpha$$$$F_X(x)=0x<\alpha$$matlab有现成函数raylpdf(X,B)求瑞利分布概率密度。

x=[0:0.01:5];
b=0.5;
p=raylpdf(x,b);
plot(x,p,'g','linewidth',2); title('Rayleigh density function pdf'); xlabel('value of x');ylabel('pdf');

概率分布函数：raylcdf(X,B)

p=raylcdf(x,b);
plot(x,p,'linewidth',2); title('Rayleigh CDF'); xlabel('value of x');
ylabel('CDF');

同理，指数分布通过老办法或现成的函数exppdf(X,MU)和expcdf(X,MU) F X ( x ) = { 1 − e − ( x − a ) / b     x ≥ a 0          x < a F_X(x) = \{\begin{array}{lr}1-e^{-(x-a)/b}~~~x \ge a\\0~~~~~~~~x < a\end{array} FX​(x)={1−e−(x−a)/b   x≥a0        x<a​ f X ( x ) = { 1 b e − ( x − a ) / b     x ≥ a 0          x < a f_X(x) = \{\begin{array}{lr}\frac{1}{b}e^{-(x-a)/b}~~~x \ge a\\0~~~~~~~~x < a\end{array} fX​(x)={b1​e−(x−a)/b   x≥a0        x<a​

plot(x,exppdf(x,0.25),'linewidth',2);xlabel('x');
 ylabel('pdf');title('exponential pdf');

plot(x,expcdf(x,0.25),'linewidth',2);xlabel('x'); 
ylabel('cdf');title('exponential cdf');

此外，matlab还有均匀分布概率密度unifpdf(X,A,B)，均匀分布概率分布unifcdf(X,A,B)等。
以上都有现成的函数，实现简单。接下来讨论离散随机变量。
二项分布
概率密度：$$f_X(x)=\sum _{k=0}^N\frac{N!}{(N-k)!k!}{P}^k(1-P{)}^{N-k}\delta (x-k)$$以下用P=0.2，N=6绘制二项分布相关函数。matlab函数用法参见：
matlab官方文档_binopdf
matlab官方文档_stem

x=0:20;
y=binopdf(x,20,0.25);
stem(x,y);title('Binomial pdf function for N=20 and P=0.25');xlabel('x');
ylabel('pdf value');
figure;
clear all;
x=0:20;
y=binocdf(x,20,0.25);
stem(x,y,':diamondr');title('Binomial CDF function for N=20 and P=0.25');xlabel('x');
ylabel('CDF value');

泊松分布
是二项分布实验次数趋近无穷时，P趋向于0，从而NP乘积等于常量b的情况。pdf和cdf写作$$f_X(x)=e^{-b}\sum _{k=0}^{\infty }\frac{b^k}{k!}\delta (x-k)$$$$F_X(x)=e^{-b}\sum _{k=0}^{\infty }\frac{b^k}{k!}u(x-k)$$matlab函数为poisspdf(x,lambda)：

clear all;
x=0:20;
y=poisspdf(x,0.7);
stem(x,y);title('Poisson pdf function for N=20 and lamda=0.7');xlabel('x');ylabel('pdf value');
figure;clear all;
x=0:20;
y=poisscdf(x,0.7);
stem(x,y,':diamondr');title('Poisson CDF function for N=20 and lamda=0.7');xlabel('x');ylabel('CDF value');

使用中心极限定理生成服从高斯分布的随机变量
随机生成10组元素服从均匀分布、居于01之间的100×100矩阵。将这10个随机矩阵相加，通过hist命令按照最大最小间隔分成10份统计落在每一份中的数目。可以看出得出的结果和高斯分布十分相似。如果将矩阵改为1000*1000，会更加近似为正态分布。

clear all;
x1=rand(100);
x2=rand(100);
x3=rand(100);
x4=rand(100);
x5=rand(100);
x6=rand(100);
x7=rand(100);
x8=rand(100);
x9=rand(100);
x10=rand(100);
y=(1/sqrt(10))*(x1+x2+x3+x4+x5+x6+x7+x8+x9+x10);
z=hist(y,10);
bar(z);title('plot of a sum variable');
xlabel("bin value");
ylabel("amplitude");

利用以上方法，可以简单地统计音乐信号中样本点的分布：
相关函数：matlab_fseek ;matlab_fread

fp=fopen('example.wav');
fseek(fp,1500000,-1);
a=fread(fp,2048);
subplot(2,1,1);plot(a);title('plot of voiced part of a signal');
xlabel('sample no.');ylabel('amplitude');
b=hist(a,100);
figure;
bar(b);
title('plot of a histogram for a piece of music');xlabel('bin value');ylabel('amplitude');

样本点分布的两个统计峰值说明样本点在这两个区间出现频次异常高，取何种区间都有两个峰，但是峰的位置差别很大。这是需要探讨的问题。

卡方检验
基础概念高中已涉及。高中时用于判断两随机变量是否相关。卡方检验的本意不是判断独立性，而是判断预测模型是否准确。用在判断独立性上，如果两者相关，就说明用两者独立的结论来模拟的结果与实际结果不符。
卡方检验的式子是：$${\chi }^2=\sum _{i=1}^N\frac{(g_i-f_i{)}^2}{f_i}$$其中f_i是预测模型得出值，g_i是实际值。容易由这个式子和独立假设推出高中常用的独立样本四格表，及常用公式$$\frac{n(ad-bc{)}^2}{(a+b)(c+d)(a+c)(b+d)}$$

K-S检验
给一组未知分布的观测值，F(x)是对应某标准分布的对应值(在例子中是落在某个窄带区间的频数，如果$$D=max|F(x)-G(x)|$$得到的结果如果比K-S表中对应的小，那么可以用已知标准分布拟合未知分布。
这里使用K-S检验法检验从上述信号得到分布如果用随机数生成的正态分布拟合的优劣。

fp=fopen('example.wav');
fseek(fp,1500000,-1);
a=fread(fp,2048);
subplot(2,1,1);plot(a);title('plot of voiced part of a signal');
xlabel('sample no.');ylabel('amplitude');
b=hist(a,100);
subplot(2,1,2);
bar(b);
title('plot of a histogram for a speech segment');xlabel('bin value');ylabel('amplitude');
x1=rand(100);
x2=rand(100);
x3=rand(100);
x4=rand(100);
x5=rand(100);
x6=rand(100);
x7=rand(100);
x8=rand(100);
x9=rand(100);
x10=rand(100);
x11=rand(100);
x12=rand(100);
y=(1/sqrt(12))*(x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12);
z=hist(y,100);
figure;
bar(z);title('plot of a sum variable generated using Central limit theorem');xlabel('bin value');ylabel('amplitude');
sum=0.0;
for i=1:100,
d(i)=abs(b(i)-z(i));
sum=sum+d(i);
end
disp(sum);

得到的距离远大于K-S表中的对应值，说明用这种正态分布拟合该段信号的分布特别不可信。所以不能用正态分布拟合这段信号。