matlab在各个学科中的应用

MATLAB在各学科中的运用

MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。[1]

MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连 接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。

MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。

学习matlab后，研究电路及自动控制系统都非常直观方便。下面就matlab在几个学科中的应用举例：

应用一 Matlab在电路中的应用 应用二 Matlab在自动控制理论中的运用 应用三 基于Matlab的通信系统仿真 应用四 Matlab在金融工程中的运用 总结

应用一 MATLAB在电路中的应用

在大二上学期，我们电气工程及其自动化专业学习了电路这门课，下面引用matlab在电路里面的应用

MATLAB在直流稳态电路中的分析及应用

设计分析

1.运用MATLAB解决数值线性代数问题及MATLAB的实现；MATLAB在“电路工作原理”中的应用；MATLAB工具箱的运用。要求选其中的一道作为课设的题目，学会简单运用MATLAB。

基本电路是电类专业非常重要的专业基本课,其中,线性含源一端口的戴维南定理及正弦稳态电路的分析计算是大家普遍反映难于理解的内容。本文以基本电路理论中典型的直流电阻电路和含有复数运算的正弦稳态电路的分析计算为例，详述了如何分别运用MATLAB语言编程的方法来对电路进行仿真分析和计算。结论表明，应用这两种方法可以使复杂电路的分析和计算变得非常快捷、方便，从而为电路分析提供了一个有效的辅助工具。

图2-1是由电压源和电阻组成的简单电路，运用回路电流法，网孔法，节点法等一些经典的电路分析法即可解决此问题。电路也可用simulink进行仿真实验，并通过数据测量等从而检测计算及编程的结果。

2.1 基本电路数学模型分析

基本电路分析的基本方法是先建立数学模型，一般是电路方程组。然后通]]]过求解方程组，得到各支路电压和电流。

图2-1 基本电路图

如图2-1所示：设三个回路的电流分别为im1，i m2，im3。对图2-1应用 网孔电流法，可列出如下方程：

网孔1： -V1 +R1im1+ R3（im1-im2） =0 网孔2： -V2+R2im2+R4im2+R3（im2-im1）=0 经化简可得：（R1+R3）im1-R3im2=V1-V2 -R3im1+（R2+R3+R4）im2=V2

整理以上方程，并写成形如AI=B的矩阵方程形式，可得：

R1+R3 -R3 im1 V1-V2 =

-R3 R2+R3+R4 im2 V2

依题，要求I1、I2、I3，根据节点电流的关系，有：

I1=im1 ； I2=im2 ； I3 =im1-im2

由以上网孔电流分析法的分析过程可得以下程序流程图：

开始 设参数，并赋予参数值（与题中所给的要求一致） 为系数矩阵中个元素赋值，并得出系数矩阵，为下面计算作准备 由AI＝B得出I矩阵的计算公式，并列出I1、I2、I3的表达式 结束

图2-2 程序流程图

2.2 Simulink 仿真分析

用simulink仿真以上电路，并观察测量数据。I1、I2、I3的测量结果再和程序的结果相比较，若相同，则说明电路仿真和程序编写多是正确的；若不同，则说明电路的仿真与程序编写至少一个有问题，认真检查程序及电路的仿真，找出错误，认真分析，纠正错误，在比较结果。

使用simulink进行仿真一般分为两步：用户首先需要在仿真模型编辑窗口中搭建好自己的模型，设置好具体模型参数和仿真参数；然后用户就可以开始仿真，simulink将根据用户搭建的模型，模型系统在用户设定条件下的具体行为。对于建模，simulink提供了一个图形化的用户的界面（GUI），用户可以像用铅笔在纸上画图一样画模型图。simulink的所有模型是分级的，因此可以通过自上而下或者自下而上来建立模型。

编程级仿真

3.1程序编写及模块分析

程序如下： clear;

R1=5;R2=6;R3=10;R4=4;V1=15;V2=10; %给定初始值 A=[R1+R3 -R3;-R3 R2+R3+R4]; %给出系数矩阵A B=[V1-V2;V2]; %给出系数矩阵B I=A\\B; %求解未知变量矩阵I im1=I(1); %网孔电流的定义 im2=I(2);

I1=im1; %求解电流I1 I2=im2; %求解电流I2 I3=im1-im2; %求解电流I3 运行结果如图3-1所示：

图3-1 M文件建模仿真结果

程序分析：程序由2.1节的设计思路分析所得：开始先给元件（与电路图相对应）赋值，因为MATLAB编程计算矩阵的，所以此赋值是必须的，其为下面的的系数矩阵赋值作准备。由以上的元件赋值可得形如AI=B矩阵方程形式的系数矩阵，这样使得A，B，I矩阵变为已知。MATLAB提供了两种除法运算：左除（＼）和右除（／），一般情况下，I=A\\B是方程I*X=B的解，而I=B/A是方程I*A=B的解。在传统的MATLAB算法中，右除是先计算矩阵的逆再相乘，而左除则不需要计算矩阵直接进行除运算。通常右除快一点，但左除可避免被除矩阵的奇异性所带来的麻烦。由上所述，由公式I=A\\B 可求得I。

以上所得I矩阵的元素有：Im1、Im2。这四个元素在矩阵I中的顺序也如此。可见I(1）＝Im1，I（2）＝Im3，由第二节的分析知I1=im1；I2=im2；I3 =im1-im2 。 程序运行后，得到I1、I2、I3的值，并显示出来。与题中要求所求的一样，程序编写正确，并正常运行，与预期一样。

3.2用Simulink 仿真电路

总电路如图3-2所示：

图3-2 Simulink 仿真电路及显示图

图3-3所示为simulink仿真电路图的电流表A3值在示波器中的演示：

图3-3 simulink仿真电路图的电流表A3值

仿真电路图说明：图3-2与实际的电路图对应，下面对以上进行简单说明：A1为测量电流I1的器件，其输出电流I1，通过display1显示；A2为测量电流I2的器件，其输出电流I2，通过display2显示；A3为测量电流I3的器件，其输出电流I3，通过Scope显示。电路其余部分与实际电路差不多，值设置也一样。

仿真结果分析：由以上电路仿真可见，电路的编程求解及仿真所得结果一样。都达到了预期的效果，可见编程及电路仿真没有问题。MATLAB与simulink的结合运用是电路求解问题得以简化，使电路求解问题省时且精炼，精简。

应用二 matlab在自动控制理论中的运用

系统频率特性的测量与分析

模拟电路图及系统结构图分别如图4－l和图4－2。

图4－1

图4－2

取R3=500kΩ，则系统传递函数为

若输入信号U1(t)=U1sinωt，则在稳态时，其输出信号为

U2(t)=U2sin(ωt+Ψ)。

改变输入信号角频率ω值，便可测得二组U2／U1和Ψ随ω变化的数值，这个变化规律就是系统的幅频特性和相频特性。 实验2．考虑二阶系统传递函数模型

试用MATLAB 绘制出不同δ和ωn 的伯德图。 实验3. 已知系统的开环传递函数为

(1) 绘制根轨迹，并根据根轨迹图求若要使系统稳定，k 的最大值；

(2) 取k=10，绘制Bode 图，根据Bode 图求系统的幅值裕量和相角裕量；

取k=10，绘制Nyquist 图，并根据Nyquist 图判断系统的稳定性；检验用上述几 种方法判断系统稳定性的一致性。 三、步骤

1．连接被测量典型环节的模拟电路。电路的输入U1接A／D、D／A卡的DAl输出，电路的输出U2接A／D、D／A卡的ADI输入。检查无误后接通电源。

2．启动计算机，在桌面双击图标[自动控制实验系统]运行软件。

3．测试计算机与实验箱的通信是否正常，通信正常继续。如通信不正常查找原因使通信正常后才可以继续进行实验。

测频率图

4．选中[实验课题→系统频率特性测量→测频率图]菜单项，鼠标单击将弹出参数设置窗口。参数设置完成后点确认等待观察波形。

测波特图

5．在测量波特图的过程中首先应选择[实验课题→系统频率特性测量→测波特图数据采样]采集信息。

6．待数据采样结束后点击[实验课题→系统频率特性测量→测波特图→图形观测]即可以在显示区内显示出所测量的波特图。

测奈氏图

7．在测量波特图的过程中首先应选择[实验课题→系统频率特性测量→测奈氏图→数据采样]采集信息。

8．待数据采样结束后点击[实验课题→系统频率特性测量→测奈氏图→图形观测]即可以在显示区内显示出所测量的波特图。

9．实验2、实验3涉及的主要命令有：Bode（）、rlocus ()、rlocfind（）、margin()、nyquist() 。

为便于比较，可用hold on 指令将多条曲线放在一个图中。进一步，为清楚起见，用legend 指令在图中加注释。用rlocus ()画出根轨迹后，需要时用rlocfind（）找出临界稳定点。 实验一：

1、根据模拟电路图，可求出U2／U1和Ψ随ω变化的数值，见下表：

由此可以画出幅频特性和相频特性，如下图所示：

实验二

取ω=6，ξ分别为0.1、0.3、0.7、1.0、1.3，绘制其波特图，可得下图：

取ξ=0.7，ω分别为2、4、6、8、10，绘制其波特图，可得下图：

实验三

1.绘制根轨迹：由[k]=rlocfind()可求出K的最大值；

2.绘制BODE图和Nyquist图：

由图可以得到;

系统的幅值裕度h=-14dB,相角裕度γ=-45.8度。相应的截止频率WC=2.25，穿越频率wx=1.14..

有奈氏判据可得，系统闭环不稳定。

应用三 基于Matlab的通信系统仿真

系统综述

利用Matlab仿真软件，完成如图1所示的一个基本的数字通信系统。信号源产生0、1等概分布的随机信号，映射到16QAM的星座图上，同时一路信号已经被分成了实部和虚部，后边的处理建立在这两路信号的基础上。实部、虚部信号分别经过平方根升余弦滤波器，再加入高斯白噪声，然后通过匹配滤波器（平方根升余弦滤波器）。最后经过采样，判决，得到0、1信号，同原信号进行比较，给出16QAM数字系统的误码。

结构框图

图1 待构建系统的框图

系统实现 随机信号的生成

利用Matlab中自带的函数randsrc来产生0、1等概分布的随机信号。源代码如下所示：

%====定义待仿真序列的维数 N global N N=320;

%====定义产生‘1’的概率为 p global p p=0.5;

%============================== %首先产生随机二进制序列 source=randsrc(1,N,[1,0;p,1-p]);

0、1等概分布的随机信号如图2所示。

0 1等概分布的信号1.51 信号幅度0.50-0.505101520253035404550 图2 0、1等概分布的随机信号波形图

星座图映射

将等概分布的0、1信号映射到16QAM星座图上。每四个bit构成一个码子，具体实现的方法是，将输入的信号进行串并转换分成两路，分别叫做I路和Q路。再把每一路的信号分别按照两位格雷码的规则进行映射，这样实际上最终得到了四位格雷码。为了清楚说明，参看表1：

表1 两位格雷码的映射规律

两位0、1码 映射后（按格雷码） 0 0 0 1 1 1 1 0 源代码如下所示：

function [y1,y2]=Qam_modulation(x) %QAM_modulation

%============================== %对产生的二进制序列进行QAM调制

%=====首先进行串并转换，将原二进制序列转换成两路信号 N=length(x); a=1:2:N; y1=x(a); y2=x(a+1);

%=====分别对两路信号进行QPSK调制

%======对两路信号分别进行2－4电平变换 a=1:2:N/2; temp1=y1(a); temp2=y1(a+1);

y11=temp1*2+temp2;

temp1=y2(a); temp2=y2(a+1);

y22=temp1*2+temp2;

%=======对两路信号分别进行相位调制 a=1:N/4;

-3 -1 1 3 y1=(y11*2-1-4)*1.*cos(2*pi*a); y2=(y22*2-1-4)*1.*cos(2*pi*a); %========按照格雷码的规则进行映射 y1(find(y11==0))=-3; y1(find(y11==1))=-1; y1(find(y11==3))=1; y1(find(y11==2))=3;

y2(find(y22==0))=-3; y2(find(y22==1))=-1; y2(find(y22==3))=1; y2(find(y22==2))=3;

得到的星座图如图3所示，图上注明了每一个点对应的01序列。

16 QAM星座图54(1 1 1 1) 3(0 1 0 0) (0 1 1 0) (1 1 1 0) 2(0 1 0 1) (0 1 1 1) (1 1 1 0) (1 1 0 1) 10(1 0 1 1) -1(0 0 0 1) (0 0 1 1) (1 0 0 1) -2(0 0 0 0) (0 0 1 0) (1 0 1 0) -3(1 0 0 0) -4-5-5-4-3-2-1012345 图3 16QAM星座图

从上边的星座图上可以清楚的看到，任意相邻的两个点之间它们对应的4个bit中只有一个有差别，也就是格雷码的特点。而采用格雷码主要目的是当信噪比较大时，也就是系统的误码率比较低的情况下，当出现一个符号错误的情况下，往往只是这个符号中的一个bit位出现了误码，因此这个情况下误码率和误bit率是4：1，这一特性在后边的误码率计算的过程中会有应用。

插值

为了能够模拟高斯白噪声的宽频谱特性，以及为了能够显示波形生成器（平方根升余弦滤波器）的效果，所以在原始信号中间添加一些0点。具体实现是分别在信号的I路和Q路中，任意相邻的两个码字之间添加7个0。源代码如下所示：

function y=insert_value(x,ratio)

%=============================== ％x是待插值的序列，ratio是插值的比例。 %两路信号进行插值

％首先产生一个长度等于ratio倍原信号长度的零向量 y=zeros(1,ratio*length(x));

％再把原信号放在对应的位置 a=1:ratio:length(y); y(a)=x;

对I路和Q路信号进行插值后的波形图如图4所示。

两路信号的波形图420-2-4010203040实部信号50607080420-2-4010203040虚部信号50607080 图4 经过插值后的两路信号波形图

波形成形（平方根升余弦滤波器）

为了避免相邻传输信号之间的串扰，多元符号需要有合适的信号波形。图1中的方波是在本地数字信号处理时常见的波形，但在实际传输时这种方波并不合适。根据奈奎斯特第一准则，在实际通信系统中一般均使接收波形为升余弦滚降信号。这一过程由发送端的基带成形滤波器和接收端的匹配滤波器两个环节共同实现，因此每个环节均为平方根升余弦滚降滤波，两个环节合成就实现了一个升余弦滚降滤波。实现平方根升余弦滚降信号的过程称为“波形成形”，通过采用合适的滤波器对多元码流进行滤波实现，由于生成的是基带信号，因此这一过程又称“基带成形滤波”。

平方根升余弦滤波器的冲激响应

基带平方根升余弦滤波器具有以下定义的理论函数

H(f)?1?1?当|f|?f(1??)2?11?fN?|f|??H(f)??sin()? 当f(1??)?|f|?f(1??) ??222f?N???当|f|?f(1??)H(f)?0??其中：fN?R1?s是奈奎斯特平率，?是滚降系数。 2Ts2下面给出平方根升余弦滤波器的冲激响应曲线，如图5所示。

平方根升余弦滤波器的冲激响应0.60.50.4Amplitude0.30.20.10-0.10510152025n (samples)30354045 图5 平方根升余弦滤波器的冲激响应曲线

从上图上不难看出来，平方根升余弦滤波器的冲激响应很显然的引入了符号间干扰（ISI）即它的冲激响应在相邻的抽样点上的值并不象升余弦滤波器那样恒为0。然而造成这一后果的原因在于，当我们引入平方根升余弦滤波器的时候，就是认为整个信道，也就是说，包括信号发送端的滤波器和信号接收端的滤波器，总体的效果是避免了符号间干扰（ISI），所以，单独看这每一个滤波器，勿庸置疑，它们都是存在着符号间干扰（ISI）的。

经过平方根升余弦滤波器后的信号

I路和Q路信号经过平方根升余弦滤波器后，成形后的波形如图6所示。 源代码如下

%x1、x2是两路输入信号，fd是信号信息位的频率，fs是信号的采样频率 function [y1,y2]=rise_cos(x1,x2,fd,fs) %生成平方根升余弦滤波器 [yf, tf]=rcosine(fd,fs, 'fir/sqrt'); %对两路信号进行滤波

[y1, to1]=rcosflt(x1, fd,fs,'filter/Fs', yf); [y2, to2]=rcosflt(x2, fd,fs,'filter/Fs', yf);

10倍载波调制

将通过成形滤波器后的信号调制到10倍于原频率的载波上。由于在仿真的过程中，只

4 滤波后的信号原始信号 20-2-4010203040实部信号506070804 滤波后的信号原始信号 20-2-4010203040虚部信号50607080 图6 通过平方根升余弦滤波器后的两路信号

能用离散的点来模拟连续信号，因而为了能够显示出一个正弦曲线，至少需要在一个正弦周期内采样到4个以上的点，这里，我们在一个周期内采10个点。假设最初的0、1信号的频率是1Hz，那么I路和Q路符号传输的频率是1/4Hz，而10倍频是建立在I路或Q路符号频率的基础上，也就是说，载频的频率是2.5Hz。按照前面的假设，那么相邻两个采样点之间的时间间隔是0.04s。而一个完整周期内的正弦波形的幅值是相同的，都是对应的这个周期内的I路和Q路线性叠加，调制后的信号为，

y(t)?I(t)cos2?fct?Q(t)sin2?fct 其中，fc为载波频率。

源代码如下

%载波调制

%x1,x2代表两路输入信号，f是输入信号的频率，hf是载波的频率 function [t,y]=modulate_to_high(x1,x2,f,hf)

%产生两个中间变量，用来存储插值后的输入信号 yo1=zeros(1,length(x1)*hf/f*10); yo2=zeros(1,length(x2)*hf/f*10); n=1:length(yo1);

%对输入信号分别进行插值，相邻的两个点之间加入9个点，且这9个点的值同第0个点的值相同 yo1(n)=x1(floor((n-1)/(hf/f*10))+1); yo2(n)=x1(floor((n-1)/(hf/f*10))+1); %生成输出输出信号的时间向量 t=(1:length(yo1))/hf*f/10; %生成载波调制信号

y=yo1.*cos(2*pi*hf*t)-yo2.*sin(2*pi*hf*t);

得到的调制到载频的信号波形如图7和图8所示，其中图7主要为了显示一个脉宽周期内的调制信号波形，图8则是显示载波信号的整体情况。

载波调制信号波形21.510.5幅度0-0.5-1-1.5-2101112131415时间1617181920 图7 载波调制信号展开图

载波调制信号波形2.521.510.5幅度0-0.5-1-1.5-2-2.50100200300时间400500600700 图8 载波调制信号整体图

加入高斯白噪声

将通过成形滤波器后的信号送到具有高斯白噪声特征的加性信道中，相当于在原信号上加入高斯白噪声。由于高斯白噪声加在了通过插值和滤波后的点上，因此在计算信噪比的时候存在一个信噪比换算的问题。当我们把仿真得到的误码率曲线同理论的误码率曲线相比较的时候，两者的信噪比的定义必须是一致的。一致包括两个方面，一是二者均为每bit符号上的信号功率和噪声功率的比值，另一个是信号的功率是指那些信息点上的平均功率，噪声也是指信息点上所对应的噪声的平均功率，但由于噪声的功率谱密度是一个定值，所以噪声的平均功率实际上就是噪声的功率谱密度。对于第二点，由于所有信号的平均功率和信息点上的信号的平均功率不同，所以需要在加入高斯噪声的时候进行纠正，具体的公式推导如下。 设SNR是最后理论计算中的信噪比，SNR'是加入高斯白噪声后的整体信号（包括插值后的点）的信噪比,

Eb是每bit信息点的平均能量，Eb'是每bit信号的平均能量，No是噪

声的平均功率，现在需要推导出SNR与SNR'的关系。

?Eb?No?SNRSNR'Eb' ???Eb'SNREb??SNR'?No即两个信噪比的比值就是平均能量的比值。 源程序如下

%对输入的两路信号加高斯白噪声，返回处理后的两路信号，信息点等效bit信噪比为snr的值 function [y1,y2]=generate_noise(x1,x2,snr) %snr1代表snr对应的符号信噪比 snr1=snr+10*log10(4);

%算出所有信号的平均功率 ss=var(x1+i*x2,1); %加入高斯白噪声

y=awgn([x1+j*x2],snr1+10*log10(ss/10),'measured'); y1=real(y); y2=imag(y);

给出加入高斯白噪声的两路信号波形。

加入高斯白噪声后的波形42幅度0-2-4010203040实部信号5060708042幅度0-2-4010203040虚部信号50607080 图9 加入高斯白噪声的两路信号波形

匹配滤波器

在数字传输系统中，滤波器是不可缺少的。滤波器的一个作用是使基带信号频谱成形，例如为了满足奈奎斯特第一准则，基带信号频谱通常采用升余弦滚降形状，这一点在波形成形部分已经有了较详细的介绍。而滤波器的另一个重要作用是在接收端限制白噪声，将信号频带外的噪声滤掉，减少它对信号正确判决的影响。为了能够使滤波器输出信噪比在信息抽样时刻的信噪比最大，所以引入了匹配滤波器。

假设匹配滤波器的频率传递函数为H(f)，时域冲激响应为h(t)。滤波器输入为发送信号与噪声的叠加，即

x(t)?S(t)?n(t)

这里，S(t)为信号，它的频谱函数为S(f)。n(t)为白色高斯噪声，其双边功率谱密度为滤波器的输出为

n0。2y(t)?[S(t)?n(t)]?h(t)

其中信号部分为

?yS(t)?S(t)?h(t)???j2?ftS(f)H(f)edf ?在t?T时刻输出的信号抽样值为

?yS(T)????S(f)H(f)ej2?fTdf

滤波器输出噪声的功率谱密度为

Fn0(f)?Fn(f)|H(f)|2

平均功率为

?N0????Fn0(f)|H(f)|2df

因此，t?T时刻的输出信噪比为

?SNR?|?S(f)H(f)ej2?fTdf|???

??2F(f)|H(f)|dfn?匹配滤波器的传递函数使SNR达到最大。在这里利用Schwartz不等式求解，最后得到传递

函数的表达式为

H(f)?KS*(f)e?2?fT

即传递函数与信号频谱的复共轭成正比。 传递函数的时域响应为

h(t)?KS(T?t)

匹配滤波器的最大输出信噪比为

|S(f)|22ESNR??df?s

n0/2n0??其中，Es为观察间隔内的信号能量。

具体到这个通信系统中，由于信号的时域响应为

?S(t)??(T)?h0(t)?h0(T?t)

其中h0(t)是平方根升余弦滤波器的冲激响应。 结合上式可以得到

h(t)?Kh0(t)

匹配滤波器实质上是一个具有与发射端的基带成形滤波器相同的滚降系数的平方根升余弦滤波器。接收端的“匹配滤波”是针对发射端的成形滤波而言，与成形滤波相匹配实现了数字通信系统的最佳接收。它与基带成形滤波器共同构成了一个奈奎斯特滤波器。 源代码同平方根升余弦滤波器的源代码相同。 信号通过匹配滤波器后的波形如图10所示

接收端经过匹配滤波器后的波形42幅度0-2-4010203040实部信号5060708042幅度0-2-4010203040虚部信号50607080 图10

经过匹配滤波器后的波形

从上边的波形可以看出来，经过匹配滤波器后的信号明显很平滑，这正好反映了低通滤波器的特性，滤掉了高频分量，为了明显的反映这一特点，将一段高斯白噪声经过匹配滤波器。波形对比如图11所示。

高斯白噪声通过匹配滤波器3210-1-2-30102030405060708090100高斯白噪声波形3210-1-20102030405060经过滤波器后的波形708090100 图11 高斯白噪声经过匹配滤波器后的波形

经过仿真，发现高斯白噪声经过一个平方根滤波器后方差保持不变。因此在加入高斯噪声时给定的信噪比需要有一定的修正，即要保证在信息点上的信噪比为给定的值。我们知道

?SSNR?10lg()，当给定snr时，需要加入的高斯噪声的功率谱密度n0?1010?S，其

n0snr[(32?32)?(32?12)?(12?32)?(12?12)]中S??10，在具体使用AWGN函数时，snr

4值是建立在当前输入信号的平均功率S'的基础上的，所以

snr'?10lg(S'S'S'S')?10lg(snr)?snr?10lg()?snr?10lg()。

?n0S101010?S采样

由于从匹配滤波器出来的信号的点数8倍于原来信息的点数，为了恢复出原信号，所以需要对该信号进行采样。从匹配滤波器出来时，首先要剔除卷积过程中冗余的点，接着抽取现在信号中的第1个，第9个，??，第8×k＋1个点，源代码如下：

function [y1,y2]=pick_sig(x1,x2,ratio) y1=x1(ratio*3*2+1:ratio:length(x1)); y2=x2(ratio*3*2+1:ratio:length(x1));

将这时的数据画到星座图上。

QAM星座图543210-1-2-3-4-5-5-4-3-2-1012345 图12

信噪比为10db时的星座图

判决解调

经过前边的匹配滤波器解调或者称为相关解调产生了一组向量，在这里就是一个一维的向

量，根据最大后验概率（MAP）准则（由于各个信号的先验概率相等，所以页可以认为是最大似然准则），得到了最小距离检测。具体在本仿真系统中，判断为各个信号的门限如表2所示。判决后得到的数据再按照格雷码的规则还原成0、1信号，最终将两路0、1信号合成一路0、1信号，用来同最初的信号一起决定误码率。

表2 判决电平对应表

判决前的信号的幅度 对应的判决后的幅度 －3 －1 1 3 A??2 ?2?A?0 0?A?2 A?2 源代码如下

function y=demodulate_sig(x1,x2) %对x1路信号进行判决 xx1(find(x1>=2))=3;

xx1(find((x1<2)&(x1>=0)))=1; xx1(find((x1>=-2)&(x1<0)))=-1; xx1(find(x1<-2))=-3; %对x2路信号进行判决 xx2(find(x2>=2))=3;

xx2(find((x2<2)&(x2>=0)))=1; xx2(find((x2>=-2)&(x2<0)))=-1; xx2(find(x2<-2))=-3;

%将x1路信号按格雷码规则还原成0、1信号 temp1=zeros(1,length(xx1)*2); temp1(find(xx1==-1)*2)=1; temp1(find(xx1==1)*2-1)=1; temp1(find(xx1==1)*2)=1; temp1(find(xx1==3)*2-1)=1;

%将x2路信号按格雷码规则还原成0、1信号 temp2=zeros(1,length(xx2)*2); temp2(find(xx2==-1)*2)=1; temp2(find(xx2==1)*2-1)=1; temp2(find(xx2==1)*2)=1; temp2(find(xx2==3)*2-1)=1; %将两路0、1信号合成一路 y=zeros(1,length(temp1)*2); y(1:2:length(y))=temp1; y(2:2:length(y))=temp2;

误码率曲线

将解调后的数据同原始数据相比较，得到该信噪比下所对应的误码率。为了得到误码率曲线，需要得到在不同的信噪比下的误码率。在仿真的过程中，假设要得到一个值得信赖的误码率数据点，至少需要在最后的数据比较的过程中得到100个错误，那么参与仿真的数据点就应该是误码率的倒数乘以100，为了提高程序的效率，首先计算出某个信噪比对应的理论的误码率，然后估计出待仿真的点数。

对于16QAM信号星座图等效为在两个正交载波上的两个PAM信号，其中每一个具有4个信号点。因为在解调器中可以将相位正交的两个信号分量完全分开，所以QAM的错误概率可以由PAM的错误概率求得。16QAM系统的正确判决概率是

Pc?(1?P4)2

式中，P4是4元PAM的错误概率，在等效QAM系统的每一个正交信号中，4元PAM具有一半的平均功率，通过适当的修改4元PAM的错误概率，可以得到

13P4?2(1?)Q(SNRs)

415其中SNRs是平均符号SNR。因此，16QAM的错误概率是

Pe?1?(1?P4)2

具体的源代码如下：

clear;

%用来仿真QAM的误bit率 snr=1:1:11;

%先来计算理论误bit率

error_theory=(1-(1-(2*(1-1/sqrt(16))*1/2*erfc(1/sqrt(2)*sqrt(3*4*10.^(snr/10)/(16-1))))).^2)/4; %用理论的误bit率来决定需要仿真的点数 N=floor(1./error_theory)*1000+100; N(find(N<5000))=5000; %开始仿真 global p;

for i=1:length(N);

%首先产生随机二进制序列

source=randsrc(1,N(i),[1,0;p,1-p]);

%对产生的二进制序列进行QAM调制

[source1,source2]=Qam_modulation(source); %插值

sig_insert1=insert_value(source1,8); sig_insert2=insert_value(source2,8);

[source1,source2]=rise_cos(sig_insert1,sig_insert2,0.25,2); %====将滤波后的信号加入高斯白噪声

[x1,x2]=generate_noise(source1',source2',snr(i)); %[x1,x2]=generate_noise(source1,source2,snr(i)); sig_noise1=x1'; sig_noise2=x2';

[sig_noise1,sig_noise2]=rise_cos(sig_noise1,sig_noise2,0.25,2); [x1,x2]=pick_sig(sig_noise1,sig_noise2,8); sig_noise1=x1; sig_noise2=x2; %解调

signal=demodulate_sig(sig_noise1,sig_noise2); %计算误bit率

error_bit(i)=length(find(signal-source)~=0)/N(i); end;

%画出图形

semilogy(snr,error_bit,'-b'); hold on

semilogy(snr,error_theory,'-r')

误码率曲线图如图 所示。从图上可以看到当信噪比小的情况下，仿真曲线和理论曲线差距略大，而随着信噪比的增大，仿真曲线越来越逼进理论曲线。简单分析不难看出，由于理论误码率曲线是建立在误符号率除以4的基础上的，而这一条件的前提是出现误符号的时候，

一个符号中只有一个bit位发生了错误，这表明误码率比较低，也就是说明信噪比比较大。所以，当信噪比比较小的时候，理论计算的误码率的值要小于仿真得到的值。

10016QAM误码率曲线仿真曲线理论曲线10-1误码率10-210-310-4123456信噪比7891011 图13 误码率曲线图

整体程序构架

前面给出的分别是每一个模块对应的函数，下面的程序表示如何将上边的各个模块连接起来。

clear;

%====定义待仿真序列的维数 N global N N=320;

%====定义产生‘1’的概率为 p global p p=0.5;

%============================== %首先产生随机二进制序列 source=randsrc(1,N,[1,0;p,1-p]);

%============================== %对产生的二进制序列进行QAM调制 [source1,source2]=Qam_modulation(source); %=============================== %画出星座图 figure(1);

plot_astrology(source1,source2);

%=============================== %两路信号进行插值

sig_insert1=insert_value(source1,8); sig_insert2=insert_value(source2,8);

%=============================== %画出两路信号的波形图 figure(2);

plot_2way(sig_insert1,sig_insert2,length(sig_insert1),0.5); title('两路信号的波形图');

%=============================== %通过低通滤波器

[sig_rcos1,sig_rcos2]=rise_cos(sig_insert1,sig_insert2,0.25,2); %=============================== %画出两路信号的波形图

figure(3);

plot_2way(sig_rcos1,sig_rcos2,length(sig_rcos1)/4,0.5); title('通过低通滤波器后两路信号波形图');

%stem_2way(sig_insert1,sig_insert2,length(sig_insert1)/4,0.5); %=============================== %====将基带信号调制到高频上

[t,sig_modulate]=modulate_to_high(sig_rcos1,sig_rcos2,0.25,2.5); figure(4);

plot(t(1:500),sig_modulate(1:500));

%=============================== %====将滤波后的信号加入高斯白噪声 snr=10;

[x1,x2]=generate_noise(sig_rcos1,sig_rcos2,snr); sig_noise1=x1'; sig_noise2=x2'; end; figure(5)

plot_2way(sig_noise1,sig_noise2,length(sig_noise1)/4,0.5); %=============================== %====经过匹配滤波器

% [x1,x2]=match_flt(sig_noise1,sig_noise2,0.25,2); % sig_match1=x1'; % sig_match2=x2';

[sig_match1,sig_match2]=rise_cos(sig_noise1,sig_noise2,0.25,2); figure(6);

plot_2way(sig_match1,sig_match2,length(sig_match1)/4,0.5); %=============================== %采样

[x1,x2]=pick_sig(sig_match1,sig_match2,8); sig_pick1=x1; sig_pick2=x2; %画出星座图 figure(7)

plot_astrology(sig_pick1,sig_pick2); %解调

signal=demodulate_sig(sig_pick1,sig_pick2); %画出误码率曲线图 figure(8) plot_snr;

讨论 信噪比修正

前边提到了在加入高斯白噪声时，需要对信噪比（SNR）进行修正。由于接收滤波器是线性的，根据随即过程理论的知识，高斯随机过程的线性变换仍然是高斯随机过程，因此滤波器的输出噪声也是高斯的。下面来计算通过平方根升余弦滤波器后的高斯噪声的方差：

?N????N0|H(f)|2df 其中|H(f)|2是升余弦滤波器的频域响应。

?有升余弦滤波器的定义可以知道

???H0(f)df?1，所以N?N0。

也就是通过平方根升余弦滤波器后的高斯白噪声的平均功率保持不变。

仿真曲线

在仿真的过程中，首先不考虑两个滤波器的影响，直接在映射到星座图上的信号加入高

10016QAM误码率曲线仿真曲线理论曲线蒙特卡罗曲线10-1误码率10-210-310-4123456信噪比7891011 图14

3条曲线比较

斯白噪声，然后解调，同发送端进行比较，最后得到误码率曲线，在图 中用绿色表示。这条曲线同理论公式推导的条件一致，因此它应当同理论曲线（图中用红色表示）吻合。而如图 所示，在信噪比较小的情况下，二者的差距比较大。这是由于理论公式中近似的认为在采用格雷码编码的前提下，误符号率和误bit率之间满足4：1的关系，而在信噪比较小的情况下，一个符号可能错误会超过1个，这样误符号率和误bit率之间的比值将会小于4：1，所以理论计算的误bit率的值要偏小。这也说明，在后边加入两个滤波器的系统中，我们判断误码率曲线是否正确的依据应该是这条绿色的曲线，而不应该是“理论曲线”。 而加入滤波器的系统经过信噪比修正后得到的误码率曲线同蒙特卡罗曲线基本吻合，也说明了仿真的正确和引入修正的必要。

仿真终止条件

在做仿真的时候，考虑到不同的信噪比对应不同的误码率，因此需要仿真的点数也应该相应的变化。由于16QAM系统已经有了比较精确的理论误码率公式，所以可以先由给定的信噪比得到理论的误码率，然后再由这个误码率得到需要仿真的点数。在点数少于一定域值的情况下，直接用域值代替。这样在很大程度上节省了时间，同时也保证了仿真得到的误码率的可靠性。

升余弦滤波器

由于在本系统中直接采用了Matlab中communication工具箱中的rcosine和rcosflt两个函数来实现平方根升余弦滤波器，所以对这两个函数的一些属性也有必要讨论一下。 其中rcosine函数是用来产生数字滤波器的参数，而rcosflt函数调用rcosine函数产生的结果对序列进行滤波。这里需要注意的是，在rcosflt函数的输入参数中如不特别声明，它将默认的对输入信号进行过采样（Upsampling），过采样率由输入参数中的信号频率（Fd）和采样频率（Fs）决定。在经过成型滤波器（第一个平方根升余弦滤波器）时，需要对原信号进行过采样，这时可以通过rcosflt函数在完成波形生成的同时进行8倍过采样，当然也可以先进行过采样后，再将得到的结果送入到rcosflt函数中，同时声明不需要进行过采样。而在经过匹配滤波器的时候，rcosflt函数不要再过采样。

另外还有一点，在rcosine函数中，有一个参数是延时（Delay），这个参数的值就是平方根升余弦滤波器时域响应中边瓣的个数。默认值是3，同时滚降系数是0.5。当滚降系数不是0.5的情况下，延时（Delay）的值需要适当的变化。例如当滚降系数小于0.5时，Delay的值应该大于3，以保证能尽可能多的将旁瓣的效果加入的滤波器中，当然随着Delay的变大，数据处理的工作量也会增大。

应用四 Matlab在金融工程中的运用

MATLAB金融计算与金融数据处理 Matlab金融工具箱模块 1. Financial Toolbox

Matlab自带金融工具箱，具有下列功能：

固定收益计算 日期数据处理 资产均值－方差分析 时间序列分析 有价证卷的收益和价格 统计分析 定价和灵敏度分析 年金和现金流计算 抵押支持债卷

Financial Derivatives Toolbox 是金融衍生产品工具箱，用于固定收益、金融衍生物以及风险投资评估分析，也可用于各种金融衍生物定价策略以及敏感度分析。 2. Financial Derivatives Toolbox 3. Financial Time Series Toolbox

Financial Time Series Toolbox 用于分析金融市场的时间序列数据。金融数据是时间序列数据，例如股票价格或每天的利息波动，可以用该工具箱进行更加直观的数据管理。该工具箱支持下列功能：

技术分析函数分析投资。可视化金融时间序列的对象；

提供两种创建金融时间序列的对象（用构造器和转换文本文件）；

Fixed－Income Toolbox扩展了Matlab在金融财经方面的应用，可以用固定收益模型进行计算，例如定价、收益和现金流动等有价证券的固定收益计算。支持的固定收益类型包括有价证券抵押回报、社会债卷和保证金等。该工具箱还能够处理相应金融衍生物的计算，支持抵押回收有价证券、国债和可转换债卷等的计算。

Garch Toolbox 提供了一个集成计算环境，允许对单变量金融时序数据的易变性进行建模。 Garch Toolbox使用一个广义ARMAX/GARCH复合模型对带有条件异方差的金融时序数据进行仿真、预测和参数识别。 Garch Toolbox提供了基本工具为单变量广义自回归条件异方差GARCH（Generalized Auto Regressive Conditional Heteroskedasticity）易变性进行建模。

Garch Toolbox采用单变量GARCH模型对金融市场中的变化性进行分析。 4. Fixed－Income Toolbox 5. Garch Toolbox

上述工具箱基本上囊括了通常的金融计算，适用于金融学术研究，特别适合金融实务工作者进行金融计算。 Financial Toolbox提供了一个基于Matlab的财务分析支撑环境，可以完成许多种财务分析统计任务；从简单计算到全面的分布式应用，财务工具箱都能够用来进行证卷定价、资产组合收益分析、偏差分析和优化业务量等工作。 金 融 数 据 统 计

本讲主要介绍了统计学的基本原理和基本统计量。要求读者掌握均匀分布、正态分布随机数生成办法，学会常用的统计绘图命令，掌握回归的方法，学会运用主成份、因子分析金融问题。

一、随机模拟基本原理

1977年，菲利浦.伯耶勒（Phelim Boyle）提出了模拟方法求解金融资产定价问题。其想法是假设资产价格分布是随机波动，如果知道了这个波动过程，就可以模拟不同的路径；每做完一次模拟，就产生一个最终资产价值，在进行若干次这样的过程，那么所得到的结果就是一个最终资产价值分布，从这个分布中可以得到期望的资产价格。 （一） 随机数生成函数

在Matlab中 unidrnd 函数可以生成1～N的均匀分布随机数。其调用方式为： R＝unidrnd( N ) 随机数矩阵

确定输出随机矩阵R的行数 生成在1～N之间的一个随机数 1. 均匀分布随机数生成函数 R＝unidrnd( N,m ) R＝unidrnd( N,m,n )

确定输出随机矩阵R的列数 输出方阵 >>unifrnd(0,1) ans=

0.4565

如果需要生成服从连续均匀分布的随机数，则可以调用 unifrnd 函数，其调用方式为： R＝unifrnd( A,B )

A，B是随机数的下界与上界 如：生成一个0～1之间随机数： 2. 生成服从连续均匀分布的随机数 R＝unifrnd( A,B,m ) R＝unifrnd( A,B,m,n ) m，n表示随机数的维数

下面介绍两种方法生成1～2之间随机矩阵K，K为5行6列矩阵。 方法1 方法2

>>unifrnd(1,2,[5,6]) ans =

1.9334 1.1338 1.5751 1.0129 1.6124 1.5869 1.6833 1.2071 1.4514 1.3840 1.6085 1.0576 1.2126 1.6072 1.0439 1.6831 1.0158 1.3676 1.8392 1.6299 1.0272 1.0928 1.0164 1.6315 1.6288 1.3705 1.3127 1.0353 1.1901 1.7176 >>unifrnd(1,2,5,6) ans =

1.6927 1.1536 1.5548 1.2731 1.9084 1.6408 1.0841 1.6756 1.1210 1.2548 1.2319 1.1909 1.4544 1.6992 1.4508 1.8656 1.2393 1.8439 1.4418 1.7275 1.7159 1.2324 1.0498 1.1739 1.3533 1.4784 1.8928 1.8049 1.0784 1.1708 R＝normrnd( mu,sigma )

正态分布的均值 随机矩阵R的行数 正态分布的方差 3. 生成正态分布的随机数 R＝normrnd( mu,sigma,m ) R＝normrnd( mu,sigma,m,n ) 随机矩阵R的列数 调用方式为： >> normrnd(0,1) ans=

-0.4326

如：生成均值为0，方差为1正态分布的随机数，可用命令

下面用两种方法生成均值为0，方差为1的正态分布矩阵，矩阵为5行6列。 方法1 方法2

>>normrnd(0,1,[5,6]) ans =

-0.3179 0.7310 -0.2556 0.1184 0.7990 -1.0078 1.0950 0.5779 -0.3775 0.3148 0.9409 -0.7420 -1.8740 0.0403 -0.2959 1.4435 -0.9921 1.0823 0.4282 0.6771 -1.4751 -0.3510 0.2120 -0.1315 0.8956 0.5689 -0.2340 0.6232 0.2379 0.3899 >>normrnd(0,1,5,6) ans =

0.0880 0.7812 -2.2023 0.0215 -1.0559 -1.1283 -0.6355 0.5690 0.9863 -1.0039 1.4725 -1.3493 -0.5596 -0.8217 -0.5186 -0.9471 0.0557 -0.2611 0.4437 -0.2656 0.3274 -0.3744 -1.2173 0.9535 -0.9499 -1.1878 0.2341 -1.1859 -0.0412 0.1286 4. 特定分布随机数发生器

在Matlab中有统一格式随机数发生器，函数名称为random，可以生成许多服从不同分布的随机数。

调用方式

y=random('name' , A1 , A2 , A3 , m , n) 输出参数

name 说明随机分布的类型，具体如下表所列。 类 贝 塔 二项分布 卡 方 指数分布 别 分布 简写 类 均匀分布 泊松分布 T 别 布 分布 简写 Uniform unif Poisson poiss 分 正态分布 非中心F分布 非中心T分布 Noncentral F Noncentral T ncf nct Beta beta F 分布 伽 码 对数正态 binomial Chisquare Exponential F bino chi2 exp f Gamma Lognormal gam logn T Normal t norm 特定分布的参数表

>> a = random ('Normal' , 0 , 1 , 3 , 2) a =0.6565 -0.2624 -1.1678 -1.2132 -0.4606 -1.3194

下面用random函数生成3行2列的正态分布随机数矩阵，正态分布的均值为0、方差为1。 5. 多元正态分布的随机数

多元正态分布的随机数可以用如下形式表示: 式中： μ是均值向量， ∑是协方差矩阵。 Xi ～ N(μ , ∑)

mu 均值 sigma 协方差 cases 样本个数

在Matlab中可使用mvnrnd函数生成多元正态分布函数。 调用方式

R= mvnrnd(mu , sigma)

R= mvnrnd(mu , sigma , cases) 输入参数

>>mu = [2 3] ; %均值 >>SIGMA=[1 1.5;1.5 3]; %协方差矩阵 >>r=mvnrnd(mu,SIGMA,100); %生成100个随机样本 >> plot(r(:,1),r(:,2),'+')

下面生成一个多元正态分布的例子。 样本的散点图如右所示： 二元正态分布的散点图

（二）多元正态分布密度函数 >> mu = [1 -1];

>> Sigma = [0.9 0.4 ; 0.4 0.3]; >> X = [ 2 1 ];

>> p = mvnpdf(X , mu , Sigma) p =1.3828 e-005

多元正态分布的密度函数是 mvnpdf。 调用方式

mvnpdf(X , mu , Sigma) 下面是一个例子。

>> mu = [1 -1];

>> Sigma = [0.9 0.4 ; 0 .4 0.3]; >> X = [2 1];

>> f = mvncdf (X , mu , Sigma) f =0.8541

F (x，y) ＝ P(X< x，Y< y)

如果计算分布函数，则X、Y为二元随机正态分布，分布函数 F(x,y) 的定义如下： 调用方式为：

p=mvncdf(X , mu , SIGMA) 下面举一个例子。

可以看出：对于随机变量 X，Y，有

P(X < 2 , Y < 1) = 0.8541 也即 X < 2且Y < 1 的概率为0.8541。 二、随机变量的数字特征

A 如果A为向量，则返回值为该向量的平均值； 如果A是矩阵，则返回值是每列的平均值。 dim dim=1（默认）表示每列平均，dim=2表示每 行平均。 （一）计算平均值 调用方式

M = mean (A)

mean(A , dim) 输入参数

在Matlab中计算几何平均数的函数为 geomean；计算调和平均数的函数是 harmmean 函数，调和平均数的计算公式是

注意样本数据不能为 0 。 >> a = [1 2 ; 3 4]; a =

1 2 3 4 下面是一个例子。 >> mean ( a ) ans =

2 3 >> mean( a , 2) ans =

1.500 0 3.500 0

（二） 剔除异常值后的平均值

X 样本观察矩阵。

percent 剔除比率，例如percent＝10表示同时剔除最大

的5%和最小的5％观察值 。 dim dim＝1（默认）表示对每列求平均值，dim＝2 表示对每行求平均值。 有时观察数据中有异常大或异常小的值，这些异常值会对平均值产生影响，需要去掉异常值。例如在体操比赛中，去掉一个最高分和最低分，然后计算运动员的最后得分。在Matlab中也有剔除异常值后的平均数。 调用方式

M = trimmean ( X , percent )

M = trimmean ( X , percent , dim ) 输入参数

>> x = rand( 1 , 20 ) >> trimmean （x , 10） ans =

0.5145 （三）计算中位数

A 样本观测矩阵

dim dim＝1（默认）表示对每列求中位数，dim＝2表示对每行求中位数 剔除10％的异常值之后的平均数为0.5145 。 调用方式

M ＝ median （A）

M ＝ median （A ，dim） 输入参数

有时数据中出现NaN，在计算中位数时需要忽略NaN，这时需要调用nanmedian函数。 （四）计算方差与标准差 A 样本值

flag 0（默认值）表示方差计算公式为 1表示方差计算公式为

一般说来，资产组合的风险越大，方差越大，波动性越大。方差由于其简单、直观以及良好的统计性质使其成为风险的代名词。

在Matlab中计算方差、标准差的函数分别是 Var、Std。 方差调用方式 Var( A )

Var( A ，flag ) 标准差调用方式 Std（A）

Std（A , flag） 输入参数

（五） 计算样本的百分位数 >> x = rand( 1 , 20 ); >> prctile ( x , 20 ) ans =

0.1688 调用方式

Y ＝ prctile（X，p，dim） 输入参数

X 观察值

p 计算大于p％值的数 dim 同上 下面是一个例子 （六）计算样本极差 r = range ( q )

r = range ( q,dim )

极差就是样本极大值与极小值的差，反映样本的离散程度。 调用方式

>> x = rand ( 1 , 20 ); >> range ( x ) ans =

0.8404

下面是一个例子

（七）计算偏度与峰度

方差作为风险的度量指标并不是完整的。比如讲，两种资产收益分布的均值和方差都是相同的，但是一种资产收益是左偏的，另一种是右偏的。对于风险而言，相对于大概率和小损失人们更加厌恶小概率大损失的情况，后一种情况给人们带来的痛苦远大于第一情况。从这个意义上讲，收益分布左偏的资产的风险水平要小于右偏的资产。此时，方差作为风险的度量指标就不是完全的，还要考虑峰度、偏度等指标。 偏度和峰度是两个高阶的统计量。计算偏度的目的在于考察组合收益率水平是否是对称分布，也就是组合产生亏损与获得盈利的概率如何；峰度是考察组合的收益率情况是否接近正态分布，如果组合的收益率存在尖峰厚尾的分布特征，则说明组合产生亏损和盈利的概率偏大，也就是在一定程度上认为组合收益率出现极端性的可能性偏大，这种组合的收益率稳定性是比较差的。

正态分布的峰度等于 3，大于3表示尖峰，小于3表示分布比较均匀。股票市场收益率的时间序列大都为尖峰肥尾。 偏度的计算公式为

式中：μ， σ分别为样本 x 的均值与方差。

如果 skewness＝0，则表示分布形态与正态分布偏度相同；如果 skewness > 0，则表示正偏差数值较大，长尾巴拖在右边；如果skeqness<0，则表示负偏差数值较大，长尾巴拖在左边。 峰度的计算公式为

>> x = rand( 1 , 20 ); >> skewness( x ) ans =

-0.0487 1. 计算偏度 调用方式

Y = skewness ( A )

Y = skewness( A , flag )

输入参数

A 样本值

flag 偏度的计算方式，0（默认）为无偏估计， 1为有偏估计 下面是一个例子。 k = kurtosis ( X )

k = kurtosis ( X, flag )

k = kurtosis ( X, flag , dim ) 2. 计算峰度 调用方式

X 样本观察矩阵

flag 峰度的计算方式，0（默认）为无偏估计， 1为有偏估计

dim dim＝1（默认）表示对每列求平均， dim＝2表示对每行求平均 输入参数

>> x = rand( 1 , 20 ) >> kurtosis ( x ) ans =

1.4407

下面是一个例子。 （八）计算绝对离差

绝对离差是以偏差绝对数来衡量决策方案的风险。在期望值相同的情况下，绝对离差越大，风险越大；绝对离差越小，风险越小。 调用方式

Y = mad( X )

Y = mad( X , n ) 输入参数

X 观察值

n 绝对偏差计算方式

n＝0（默认）计算公式为mean(abs(X－mean(X))) n＝1计算公式为median(abs(X－median(X))) >> x = rand ( 1 , 20 ) >> mad ( x ) ans =

0.1750

下面是一个绝对离差的例子。 （九）计算中心矩

数理统计中经常用到中心矩的概念，k阶中心矩的计算公式为

可以看出 1 阶中心矩为 0，如果观察值是矩阵则以每列为样本计算中心矩。 X 观察样本值

order 中心矩的阶数，必须为正整数 调用方式

M = moment ( X , order ) 输入参数

>> X = rand ( [ 6 5 ] ) X =

0.4154 0.9901 0.3200 0.4399 0.1338 0.3050 0.7889 0.9601 0.9334 0.2071 0.8744 0.4387 0.7266 0.6833 0.6072 0.0150 0.4983 0.4120 0.2126 0.6299 0.7680 0.2140 0.7446 0.8392 0.3705 0.9708 0.6435 0.2679 0.6288 0.5751 >> m = moment ( X , 3 ) m =

-0.0113 0.0014 0.0032 -0.0058 -0.0023 下面计算样本的 3 阶矩。 （十）计算协方差和相关系数

协方差是一个用于衡量投资组合任意两个资产相关性的统计指标。当协方差为正值时，表示两种资产的收益率呈同方向变动；协方差为负值时，表示两种资产的收益率呈相反方向变化；协方差等于 0 时，表示两种资产不存在相关性。

相关系数总是在 -1～1 之间的范围内变动，-1 表示完全负相关（反向），1 表示完全正相关（同向），0 则表示不相关。Matlab计算协方差、相关系数的函数分别是cov 和 corrcoef 。 协方差 调用方式 C = cov ( X )

C = cov( x , y ) 下面是一个例子

>> A = [ -1 1 2 ; -2 3 1 ; 4 0 3]; >> cov ( A ) ans =

10.333 -4.1667 3.0000 -4.1667 2.3333 -1.5000 3.0000 -1.5000 1.0000 X 观察值矩阵 Y 观察向量

param1 参数 1 ，参数的值如下：

alpha表示置信度，在 0～1 之间 val1 参数 1 的值 param2 参数 2

val2 参数 2 的值 2. 相关系数 调用方式

R = corrcoef ( X )

R = corrcoef ( x , y )

[ R , P ] = corrcoef ( X , 'param1' , val1 , 'param2', val2 , … ) 输入参数

R 相关系数矩阵

P 每个相关系数的概率矩阵 输出参数

>> x = rand ( 30 , 4 ); >> x (: , 4) = sum( x , 2 ); >> [ r , p ] = corrcoef ( x ) >> [ i , j ] = find ( p < 0.05 ); >> [ i , j ]

下面是一个计算相关系数的例子。 r =

1.0000 0.1412 -0.1954 0.4993 0.1412 1.0000 -0.1312 0.5848 -0.1954 -0.1312 1.0000 0.3729 0.4993 0.5848 0.3729 1.0000 p =

1.0000 0.4566 0.3008 0.0050 0.4566 1.0000 0.4896 0.0007 0.3008 0.4896 1.0000 0.0424 0.0050 0.0007 0.0424 1.0000 ans =

4 1 4 2 4 3 1 4 2 4 3 4 三、 统 计 绘 图

>> x = [ 1 2 3 5 7 3 3.4 ] x =

1.0000 2.0000 3.0000 5.0000 7.0000 >> tabulate ( x )

Value Count Percent

1 1 14.29% 2 1 14.29% 3 2 28.57% 3.4 1 14.29% 5 1 14.29% 7 1 14.29% （一）样本频率分布图

样本频率分布图函数是 tabulate 。

下面调用 cdfplot 函数绘出 x 的分布图。 >>cdfplot ( x ) 向量 x 的分布图

（二）最小二乘拟合图

3.0000 3.4000

在 Matlab 中绘制最小二乘拟合图的命令是 lsline ，下面是一个例子。 >>x=rand ( 1 , 20 ) >>x=cumsum ( x ) >>plot ( x , '+' ) >>lsline

最小二乘拟合图

（三）正态分布概率图

有时需要判断样本数据是否服从正态分布，normplot 函数用图的形式给出直观的判断。如果数据点越靠近直线则分布越近似正态分布，越远则越偏离正态分布。 >>x=normrnd(0,1,20,1) >>plot(x,'+') >>normplot(x) 正态分布拟合图 从图中可以看出，数据点基本上是直线，所以符合正态分布。如果判断数据是否服从 Weibull 分布则可以对生成的数据用 wblplot 函数进行判断。下图给出了Weibull 分布拟合的结果。 从图中看出对于数据较小、较大的点偏离较大，数据不服从Weibull 分布。 Weibull分布拟合图 （四）样本密度图 >>randn('seed',0); >>x=randn(1,20); >>x=cumsum(x) >>capaplot(x,[0,10])

在Matlab中绘出样本数据的密度图函数为 capaplot。 样本的密度示意图如右图所示。 样本的密度示意图 （五）频率直方图

Y 观察值。如果 Y 是一个向量，则画出一个频率图; 如果 Y 是一个 m×p 阶矩阵，则对 Y 每一列分别 作频率图

nbins 频率图分成 nbins 等分的区间段，默认值为 10 。 调用方式

n = hist ( Y )

n = hist ( Y , nbins )

[ n , xout ] = hist ( Y , nbins ) 输入参数

n 样本落在区间段的频率 xout 区间断的刻度 输出参数

下面是一个例子。 >>randn('seed',0) >>x=randn(1,200) ; >>hist(x)

频率直方图如右图所示

>> hist ( x , min ( x ) : 0.3 : max ( x ) ) ;

其中 min ( x ) : 0.3 : max ( x ) 表示频率图 X 轴的刻度起点是 x 最小元素，终点是 x 中最大元素，刻度间隔 0.3。规定刻度间隔的频率直方图如下图所示。 如果在频率图的基础上加上正态分布拟合图，则可以用 histfit 函数。 >>randn('seed',0) >>x=randn(1,20) >>histfit(x)

带有密度函数的频率直方图如右图所示。 （六）盒 图

X 样本观察值 G 各组的名称 Param1 参数 1 的名称 val1 参数 1 的值 在 Matlab 中绘制样本数据的盒图函数是 boxplot。 调用方式 boxplot ( X )

boxplot ( X , G )

boxplot ( X , 'Param1',val1,'Param2',val2,…) 输入参数

各参数的名称和内容如下表所列 参 数 名 参 数 值 称 notch symbol orientation whisker labels colors widths positions grouporder ?on?生成有缺口的盒图，?off?生成矩形盒图 线条类型，默认值为?r +? ?vertical?（默认值）为垂直型盒图，?horizontal?为水平型盒图 盒图须线的长度，默认值＝1.5×四分位间距 盒图行列的名称标签 线条颜色 盒图的宽度，默认值＝0.5 盒图的位置，默认值为 1：n 组的次序

>>x1=normrnd(5,1,100,1); >>x2=normrnd(6,1,100,1); >>boxplot([x1,x2],'notch','on')

盒子的上下两条线分别为样本的 25％和 75％分位数，盒子的上下底之间的距离为四分位的间距。

盒子的中间线为样本中值，如果样本中值不在盒子的中间，表示存在一定的偏度。盒子的上方和下方有两条虚线，显示了样本的范围，野值（异常值）位于超过盒子顶端、底端 1.5 倍的四分位数。

含有缺口的盒图中齿形缺口表示样本中值的置信区间。 图 中的内容说明如下： 正态分布盒图如右图所示。 四、多元线性回归分析

在金融上常常需要对金融、经济数据进行回归，其中最简单的是多元线性回归。 （一）多元线性回归 b=regress(Y,X)

[b,bint]=regress(Y,X) [b,bint,r]=regress(Y,X) [b,bint,r,rint]=regress(Y,X)

[b,bint,r,rint,stats]=regress(Y,X)

[b,bint,r,rint,stats]=regress(Y,X,alpha)

假设因变量 Y 和自变量 X 之间服从以下的线性模型：

式中：Y 是因变量的观察值，X是自变量回归矩阵，β 是参数向量，εβ的最小二乘解是 调用方式

X 自变量观察值，注意如果模型中有常数项，则 X 的 第一列所有元素为 1 。 Y 因变量观察值向量 alpha 参数的置信度 输入参数

例1 首先按照下面模型生成一系列随机数，然后回归。 b β 的估计值，注意 b 中已经包含了常数项 bint β 的置信区间 r 残值

rint 残值的置信区间 stats R2、F、概率 p 输出参数

>>[b,bint,r,rint,stats]=regress(Y',[ones(10,1),X'],0.05) 下面生成一组随机数 >>X= 1:10;

>>Y=0.1+0.4*X + normrnd(0,0.1,1,10); 下面估计 β : rint =

-0.1794 0.1149 -0.1435 0.1764 -0.1625 0.1737 -0.2005 0.1417 -0.1795 0.1712 -0.1046 0.2290 -0.0569 0.2470 -0.2559 -0.0402 -0.0667 0.2238 -0.1889 0.1008

是白噪声。 stats =

1.0e+003 *

0.0010 2.2837 0.0000 0.0000 b =

0.1303 0.3953 bint =

0.0120 0.2487 0.3762 0.4144 r =

-0.0323 0.0165 0.0056 -0.0294 -0.0041 0.0622 0.0950 -0.1480 0.0785 -0.0440

从 b 的估计值可以得知常数项和一次项的系数分别为0.1303, 0.3953。在 0.05 置信水平下常数项系数估计区间为[0.0120 0.2487]，X的系数置信区间为[0.3762 0.4144]。由于样本数量非常少，参数估计并不稳定。下图是残差及其置信区间图。 >>rcoplot(r,rint)

（二）多元正态回归

在 Matlab 中 mvnrmle 函数可以进行多元正态回归，假设Yk 为随机变量，其分布如下： 式中：N ( g , g )为多元正态分布。 调用方式

[Parameters,Covariance,Resid,Info]=

mvnrmle(Y,Design,MaxIterations,TolParam,Tol0bj,Covar0) Parameters 参数 Covariance 协方差 Resid 残差

Info 估计过程的相关系数

Y 观察值矩阵，Yn×k中n是样本的个数，k是资产的数目 Design 自变量单元变量矩阵，如果 Y 只有一个资产时，Design是

一个矩阵，如果 Y 中的资产个数大于一个时，Design 是 一个单元向量，每个元素都是一个矩阵。Y 的 k 列和 Design 第 k 个元素中的矩阵进行回归 MaxIterations TolParam Tol0bj Covaro 输出参数

输入参数

（三）估计多元正态分布每个资产的标准差

Data 观察值矩阵，Yn×k 中n样本的个数，k是资产的数目 Design 自变量单元变量矩阵，如果 Y 只有一个资产时，Design

是一个单元向量，如果Y含有多于一个资产时，Design 是一个单元变量矩阵 Covariance 回归时的残值 调用方式

[StdParameters, StdCovariance]=

mvnrstd(Data, Design, Covariance) 输入参数

StdParameters 每个资产的标准差 StdCovariance 协方差 输出参数

（四）岭 回 归

线性回归中参数估计 ，如果观察值 X 存在自相关性，则 XTX是奇异矩阵，估计值就会出现非常大的误差，这时矩阵XTX 需要加上一个对角元素是常数 k 的单位阵，即 。Matlab提供了岭回归 ridge 函数求解该问题。

b 模型估计参数 调用方式

b1 = ridge ( Y , X , k )

b0 = ridge ( Y , X , k，0 ) 输入参数

Y 因变量观察值 X 自变量观察值

k k 表示控制系数，可以根据需 要进行选择。 输出参数 k = 0:0.01:1;

b = ridge(heat, ingredients, k); plot(k, b');

xlabel('Ridge parameter'); ylabel('Standardized coef.'); title('Ridge Trace for Hald Data') legend('x1','x2','x3','x4');

例2 对 hald 文件中的数据进行岭回归。 >>load hald

查看工作区中的变量。 >>who

Your variables are :

hald heat ingredients 五、主成分分析

主成分分析是在各个变量之间相关关系研究的基础上，用较少的新变量代替原来较多的变量，而且使这些较少的新变量尽可能多地保留原来较多的变量所反映的信息。

（一）主成分分析基本原理

首先对样本进行标准化处理，为简单起见，标准化后的样本仍记为X1 , X2 , X3 , ? , Xp 。 设样本为 X1 , X2 , X3 ,? , Xp，其对应的样本均值为

对应的标准差为 S1 ，S2，S3 , ? , Sp 。

设 F1 , F2 , F3 , ? , Fp 是主成分，也即是 X1 , X2 , X3 , ? , Xp 的线性表示，同时满足下面的条件：

主成分是原样本的正交变换。 各主成分之间互不相关。 主成分的总方差不变。

主成分按方差从大到小排序。 主成分具有如下性质：

这一性质说明，主成分是原变量的线性组合，是对原变量信息的一种改良；主成分不增加总信息量，也不减少总信息量。

设 为主成分的特征值，前 k 个方差累积贡献率为 一般当累积贡献率大于 85％ 时不再增加新的主成分。

保留多少个主成分取决于保留部分的累积方差在总方差中占的百分比（即累计贡献率），它标志着前几个主成分概括的信息的多寡。在实践中，粗略规定一个百分比就可以决定保留几个主成分，如果多留一个主成分，但累积方差增加无几，便不再多留。 （二）主成分分析函数

COEFF 主成分系数 SCORE 新坐标系

latent X的协方差矩阵的特征值 tsquare Hotelling 统计量的值

在Matlab中提供了两个主成分分析函数princomp和pcacov。 [COEFF,SCORE]=princomp(X)

[COEFF,SCORE,latent]=princomp(X)

[COEFF,SCORE,latent,tsquare]=princomp(X) 输入参数

X 观察变量 输出参数 调用方式

>>corrcoef ( ingredients ) ans =

1.0000 0.2286 -0.8241 -0.2454 0.2286 1.0000 -0.1392 -0.9730 -0.8241 -0.1392 1.0000 0.0295 -0.2454 -0.9730 0.0295 1.0000

例3 用Matlab自带数据进行主成分分析。Matlab中自带了数据文件 hald ，可以直接调用进行主成分分析。hald 数据考虑影响温度的 4 个因素，温度保存在 heat 变量中，4 个因素的观察值保存在 ingredients 中。由于 4 个因素之间存在相关性，无法直接回归，为解决这个问题，首先进行主成分分析，生成四个主成分变量，主成分之间互不相关，而且和观察值的信息是同样的。

第一步：载入数据，考察变量之间的相关性。

>> load hald %载入Matlab自带的数据文件

考察相关性

发现自变量 2 与变量 3 之间的高度相关，所以需要剔除相关性。 第二步：主成分分析。

>> [pc,score,latent,tsquare] = princomp(ingredients) 主成分系数 pc =

0.0678 0.6460 -0.5673 0.5062 0.6785 0.0200 0.5440 0.4933 -0.0290 -0.7553 -0.4036 0.5156 -0.7309 0.1085 主成分的方差贡献率 score =

-36.8218 6.8709 -29.6073 -4.6109 12.9818 4.2049 -23.7147 6.6341 0.5532 4.4617 10.8125 3.6466 32.5882 -8.9798 -22.6064 -10.7259 9.2626 -8.9854 3.2840 14.1573 -9.2200 -12.3861 25.5849 2.7817 26.9032 2.9310 协方差的特征值 latent = 517.7969 67.4964 12.4054 0.2372 tsquare = 5.6803 3.0758 6.0002 2.6198 3.3681 0.5668 3.4818 3.9794 2.6086 7.4818 4.1830 2.2327

0.4684 4.5909 2.2476 -0.9022 -1.8547 6.0874 -0.9130 1.6063 -3.2365 0.0169 -7.0465 -3.4283 0.3867 2.4455 0.4844 0.3967 -0.3958 -1.1261 -0.3786 0.1424 -0.1350 0.0818 0.3243 -0.5437 0.3405 0.4352 0.4468 0.4116 2.7216

由此可以得到4个主成分如下：

从特征值可以看出前面两个主成分可以很好的解释98%的方差： (517.796 9+67.496 4)/(517.7969+67.4964+12.4054+0.2372) =98%

>>covx=cov(ingredients);

>>[pc,latent,explain]=pcacov(covx)

采用 pcacov 函数计算主成分的结果同 princomp 函数的结果是一样的。 >> load hald

然后计算观察变量 ingredients 协方差。 >> new = ingredients * pc

在确定主成分后对主成分进行回归。

将生成的 4 个主成分，保存在变量 new 中，变量 new 中的每列就是一个主成分。 >> regress(heat,new(:,1:2)) ans =

2.1843 1.0894

验证主成分之间的相关性： >> corrcoef(new) ans =

1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 发现主成分之间是不相关的。

第三步：用前两个主成分进行回归。 >> 2.1843 * pc(:,1) + 1.0894 * pc(:,2) ans =

0.8519 1.5039 -0.8862 -1.4783

这样温度和主成分之间的关系如下： 还原成线性回归形式：

这样温度和自变量之间的关系如下： 第四步：验证主成分分析优点。

下面计算 heat 和 ingredients 中各个分量的相关系数。 >> corr(heat, ingredients(:,1)) ans =

0.7307

>> corr(heat, ingredients(:,2)) ans =

0.8163

>> corr(heat, ingredients(:,3))

ans =

-0.5347

>> corr(heat, ingredients(:,4)) ans =

0.8213

相关系数表明前两个自变量和因变量之间是正相关，后两个自变量和因变量之间负相关。式(1)和实际结果相印证。

从式(2)可以看出自变量和因变量之间都是正相关，显然和实际结果不相印证。 自变量和因变量之间的函数关系如下： >> c0 = regress(heat,ingredients) c0 =

2.1930 1.1533 0.7585 0.4863

如果不用主成分回归，而是直接对自变量和因变量进行回归。 六、因子分析

因子分析法（factor analysis）是一种用来分析隐藏在表象背后的因子作用的一类统计模型和方法，它起源于心理度量学，最初是研究如何使用少数几个变量来解释众多原始变量，同时又尽量避免信息丢失的多元统计分析方法。在实际问题的分析过程中，常用因子分析去除重叠信息，将原始的众多指标综合成较少的几个因子变量来分析。 1. 因子分析法相关概念

设样本为 , 其对应的样本均值为

, 对应的标准差为 。 因子分析的核心问题有两个，一是确定因子的个数，二是对因子变量进行命名。

首先对样本进行标准化处理，为简单起见，标准化后的样本仍记为 。 因子分析的数学模型如下： X = f X +s

式中：X是经过正交处理的样本值，f 是负荷矩阵，表示公共因子部分，s是其他特殊因子。要求f，s之间不相关而且不可观察。 lambda 负荷矩阵估计值

psi 特殊因素矩阵估计值 T 因素负荷旋转矩阵 stats 假设检验

F F是n×m维因子的分矩阵 2. 因子分析函数调用方式 调用方式

lambda = factoran(X,m)

[ lanmba , psi ] = factoran(X,m) [ lanmba , psi , T ] = factoran(X,m)

[ lanmba , psi , T , stats ] = factoran(X,m) [ lanmba , psi , T , stats , F] = factoran(X,m) 输入参数

X 观察值矩阵，每列属于同一个变量 m 公共因素的个数 输出参数

下面用一个例子说明

例4 在Matlab中自带了10个公司100周的收益率数据，10家公司分成3个行业，4家属于科技公司，3家属于金融公司，另外3家属于零售类公司。一般来讲，同一行业内上市公司同质性强，走势也相同，不同行业内之间股票受不同风险因素影响，走势差别较大。下面验证该结论。 >> load stockreturns

>> m=3; ％考虑3个行业的因素

>>[LoadingsPM,specificVarPM] =factoran(stocks,m,'rotate','promax') 从上述结果可以明显看出第1到第4家公司属于同一类公司，第5家到第7家属于同一类公司，第8家到第10家属于同一类公司，说明不同行业间的股票收益率存在不同。 LoadingsPM =

0.9452 0.1214 -0.0617 0.7064 -0.0178 0.2058 0.3885 -0.0994 0.0975 0.4162 -0.0148 -0.1298 0.1021 0.9019 0.0768 0.0873 0.7709 -0.0821 -0.1616 0.5320 -0.0888 0.2169 0.2844 0.6635 0.0016 -0.1881 0.7849 -0.2289 0.0636 0.6475 specificVarPM = 0.0991 0.3431 0.8097 0.8559 0.1429 0.3691 0.6928 0.3162 0.3311 0.6544 七、方差分析

方差分析主要用来检验两个以上样本的平均值差异的显著程度，由此判断样本究竟是否抽自具有同一均值的总体。 （一）单因素方差分析

一般地，假定所检验的结果受同一因素 A 的影响，可以取 k个不同的水平 1, 2, 3, ?, k。对于因素的每一个水平 i 都进行 n 次试验，结果分别为 , 一般把这一组样本记作Xi。假定 ，即对于因素的每一个水平，所得到的结果都服从正态分布，且方差相等。 在方差分析中，通常把对实验结果发生影响和起作用的自变量称为因素。如果分析一个因素

对于试验结果的影响和作用，就称为单因素方差分析。 （二）方差分析步骤

方差分析一方面确定因素的不同水平下均值之间的方差，作为所有试验数据所组成的全部总体方差的一个估计值；另一方面，再考虑同一水平下不同试验数据对于这一水平的均值的方差。由此，计算出全部数据的总体方差的第二个估计值；最后，比较上述两个估计值。如果这两个方差的估计值比较接近就说明因素的差异并不大，则接受零假设；否则接受备择假设。根据上述思路可以得到方差分析的方法和步骤。 1. 提出假设 H0:μ1=μ2=?=μk，即因素的不同水平对试验结果无显著影响。

H1:不是所有的μi都相等(i=1,2,?,k)，即因素的不同水平对试验结果有显著影响。 2. 方差分解

下面先定义总离差平方和为各样本观察值与总均值的离差平方和。记 将总离差平方和分解为两部分： N = nk为样本观察值总数 记

表示同一样本组内，由于随机因素影响所产生的离差平方和，简称组内平方和。 记

表示不同的样本组之间，由于变异因素的不同水平影响所产生的离差平方和，简称为组间平方和。

对应于 SST，SSR 和 SSE 的自由度分别为 n-1，k-1，n-k。 由此可以得到： SST = SSR + SSE 当原假设 H0：μ1=μ2=?=μk 成立时，E(MSE)=E(MSR)

=σ2。此时MSR较小，F值也较小。反之H0不成立时，MSR较大，F值也较大。对于给定的显著性水平α查F分布表得到 F1- α(k-1,n-k)。 如果F > F1- α(k-1,n-k) ，则原假设不成立，即 k 个组的总体均值之间有显著的差异，就拒绝H0。若F ≤ F1- α(k-1,n-k) ，则原假设成立，即k个组的总体均值之间没有显著的差异，就接受H0。 3. F检验

将SSR和SSE分别除以自由度，即得各自的均方差： 组间均方差 MSR = SSR/(k-1) 组内均方差 MSE = SSE/(n-k) 检验统计量

由于方差分析表（见下表）概括了方差分析的统计量之间的关系，在进行方差分析时就可以直接按照方差分析表来逐行、逐列地计算出有关的统计量，最后得到检验量 F的值，并把这一 F 值与查表所得到的一定显著性水平下的 F检验的临界值进行比较，然后作出接受或拒绝原假设的结论。 4. 方差分析表

上述方差分析的方法可以用一张标准形式的表格来实现，称为方差分析表。方差分析表分为五列：第一列表示方差的来源；第二列表示方差的离差的平方和；第三列表示自由度；第四列表示均方差；第五列表示统计检验量 F。表格又分为三行：第一行是组间的方差SSR和均方差 MSR，表示因素的不同水平的影响所产生的方差，其值作为计算统计检验量 F 时的分子；第二行是组内方差 SSE 和均方差 MSE，表示随机误差所引起的方差，其值作为计算统计检验量 F 的分母；第三行是检验行，表示总的方差 SST。 单因素方差分析表

方差来离差平自由度 均方统计检源 组间差 组内差 总方差 方和 SSR SSE SST k-1 n-k n-1 差 MSR MSE 验量 F F = MSRMSE / （三）单因素方差分析函数 p = anova1(X)

p = anova1(Y,group)

p = anova1( Y,group,'displayopt' )

[p,table] = anova1( Y,group,'displayopt' )

[p,table,stats] = anova1( Y,group,'displayopt' )

单因素方差分析是比较两组和多组样本的均值，假设各组变量之间相互独立，方差相等，且服从正态分布。原假设是各组均值全部相等。 调用方式 输入参数

X 样本观察值，要求各列均为彼此独立的样本 Y 观察值向量

group 组别，Y中每个元素所属的类别 displayopt 取值为“off”与“on”，分别表示掩藏与显示方

差分析表图和盒图。盒图上下线为25%和75% 分位数，中间线为中位数。 输出参数

p 各列均值相等的概率 table 方差分析表 stats 统计量 ST达声 恒瑞医药 日 期 价 格 收益率 日 期 价 格 收益率 2006 – 8 2.48 - 14 7 2006 – 8 - 13.191 14 2006 – 8 2.41 - 0.030 96 2006 – 8 - 13.267 - 15 15 0.005 762 2006 – 8 – 16 2.5 0.373 44 2006 – 8 13.066 - 0.015 – 16 15 2006 – 8 2.48 – 17 - 0.008 2006 – 8 13.167 0.007 73 – 17 0.002 202 0 2006 – 8 2.43 - 0.020 16 2006 – 8 13.196 – 18 – 18 2006 – 8 2.41 -0.008 23 2006 – 8 13.196 – 21 2006 – 8 – 22 2006 – 8 – 23 2.4 0 – 21 2.4 - 0.004 15 2006 – 8 13.309 – 22 0.008 563 2006 – 8 14.185 0.065 82 – 23 2006 – 8 2.37 - 0.012 5 2006– 8 – 14.143 - 0.002 – 24 24 96 0.000 424 0.009 612 0 2006 – 8 2.42 0.021 097 2006 – 8 14.149 – 25 – 25 2006 – 8 2.48 0.024 793 2006 – 8 14.285 – 28 – 28 2006 – 8 2.49 0.004 032 2006 – 8 14.285 – 29 – 29

例5 恒瑞医药（600276）和ST达声（000007）的股价如下表所列。 恒瑞医药和ST达声价格收益表

>>group={'stds','stds','stds','stds','stds','stds','stds','stds','stds','stds','stds','hryy','hryy','hryy','hryy','hryy','hryy','hryy','hryy','hryy','hryy','hryy'} 下面检验二者的收益率是否相等

>> rate = [ - 0.03096 0.037344 -0.008 - 0.02016 - 0.00823 - 0.00415...

0 - 0.0125 0.021097 0.024793 0.004032 0.005762 ...

- 0.01515 0.00773 0.002202 0 0.008563 0.06582 ...

- 0.00296 0.000424 0.009612 0 ] >> [p,table,stat] = anova1(rate,group,'on') p =

0.2412 table =

'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'

'Groups' [8.7401e-004] [ 1] [8.7401e-004] [1.4976] [0.2412]

'Error' [ 0.0082] [14] [5.8360e-004] [] []

'Total' [ 0.0090] [15] [] [] [] stat =

gnames: {2x1 cell} n: [11 5] source: 'anova1'

means: [3.4636e-004 0.0163] df: 14 s: 0.0242

从结果可知p=0. 2412>0.05，所以不能拒绝原假设，可以认为二者在此期间的收益率相等。其原因可能是以日收益率为样本，间隔时间太短，不足以反映股票之间的差别。 ST达声和恒瑞医药方差分析表如右上图所示 从图中可以看出方差分析表。 第一列显示数据源(source)。

第二列显示各类数据的平方和(SS)。 第三列显示各类数据相应的自由度(df)。 第四列显示均方差MS。

第五列显示Friendman检验的F统计量(F)。 第六列显示F统计量对应的概率值。

ST达声和恒瑞医药收益率盒图如右下图所示 （四）双因素方差分析

设因素 A 有a 个水平，因素 B 有 b 个水平，试验的重复次数记作 n。记 Xijk为在因素 A 的第 i 个水平、因素 B 的第 j 个水平下进行第 k 次试验时的观察值（i = 1,2,?,a; j =

1,2,?,b; k = 1,2,?,n）。记

前面所研究的是试验结果仅受一个因素影响的情形。要求检验的是当因素取两个不同水平时对结果所产生的影响是否显著。但在实践中，某种试验结果往往受到两个或两个以上因素的影响。双因素方差分析的基本思想与单因素方差分析基本相同。首先分别计算出总变差、各个因素的变差以及随机误差的变差；其次根据各变差相应的自由度求出均方差；最后计算出 F 值并作 F 检验。

为在因素 A 的第 i 个水平、因素 B 的第 j 个水平下进行各次重复试验的所有观察值的总和。

（i = 1,2,?,a; j = 1,2,?,b） 为在因素A的第i个水平、因素B的第j个水平下进行各次重复试验的所有观察值的平均值。 式中：N＝abn 是所有观测值的总数， 是所有观察值的平均值。 对于因素 A：H0：因素 A 的各个水平的影响无显著差异。

H1：因素 A 的各个水平的影响有显著差异。 对于因素 B：H0：因素 B 的各个水平的影响无显著差异。

H1：因素B 的各个水平的影响有显著差异。 对于因素 AB 的交互作用：

H0：因素 AB 的各个水平的影响无显著差异。 H1：因素 AB 的各个水平的影响有显著差异。 利用上面所引入的符号，可以得到有交互作用的两因素方差分析的步骤如下： 建立假设

由于两因素有交互影响，因此除了分别检验两因素单独对试验结果的影响外，还必须检验两因素交互作用的影响是否显著。 2. 离差平方和的分解

总离差平方和 SST 的自由度为N – 1 。

有交互作用的两因素方差分析的总离差平方和可以分解为 4 项： 表示因素间交互作用的离差平方和，自由度为

(N-1)-(a-1)-(b-1)-(n-1)ab=(a-1)(b-1) 分别记

为因素 A 的离差平方和，自由度为 a – 1 。 为因素 B 的离差平方和，自由度为 b – 1 。 表示随机误差的离差平方和，自由度为

N-ab =abn-ab=ab(n-1)。 方差来源 离差平方和 自由度 因素A 因素B SSA SSB a－1 b－1 均方差 统计检验量F MSA＝SSA/（a－1） MSB＝SSB/（b－1） 交互作用 SSAB 误差E 总方差 SSE SST (a－1)(b－1) MSAB＝SSAB/ (a－1)(b－1) n－a b n－1 MSE＝SSE/（N－ab）

有交互影响的双因素方差分析表 编制方差分析表，进行 F 检验

本文来源：https://www.bwwdw.com/article/qm4w.html