mxsb.m

统计学习的Matlab程序:包括各种假设检验,模型识别和线性回归等.

function [p,Q]=mxsb(Z)
%MXSB 模型识别——确定线性模型的类别和阶数
%《随机过程》(汪荣鑫编)P132.
%Z为时间序列的向量,p为自相关系数向量，Q为偏相关系数向量
%蓝色线是自相关系数向量p的图像；红色线是偏相关系数向量Q的图像

%By：Ji Lin
%Email: [e-mail]
%Blog: http://linji526.spaces.live.com
%Date: 2007，2008/08/18

if nargin < 1
    Z=[15600,8960,10400,10600,10800,9880,9850,10900,8810,9960,12200,7510,...
        8640,6380,6810,8820,14400,7440,7240,6430,11000,7340,9260,5290,9130,...
        7480,6980,9650,7260,8750,9900,7310,9040,7310,8850,7840,10700,6190,...
        9610,7580,9990,6150,8250,6030,8980,6180,9630,9490,2340,11100,5090,...
        10900,6490,12600,6640,7430,6760,10000,9300];
end


W=Z-mean(Z);
n=length(Z);K=15;
r=ones(1,K);p=ones(1,K);
r0=W(1:n)*W(1:n)'/n;
for k=1:K
    r(k)=W(1:n-k)*W(1+k:n)'/n;
    p(k)=r(k)/r0;
end
q=ones(K,K);Q=ones(1,K);
q(1,1)=p(1);Q(1)=q(1,1);
for k=1:K-1
    q(k+1,k+1)=(p(k+1)-p(k:-1:1)*q(k,1:k)')/(1-p(1:k)*q(k,1:k)');
    for j=1:k
        q(k+1,j)=q(k,j)-q(k+1,k+1)*q(k,k-j+1);        
    end
    Q(k+1)=q(k+1,k+1);
end
plot(0:K,[1,p])
hold on
plot(0:K,[1,Q],'r:p')
legend('p(k)','Q(k)')
title('模型识别——确定线性模型的类别和阶数')
xlabel('k')
ylabel('p(k),Q(k)')
grid on