matlab在各个学科中的应用

MATLAB在各学科中的运用MATLAB是由美国MATHWORKS公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、FORTRAN）的编辑模式，代表了当今国际科学计算软件的先进水平。1MATLAB和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语言编程的方法来对电路进行仿真分析和计算。结论表明，应用这两种方法可以使复杂电路的分析和计算变得非常快捷、方便，从而为电路分析提供了一个有效的辅助工具。图21是由电压源和电阻组成的简单电路，运用回路电流法，网孔法，节点法等一些经典的电路分析法即可解决此问题。电路也可用SIMULINK进行仿真实验，并通过数据测量等从而检测计算及编程的结果。21基本电路数学模型分析基本电路分析的基本方法是先建立数学模型，一般是电路方程组。然后通过求解方程组，得到各支路电压和电流。图21基本电路图如图21所示设三个回路的电流分别为IM1，IM2，IM3。对图21应用网孔电流法，可列出如下方程网孔1V1R1IM1R3（IM1IM2）0网孔2V2R2IM2R4IM2R3（IM2IM1）0经化简可得（R1R3）IM1R3IM2V1V2R3IM1（R2R3R4）IM2V2整理以上方程，并写成形如AIB的矩阵方程形式，可得R1R3R3IM1V1V2R3R2R3R4IM2V2依题，要求I1、I2、I3，根据节点电流的关系，有I1IM1；I2IM2；I3IM1IM2由以上网孔电流分析法的分析过程可得以下程序流程图设参数，并赋予参数值（与题中所给的要求一致）为系数矩阵中个元素赋值，并得出系数矩阵，为下面计算作准备由AIB得出I矩阵的计算公式，并列出I1、I2、I3的表达式开始结束图22程序流程图22SIMULINK仿真分析用SIMULINK仿真以上电路，并观察测量数据。I1、I2、I3的测量结果再和程序的结果相比较，若相同，则说明电路仿真和程序编写多是正确的；若不同，则说明电路的仿真与程序编写至少一个有问题，认真检查程序及电路的仿真，找出错误，认真分析，纠正错误，在比较结果。使用SIMULINK进行仿真一般分为两步用户首先需要在仿真模型编辑窗口中搭建好自己的模型，设置好具体模型参数和仿真参数；然后用户就可以开始仿真，SIMULINK将根据用户搭建的模型，模型系统在用户设定条件下的具体行为。对于建模，SIMULINK提供了一个图形化的用户的界面（GUI），用户可以像用铅笔在纸上画图一样画模型图。SIMULINK的所有模型是分级的，因此可以通过自上而下或者自下而上来建立模型。编程级仿真31程序编写及模块分析程序如下CLEARR15R26R310R44V115V210给定初始值AR1R3R3R3R2R3R4给出系数矩阵ABV1V2V2给出系数矩阵BIAB求解未知变量矩阵IIM1I1网孔电流的定义IM2I2I1IM1求解电流I1I2IM2求解电流I2I3IM1IM2求解电流I3运行结果如图31所示图31M文件建模仿真结果程序分析程序由21节的设计思路分析所得开始先给元件（与电路图相对应）赋值，因为MATLAB编程计算矩阵的，所以此赋值是必须的，其为下面的的系数矩阵赋值作准备。由以上的元件赋值可得形如AIB矩阵方程形式的系数矩阵，这样使得A，B，I矩阵变为已知。MATLAB提供了两种除法运算左除（）和右除（），一般情况下，IAB是方程IXB的解，而IB/A是方程IAB的解。在传统的MATLAB算法中，右除是先计算矩阵的逆再相乘，而左除则不需要计算矩阵直接进行除运算。通常右除快一点，但左除可避免被除矩阵的奇异性所带来的麻烦。由上所述，由公式IAB可求得I。以上所得I矩阵的元素有IM1、IM2。这四个元素在矩阵I中的顺序也如此。可见I1）IM1，I（2）IM3，由第二节的分析知I1IM1；I2IM2；I3IM1IM2。程序运行后，得到I1、I2、I3的值，并显示出来。与题中要求所求的一样，程序编写正确，并正常运行，与预期一样。32用SIMULINK仿真电路总电路如图32所示图32SIMULINK仿真电路及显示图图33所示为SIMULINK仿真电路图的电流表A3值在示波器中的演示图33SIMULINK仿真电路图的电流表A3值仿真电路图说明图32与实际的电路图对应，下面对以上进行简单说明A1为测量电流I1的器件，其输出电流I1，通过DISPLAY1显示；A2为测量电流I2的器件，其输出电流I2，通过DISPLAY2显示；A3为测量电流I3的器件，其输出电流I3，通过SCOPE显示。电路其余部分与实际电路差不多，值设置也一样。仿真结果分析由以上电路仿真可见，电路的编程求解及仿真所得结果一样。都达到了预期的效果，可见编程及电路仿真没有问题。MATLAB与SIMULINK的结合运用是电路求解问题得以简化，使电路求解问题省时且精炼，精简。应用二MATLAB在自动控制理论中的运用系统频率特性的测量与分析模拟电路图及系统结构图分别如图4L和图42。图41图42取R3500K，则系统传递函数为若输入信号U1TU1SINT，则在稳态时，其输出信号为U2TU2SINT。改变输入信号角频率值，便可测得二组U2U1和随变化的数值，这个变化规律就是系统的幅频特性和相频特性。实验2考虑二阶系统传递函数模型试用MATLAB绘制出不同和N的伯德图。实验3已知系统的开环传递函数为1绘制根轨迹，并根据根轨迹图求若要使系统稳定，K的最大值；2取K10，绘制BODE图，根据BODE图求系统的幅值裕量和相角裕量；取K10，绘制NYQUIST图，并根据NYQUIST图判断系统的稳定性；检验用上述几种方法判断系统稳定性的一致性。三、步骤1连接被测量典型环节的模拟电路。电路的输入U1接AD、DA卡的DAL输出，电路的输出U2接AD、DA卡的ADI输入。检查无误后接通电源。2启动计算机，在桌面双击图标自动控制实验系统运行软件。3测试计算机与实验箱的通信是否正常，通信正常继续。如通信不正常查找原因使通信正常后才可以继续进行实验。测频率图4选中实验课题系统频率特性测量测频率图菜单项，鼠标单击将弹出参数设置窗口。参数设置完成后点确认等待观察波形。测波特图5在测量波特图的过程中首先应选择实验课题系统频率特性测量测波特图数据采样采集信息。6待数据采样结束后点击实验课题系统频率特性测量测波特图图形观测即可以在显示区内显示出所测量的波特图。测奈氏图7在测量波特图的过程中首先应选择实验课题系统频率特性测量测奈氏图数据采样采集信息。8待数据采样结束后点击实验课题系统频率特性测量测奈氏图图形观测即可以在显示区内显示出所测量的波特图。9实验2、实验3涉及的主要命令有BODE（）、RLOCUS、RLOCFIND（）、MARGIN、NYQUIST。为便于比较，可用HOLDON指令将多条曲线放在一个图中。进一步，为清楚起见，用LEGEND指令在图中加注释。用RLOCUS画出根轨迹后，需要时用RLOCFIND（）找出临界稳定点。实验一1、根据模拟电路图，可求出U2U1和随变化的数值，见下表由此可以画出幅频特性和相频特性，如下图所示实验二取6，分别为01、03、07、10、13，绘制其波特图，可得下图取07，分别为2、4、6、8、10，绘制其波特图，可得下图实验三1绘制根轨迹由KRLOCFIND可求出K的最大值；2绘制BODE图和NYQUIST图由图可以得到系统的幅值裕度H14DB,相角裕度458度。相应的截止频率WC225，穿越频率WX114有奈氏判据可得，系统闭环不稳定。应用三基于MATLAB的通信系统仿真系统综述利用MATLAB仿真软件，完成如图1所示的一个基本的数字通信系统。信号源产生0、1等概分布的随机信号，映射到16QAM的星座图上，同时一路信号已经被分成了实部和虚部，后边的处理建立在这两路信号的基础上。实部、虚部信号分别经过平方根升余弦滤波器，再加入高斯白噪声，然后通过匹配滤波器（平方根升余弦滤波器）。最后经过采样，判决，得到0、1信号，同原信号进行比较，给出16QAM数字系统的误码。结构框图图1待构建系统的框图系统实现随机信号的生成利用MATLAB中自带的函数RANDSRC来产生0、1等概分布的随机信号。源代码如下所示定义待仿真序列的维数NGLOBALNN320定义产生1的概率为PGLOBALPP05首先产生随机二进制序列SOURCERANDSRC1,N,1,0P,1P0、1等概分布的随机信号如图2所示。051015202530354045500500511501两两两两两两两两两两两图20、1等概分布的随机信号波形图星座图映射将等概分布的0、1信号映射到16QAM星座图上。每四个BIT构成一个码子，具体实现的方法是，将输入的信号进行串并转换分成两路，分别叫做I路和Q路。再把每一路的信号分别按照两位格雷码的规则进行映射，这样实际上最终得到了四位格雷码。为了清楚说明，参看表1表1两位格雷码的映射规律两位0、1码映射后（按格雷码）003011111103源代码如下所示FUNCTIONY1,Y2QAM_MODULATIONXQAM_MODULATION对产生的二进制序列进行QAM调制首先进行串并转换，将原二进制序列转换成两路信号NLENGTHXA12NY1XAY2XA1分别对两路信号进行QPSK调制对两路信号分别进行24电平变换A12N/2TEMP1Y1ATEMP2Y1A1Y11TEMP12TEMP2TEMP1Y2ATEMP2Y2A1Y22TEMP12TEMP2对两路信号分别进行相位调制A1N/4Y1Y112141COS2PIAY2Y222141COS2PIA按照格雷码的规则进行映射Y1FINDY1103Y1FINDY1111Y1FINDY1131Y1FINDY1123Y2FINDY2203Y2FINDY2211Y2FINDY2231Y2FINDY2223得到的星座图如图3所示，图上注明了每一个点对应的01序列。543210123455432101234516QAM两两两0000001010101000100110110011000101010100011011101111011111101101图316QAM星座图从上边的星座图上可以清楚的看到，任意相邻的两个点之间它们对应的4个BIT中只有一个有差别，也就是格雷码的特点。而采用格雷码主要目的是当信噪比较大时，也就是系统的误码率比较低的情况下，当出现一个符号错误的情况下，往往只是这个符号中的一个BIT位出现了误码，因此这个情况下误码率和误BIT率是41，这一特性在后边的误码率计算的过程中会有应用。插值为了能够模拟高斯白噪声的宽频谱特性，以及为了能够显示波形生成器（平方根升余弦滤波器）的效果，所以在原始信号中间添加一些0点。具体实现是分别在信号的I路和Q路中，任意相邻的两个码字之间添加7个0。源代码如下所示FUNCTIONYINSERT_VALUEX,RATIOX是待插值的序列，RATIO是插值的比例。两路信号进行插值首先产生一个长度等于RATIO倍原信号长度的零向量YZEROS1,RATIOLENGTHX再把原信号放在对应的位置A1RATIOLENGTHYYAX对I路和Q路信号进行插值后的波形图如图4所示。0102030405060708042024两两两两0102030405060708042024两两两两两两两两两两两两图4经过插值后的两路信号波形图波形成形（平方根升余弦滤波器）为了避免相邻传输信号之间的串扰，多元符号需要有合适的信号波形。图1中的方波是在本地数字信号处理时常见的波形，但在实际传输时这种方波并不合适。根据奈奎斯特第一准则，在实际通信系统中一般均使接收波形为升余弦滚降信号。这一过程由发送端的基带成形滤波器和接收端的匹配滤波器两个环节共同实现，因此每个环节均为平方根升余弦滚降滤波，两个环节合成就实现了一个升余弦滚降滤波。实现平方根升余弦滚降信号的过程称为“波形成形”，通过采用合适的滤波器对多元码流进行滤波实现，由于生成的是基带信号，因此这一过程又称“基带成形滤波”。平方根升余弦滤波器的冲激响应基带平方根升余弦滤波器具有以下定义的理论函数0|2SIN121FHFFFN1|FFF当当当其中是奈奎斯特平率，是滚降系数。21SNRTF下面给出平方根升余弦滤波器的冲激响应曲线，如图5所示。051015202530354045010010203040506NSAMPLESAMPLITUDE两两两两两两两两两两两两两两图5平方根升余弦滤波器的冲激响应曲线从上图上不难看出来，平方根升余弦滤波器的冲激响应很显然的引入了符号间干扰（ISI）即它的冲激响应在相邻的抽样点上的值并不象升余弦滤波器那样恒为0。然而造成这一后果的原因在于，当我们引入平方根升余弦滤波器的时候，就是认为整个信道，也就是说，包括信号发送端的滤波器和信号接收端的滤波器，总体的效果是避免了符号间干扰（ISI），所以，单独看这每一个滤波器，勿庸置疑，它们都是存在着符号间干扰（ISI）的。经过平方根升余弦滤波器后的信号I路和Q路信号经过平方根升余弦滤波器后，成形后的波形如图6所示。源代码如下X1、X2是两路输入信号，FD是信号信息位的频率，FS是信号的采样频率FUNCTIONY1,Y2RISE_COSX1,X2,FD,FS生成平方根升余弦滤波器YF,TFRCOSINEFD,FS,FIR/SQRT对两路信号进行滤波Y1,TO1RCOSFLTX1,FD,FS,FILTER/FS,YFY2,TO2RCOSFLTX2,FD,FS,FILTER/FS,YF10倍载波调制将通过成形滤波器后的信号调制到10倍于原频率的载波上。由于在仿真的过程中，只0102030405060708042024两两两两0102030405060708042024两两两两两两两两两两两两两两两两两两两两两两两两图6通过平方根升余弦滤波器后的两路信号能用离散的点来模拟连续信号，因而为了能够显示出一个正弦曲线，至少需要在一个正弦周期内采样到4个以上的点，这里，我们在一个周期内采10个点。假设最初的0、1信号的频率是1HZ，那么I路和Q路符号传输的频率是1/4HZ，而10倍频是建立在I路或Q路符号频率的基础上，也就是说，载频的频率是25HZ。按照前面的假设，那么相邻两个采样点之间的时间间隔是004S。而一个完整周期内的正弦波形的幅值是相同的，都是对应的这个周期内的I路和Q路线性叠加，调制后的信号为，其中，为载波频率。TFTFTYCC2SIN2OSCF源代码如下载波调制X1,X2代表两路输入信号，F是输入信号的频率，HF是载波的频率FUNCTIONT,YMODULATE_TO_HIGHX1,X2,F,HF产生两个中间变量，用来存储插值后的输入信号YO1ZEROS1,LENGTHX1HF/F10YO2ZEROS1,LENGTHX2HF/F10N1LENGTHYO1对输入信号分别进行插值，相邻的两个点之间加入9个点，且这9个点的值同第0个点的值相同YO1NX1FLOORN1/HF/F101YO2NX1FLOORN1/HF/F101生成输出输出信号的时间向量T1LENGTHYO1/HFF/10生成载波调制信号YYO1COS2PIHFTYO2SIN2PIHFT得到的调制到载频的信号波形如图7和图8所示，其中图7主要为了显示一个脉宽周期内的调制信号波形，图8则是显示载波信号的整体情况。10111213141516171819202151050051152两两两两两两两两两两两两图7载波调制信号展开图010020030040050060070025215105005115225两两两两两两两两两两两两图8载波调制信号整体图加入高斯白噪声将通过成形滤波器后的信号送到具有高斯白噪声特征的加性信道中，相当于在原信号上加入高斯白噪声。由于高斯白噪声加在了通过插值和滤波后的点上，因此在计算信噪比的时候存在一个信噪比换算的问题。当我们把仿真得到的误码率曲线同理论的误码率曲线相比较的时候，两者的信噪比的定义必须是一致的。一致包括两个方面，一是二者均为每BIT符号上的信号功率和噪声功率的比值，另一个是信号的功率是指那些信息点上的平均功率，噪声也是指信息点上所对应的噪声的平均功率，但由于噪声的功率谱密度是一个定值，所以噪声的平均功率实际上就是噪声的功率谱密度。对于第二点，由于所有信号的平均功率和信息点上的信号的平均功率不同，所以需要在加入高斯噪声的时候进行纠正，具体的公式推导如下。设是最后理论计算中的信噪比，是加入高斯白噪声后的整体信号（包括插值后SNRSNR的点）的信噪比,是每BIT信息点的平均能量，是每BIT信号的平均能量，是噪EBEBNO声的平均功率，现在需要推导出与的关系。SNRSNROEBEB即两个信噪比的比值就是平均能量的比值。源程序如下对输入的两路信号加高斯白噪声，返回处理后的两路信号，信息点等效BIT信噪比为SNR的值FUNCTIONY1,Y2GENERATE_NOISEX1,X2,SNRSNR1代表SNR对应的符号信噪比SNR1SNR10LOG104算出所有信号的平均功率SSVARX1IX2,1加入高斯白噪声YAWGNX1JX2,SNR110LOG10SS/10,MEASUREDY1REALYY2IMAGY给出加入高斯白噪声的两路信号波形。0102030405060708042024两两两两两两0102030405060708042024两两两两两两两两两两两两两两两两两图9加入高斯白噪声的两路信号波形匹配滤波器在数字传输系统中，滤波器是不可缺少的。滤波器的一个作用是使基带信号频谱成形，例如为了满足奈奎斯特第一准则，基带信号频谱通常采用升余弦滚降形状，这一点在波形成形部分已经有了较详细的介绍。而滤波器的另一个重要作用是在接收端限制白噪声，将信号频带外的噪声滤掉，减少它对信号正确判决的影响。为了能够使滤波器输出信噪比在信息抽样时刻的信噪比最大，所以引入了匹配滤波器。假设匹配滤波器的频率传递函数为，时域冲激响应为。滤波器输入为发送信号FHTH与噪声的叠加，即TNSTX这里，为信号，它的频谱函数为。为白色高斯噪声，其双边功率谱密度为FSTN。滤波器的输出为20THNTSY其中信号部分为DFEFHSTTTJS2在时刻输出的信号抽样值为TFEFYTJS2滤波器输出噪声的功率谱密度为2|0FHFFNN平均功率为DFFNN20|0因此，时刻的输出信噪比为TTDFHFFESRNFTJ2|匹配滤波器的传递函数使达到最大。在这里利用SCHWARTZ不等式求解，最后得到传SNR递函数的表达式为FTEKSF2即传递函数与信号频谱的复共轭成正比。传递函数的时域响应为TTH匹配滤波器的最大输出信噪比为002/|NEDFSNRS其中，为观察间隔内的信号能量。SE具体到这个通信系统中，由于信号的时域响应为0TTHTTS其中是平方根升余弦滤波器的冲激响应。0H结合上式可以得到0TKT匹配滤波器实质上是一个具有与发射端的基带成形滤波器相同的滚降系数的平方根升余弦滤波器。接收端的“匹配滤波”是针对发射端的成形滤波而言，与成形滤波相匹配实现了数字通信系统的最佳接收。它与基带成形滤波器共同构成了一个奈奎斯特滤波器。源代码同平方根升余弦滤波器的源代码相同。信号通过匹配滤波器后的波形如图10所示0102030405060708042024两两两两两两0102030405060708042024两两两两两两两两两两两两两两两两两两两两图10经过匹配滤波器后的波形从上边的波形可以看出来，经过匹配滤波器后的信号明显很平滑，这正好反映了低通滤波器的特性，滤掉了高频分量，为了明显的反映这一特点，将一段高斯白噪声经过匹配滤波器。波形对比如图11所示。01020304050607080901003210123两两两两两两两两两两两两0102030405060708090100210123两两两两两两两两两两两两两两两两图11高斯白噪声经过匹配滤波器后的波形经过仿真，发现高斯白噪声经过一个平方根滤波器后方差保持不变。因此在加入高斯噪声时给定的信噪比需要有一定的修正，即要保证在信息点上的信噪比为给定的值。我们知道，当给定SNR时，需要加入的高斯噪声的功率谱密度，其LG10NSNRSNSNR100中，在具体使用AWGN函数时，1043132222SNR值是建立在当前输入信号的平均功率的基础上的，所以S。10LGLG1010LGL10SSNRSNRNSSRSNR采样由于从匹配滤波器出来的信号的点数8倍于原来信息的点数，为了恢复出原信号，所以需要对该信号进行采样。从匹配滤波器出来时，首先要剔除卷积过程中冗余的点，接着抽取现在信号中的第1个，第9个，第8K1个点，源代码如下FUNCTIONY1,Y2PICK_SIGX1,X2,RATIOY1X1RATIO321RATIOLENGTHX1Y2X2RATIO321RATIOLENGTHX1将这时的数据画到星座图上。5432101234554321012345QAM两两两图12信噪比为10DB时的星座图判决解调经过前边的匹配滤波器解调或者称为相关解调产生了一组向量，在这里就是一个一维的向量，根据最大后验概率（MAP）准则（由于各个信号的先验概率相等，所以页可以认为是最大似然准则），得到了最小距离检测。具体在本仿真系统中，判断为各个信号的门限如表2所示。判决后得到的数据再按照格雷码的规则还原成0、1信号，最终将两路0、1信号合成一路0、1信号，用来同最初的信号一起决定误码率。表2判决电平对应表判决前的信号的幅度对应的判决后的幅度2A30113源代码如下FUNCTIONYDEMODULATE_SIGX1,X2对X1路信号进行判决XX1FINDX123XX1FINDX101XX1FINDX12XX2FINDX201XX2FINDX22均值SIGMA115153协方差矩阵RMVNRNDMU,SIGMA,100生成100个随机样本PLOTR,1,R,2,下面生成一个多元正态分布的例子。样本的散点图如右所示二元正态分布的散点图（二）多元正态分布密度函数MU11SIGMA09040403X21PMVNPDFX,MU,SIGMAP13828E005多元正态分布的密度函数是MVNPDF。调用方式MVNPDFX,MU,SIGMA下面是一个例子。MU11SIGMA09040403X21FMVNCDFX,MU,SIGMAF08541FX，YPXA1234A1234下面是一个例子。MEANAANS23MEANA,2ANS1500035000（二）剔除异常值后的平均值X样本观察矩阵。PERCENT剔除比率，例如PERCENT10表示同时剔除最大的5和最小的5观察值。DIMDIM1（默认）表示对每列求平均值，DIM2表示对每行求平均值。有时观察数据中有异常大或异常小的值，这些异常值会对平均值产生影响，需要去掉异常值。例如在体操比赛中，去掉一个最高分和最低分，然后计算运动员的最后得分。在MATLAB中也有剔除异常值后的平均数。调用方式MTRIMMEANX,PERCENTMTRIMMEANX,PERCENT,DIM输入参数XRAND1,20TRIMMEAN（X,10）ANS05145（三）计算中位数A样本观测矩阵DIMDIM1（默认）表示对每列求中位数，DIM2表示对每行求中位数剔除10的异常值之后的平均数为05145。调用方式MMEDIAN（A）MMEDIAN（A，DIM）输入参数有时数据中出现NAN，在计算中位数时需要忽略NAN，这时需要调用NANMEDIAN函数。（四）计算方差与标准差A样本值FLAG0（默认值）表示方差计算公式为1表示方差计算公式为一般说来，资产组合的风险越大，方差越大，波动性越大。方差由于其简单、直观以及良好的统计性质使其成为风险的代名词。在MATLAB中计算方差、标准差的函数分别是VAR、STD。方差调用方式VARAVARA，FLAG标准差调用方式STD（A）STD（A,FLAG）输入参数（五）计算样本的百分位数XRAND1,20PRCTILEX,20ANS01688调用方式YPRCTILE（X，P，DIM）输入参数X观察值P计算大于P值的数DIM同上下面是一个例子（六）计算样本极差RRANGEQRRANGEQ,DIM极差就是样本极大值与极小值的差，反映样本的离散程度。调用方式XRAND1,20RANGEXANS08404下面是一个例子（七）计算偏度与峰度方差作为风险的度量指标并不是完整的。比如讲，两种资产收益分布的均值和方差都是相同的，但是一种资产收益是左偏的，另一种是右偏的。对于风险而言，相对于大概率和小损失人们更加厌恶小概率大损失的情况，后一种情况给人们带来的痛苦远大于第一情况。从这个意义上讲，收益分布左偏的资产的风险水平要小于右偏的资产。此时，方差作为风险的度量指标就不是完全的，还要考虑峰度、偏度等指标。偏度和峰度是两个高阶的统计量。计算偏度的目的在于考察组合收益率水平是否是对称分布，也就是组合产生亏损与获得盈利的概率如何；峰度是考察组合的收益率情况是否接近正态分布，如果组合的收益率存在尖峰厚尾的分布特征，则说明组合产生亏损和盈利的概率偏大，也就是在一定程度上认为组合收益率出现极端性的可能性偏大，这种组合的收益率稳定性是比较差的。正态分布的峰度等于3，大于3表示尖峰，小于3表示分布比较均匀。股票市场收益率的时间序列大都为尖峰肥尾。偏度的计算公式为式中，分别为样本X的均值与方差。如果SKEWNESS0，则表示分布形态与正态分布偏度相同；如果SKEWNESS0，则表示正偏差数值较大，长尾巴拖在右边；如果SKEQNESSXRAND1,20SKEWNESSXANS004871计算偏度调用方式YSKEWNESSAYSKEWNESSA,FLAG输入参数A样本值FLAG偏度的计算方式，0（默认）为无偏估计，1为有偏估计下面是一个例子。KKURTOSISXKKURTOSISX,FLAGKKURTOSISX,FLAG,DIM2计算峰度调用方式X样本观察矩阵FLAG峰度的计算方式，0（默认）为无偏估计，1为有偏估计DIMDIM1（默认）表示对每列求平均，DIM2表示对每行求平均输入参数XRAND1,20KURTOSISXANS14407下面是一个例子。（八）计算绝对离差绝对离差是以偏差绝对数来衡量决策方案的风险。在期望值相同的情况下，绝对离差越大，风险越大；绝对离差越小，风险越小。调用方式YMADXYMADX,N输入参数X观察值N绝对偏差计算方式N0（默认）计算公式为MEANABSXMEANXN1计算公式为MEDIANABSXMEDIANXXRAND1,20MADXANS01750下面是一个绝对离差的例子。（九）计算中心矩数理统计中经常用到中心矩的概念，K阶中心矩的计算公式为可以看出1阶中心矩为0，如果观察值是矩阵则以每列为样本计算中心矩。X观察样本值ORDER中心矩的阶数，必须为正整数调用方式MMOMENTX,ORDER输入参数XRAND65X041540990103200043990133803050078890960109334020710874404387072660683306072001500498304120021260629907680021400744608392037050970806435026790628805751MMOMENTX,3M0011300014000320005800023下面计算样本的3阶矩。（十）计算协方差和相关系数协方差是一个用于衡量投资组合任意两个资产相关性的统计指标。当协方差为正值时，表示两种资产的收益率呈同方向变动；协方差为负值时，表示两种资产的收益率呈相反方向变化；协方差等于0时，表示两种资产不存在相关性。相关系数总是在11之间的范围内变动，1表示完全负相关（反向），1表示完全正相关（同向），0则表示不相关。MATLAB计算协方差、相关系数的函数分别是COV和CORRCOEF。协方差调用方式CCOVXCCOVX,Y下面是一个例子A112231403COVAANS103334166730000416672333315000300001500010000X观察值矩阵Y观察向量PARAM1参数1，参数的值如下ALPHA表示置信度，在01之间VAL1参数1的值PARAM2参数2VAL2参数2的值2相关系数调用方式RCORRCOEFXRCORRCOEFX,YR,PCORRCOEFX,PARAM1,VAL1,PARAM2,VAL2,输入参数R相关系数矩阵P每个相关系数的概率矩阵输出参数XRAND30,4X,4SUMX,2R,PCORRCOEFXI,JFINDPI,J下面是一个计算相关系数的例子。R10000014120195404993014121000001312058480195401312100000372904993058480372910000P10000045660300800050045661000004896000070300804896100000042400050000070042410000ANS414243142434三、统计绘图X12357334X10000200003000050000700003000034000TABULATEXVALUECOUNTPERCENT1114292114293228573411429511429711429（一）样本频率分布图样本频率分布图函数是TABULATE。下面调用CDFPLOT函数绘出X的分布图。CDFPLOTX向量X的分布图（二）最小二乘拟合图在MATLAB中绘制最小二乘拟合图的命令是LSLINE，下面是一个例子。XRAND1,20XCUMSUMXPLOTX,LSLINE最小二乘拟合图（三）正态分布概率图有时需要判断样本数据是否服从正态分布，NORMPLOT函数用图的形式给出直观的判断。如果数据点越靠近直线则分布越近似正态分布，越远则越偏离正态分布。XNORMRND0,1,20,1PLOTX,NORMPLOTX正态分布拟合图从图中可以看出，数据点基本上是直线，所以符合正态分布。如果判断数据是否服从WEIBULL分布则可以对生成的数据用WBLPLOT函数进行判断。下图给出了WEIBULL分布拟合的结果。从图中看出对于数据较小、较大的点偏离较大，数据不服从WEIBULL分布。WEIBULL分布拟合图（四）样本密度图RANDNSEED,0XRANDN1,20XCUMSUMXCAPAPLOTX,0,10在MATLAB中绘出样本数据的密度图函数为CAPAPLOT。样本的密度示意图如右图所示。样本的密度示意图（五）频率直方图Y观察值。如果Y是一个向量，则画出一个频率图如果Y是一个MP阶矩阵，则对Y每一列分别作频率图NBINS频率图分成NBINS等分的区间段，默认值为10。调用方式NHISTYNHISTY,NBINSN,XOUTHISTY,NBINS输入参数N样本落在区间段的频率XOUT区间断的刻度输出参数下面是一个例子。RANDNSEED,0XRANDN1,200HISTX频率直方图如右图所示HISTX,MINX03MAXX其中MINX03MAXX表示频率图X轴的刻度起点是X最小元素，终点是X中最大元素，刻度间隔03。规定刻度间隔的频率直方图如下图所示。如果在频率图的基础上加上正态分布拟合图，则可以用HISTFIT函数。RANDNSEED,0XRANDN1,20HISTFITX带有密度函数的频率直方图如右图所示。（六）盒图X样本观察值G各组的名称PARAM1参数1的名称VAL1参数1的值在MATLAB中绘制样本数据的盒图函数是BOXPLOT。调用方式BOXPLOTXBOXPLOTX,GBOXPLOTX,PARAM1,VAL1,PARAM2,VAL2,输入参数各参数的名称和内容如下表所列参数名称参数值NOTCHON生成有缺口的盒图，OFF生成矩形盒图SYMBOL线条类型，默认值为RORIENTATIONVERTICAL（默认值）为垂直型盒图，HORIZONTAL为水平型盒图WHISKER盒图须线的长度，默认值15四分位间距LABELS盒图行列的名称标签COLORS线条颜色WIDTHS盒图的宽度，默认值05POSITIONS盒图的位置，默认值为1NGROUPORDER组的次序X1NORMRND5,1,100,1X2NORMRND6,1,100,1BOXPLOTX1,X2,NOTCH,ON盒子的上下两条线分别为样本的25和75分位数，盒子的上下底之间的距离为四分位的间距。盒子的中间线为样本中值，如果样本中值不在盒子的中间，表示存在一定的偏度。盒子的上方和下方有两条虚线，显示了样本的范围，野值（异常值）位于超过盒子顶端、底端15倍的四分位数。含有缺口的盒图中齿形缺口表示样本中值的置信区间。图中的内容说明如下正态分布盒图如右图所示。四、多元线性回归分析在金融上常常需要对金融、经济数据进行回归，其中最简单的是多元线性回归。（一）多元线性回归BREGRESSY,XB,BINTREGRESSY,XB,BINT,RREGRESSY,XB,BINT,R,RINTREGRESSY,XB,BINT,R,RINT,STATSREGRESSY,XB,BINT,R,RINT,STATSREGRESSY,X,ALPHA假设因变量Y和自变量X之间服从以下的线性模型式中Y是因变量的观察值，X是自变量回归矩阵，是参数向量，是白噪声。的最小二乘解是调用方式X自变量观察值，注意如果模型中有常数项，则X的第一列所有元素为1。Y因变量观察值向量ALPHA参数的置信度输入参数例1首先按照下面模型生成一系列随机数，然后回归。B的估计值，注意B中已经包含了常数项BINT的置信区间R残值RINT残值的置信区间STATSR2、F、概率P输出参数B,BINT,R,RINT,STATSREGRESSY,ONES10,1,X,005下面生成一组随机数X110Y0104XNORMRND0,01,1,10下面估计RINT0179401149014350176401625017370200501417017950171201046022900056902470025590040200667022380188901008STATS10E00300010228370000000000B0130303953BINT00120024870376204144R00323001650005600294000410062200950014800078500440从B的估计值可以得知常数项和一次项的系数分别为01303,03953。在005置信水平下常数项系数估计区间为0012002487，X的系数置信区间为0376204144。由于样本数量非常少，参数估计并不稳定。下图是残差及其置信区间图。RCOPLOTR,RINT（二）多元正态回归在MATLAB中MVNRMLE函数可以进行多元正态回归，假设YK为随机变量，其分布如下式中NG,G为多元正态分布。调用方式PARAMETERS,COVARIANCE,RESID,INFOMVNRMLEY,DESIGN,MAXITERATIONS,TOLPARAM,TOL0BJ,COVAR0PARAMETERS参数COVARIANCE协方差RESID残差INFO估计过程的相关系数Y观察值矩阵，YNK中N是样本的个数，K是资产的数目DESIGN自变量单元变量矩阵，如果Y只有一个资产时，DESIGN是一个矩阵，如果Y中的资产个数大于一个时，DESIGN是一个单元向量，每个元素都是一个矩阵。Y的K列和DESIGN第K个元素中的矩阵进行回归MAXITERATIONSTOLPARAMTOL0BJCOVARO输出参数输入参数（三）估计多元正态分布每个资产的标准差DATA观察值矩阵，YNK中N样本的个数，K是资产的数目DESIGN自变量单元变量矩阵，如果Y只有一个资产时，DESIGN是一个单元向量，如果Y含有多于一个资产时，DESIGN是一个单元变量矩阵COVARIANCE回归时的残值调用方式STDPARAMETERS,STDCOVARIANCEMVNRSTDDATA,DESIGN,COVARIANCE输入参数STDPARAMETERS每个资产的标准差STDCOVARIANCE协方差输出参数（四）岭回归线性回归中参数估计，如果观察值X存在自相关性，则XTX是奇异矩阵，估计值就会出现非常大的误差，这时矩阵XTX需要加上一个对角元素是常数K的单位阵，即。MATLAB提供了岭回归RIDGE函数求解该问题。B模型估计参数调用方式B1RIDGEY,X,KB0RIDGEY,X,K，0输入参数Y因变量观察值X自变量观察值KK表示控制系数，可以根据需要进行选择。输出参数K00011BRIDGEHEAT,INGREDIENTS,KPLOTK,BXLABELRIDGEPARAMETERYLABELSTANDARDIZEDCOEFTITLERIDGETRACEFORHALDDATALEGENDX1,X2,X3,X4例2对HALD文件中的数据进行岭回归。LOADHALD查看工作区中的变量。WHOYOURVARIABLESAREHALDHEATINGREDIENTS五、主成分分析主成分分析是在各个变量之间相关关系研究的基础上，用较少的新变量代替原来较多的变量，而且使这些较少的新变量尽可能多地保留原来较多的变量所反映的信息。（一）主成分分析基本原理首先对样本进行标准化处理，为简单起见，标准化后的样本仍记为X1,X2,X3,XP。设样本为X1,X2,X3,XP，其对应的样本均值为对应的标准差为S1，S2，S3,SP。设F1,F2,F3,FP是主成分，也即是X1,X2,X3,XP的线性表示，同时满足下面的条件主成分是原样本的正交变换。各主成分之间互不相关。主成分的总方差不变。主成分按方差从大到小排序。主成分具有如下性质这一性质说明，主成分是原变量的线性组合，是对原变量信息的一种改良；主成分不增加总信息量，也不减少总信息量。设为主成分的特征值，前K个方差累积贡献率为一般当累积贡献率大于85时不再增加新的主成分。保留多少个主成分取决于保留部分的累积方差在总方差中占的百分比（即累计贡献率），它标志着前几个主成分概括的信息的多寡。在实践中，粗略规定一个百分比就可以决定保留几个主成分，如果多留一个主成分，但累积方差增加无几，便不再多留。（二）主成分分析函数COEFF主成分系数SCORE新坐标系LATENTX的协方差矩阵的特征值TSQUAREHOTELLING统计量的值在MATLAB中提供了两个主成分分析函数PRINCOMP和PCACOV。COEFF,SCOREPRINCOMPXCOEFF,SCORE,LATENTPRINCOMPXCOEFF,SCORE,LATENT,TSQUAREPRINCOMPX输入参数X观察变量输出参数调用方式CORRCOEFINGREDIENTSANS10000022860824102454022861000001392097300824101392100000029502454097300029510000例3用MATLAB自带数据进行主成分分析。MATLAB中自带了数据文件HALD，可以直接调用进行主成分分析。HALD数据考虑影响温度的4个因素，温度保存在HEAT变量中，4个因素的观察值保存在INGREDIENTS中。由于4个因素之间存在相关性，无法直接回归，为解决这个问题，首先进行主成分分析，生成四个主成分变量，主成分之间互不相关，而且和观察值的信息是同样的。第一步载入数据，考察变量之间的相关性。LOADHALD载入MATLAB自带的数据文件考察相关性发现自变量2与变量3之间的高度相关，所以需要剔除相关性。第二步主成分分析。PC,SCORE,LATENT,TSQUAREPRINCOMPINGREDIENTS主成分系数PC00678064600567305062067850020005440049330029007553040360515607309010850468404844主成分的方差贡献率SCORE36821868709459090396729607346109224760395812981842049090221126123714766341185470378605532446176087401424108125364660913001350325882897981606300818226064107259323650324392626898540016905437328401415737046503405922001238613428304352255849278170386704468269032293102445504116协方差的特征值LATENT517796967496412405402372TSQUARE56803