《数据分析方法》
课 程 实 验 报 告
1 实验内容
(1)掌握回归分析的思想和计算步骤;
(2)编写程序完成回归分析的计算,包括后续的显著性检验、残差分析、Box-Cox变换等内容。
2模型建立与求解(数据结构与算法描述)
Y与X1,X2,…,
的关系可表示为Y=f(X1,X2,…,)+ε,回归分析既是
利用Y与X1,X2,…
,的观测数据,并在误差项的某些假定下确定
f(X1,X2,…,)。利用统计推断方法对所确定的函数的合理性一级由此关系
所揭示的Y与X1,X2,…,的关系作分析。
用最小二乘法计算回归参数
总离差平方和
残差平方和
回归平方和
SST=SSE+SSR R为复相关系数
Box-Cox变换
对Y做如下变换:
通过最大似然方法确定问题转化为选择,使
达到最小,其中
3.实验数据为:
习题
2.4
误差方差估计为:4.9445e+004 复相关系数:0.9595 exersice2_6
得 b=-54.9877 4.7082 0.3393 R=0.9376 因此回归方程为:y=-54.9877+4.7082X1+0.3393X2 习题2.9
复相关系数:0.8202
4.程序代码清单:
2.4
function huigui(X,Y) D=xlsread('2.4.xls'); p=3;n=15; X=D(:,1:end-1); Y=D(:,end); a=ones(n,1); X=[a,X];
B=(inv(X'*X))*X'*Y H=X*inv(X'*X)*X'; y1=mean(Y);%ƽ¾ùÖµ
y2=B(1)+B(2).*X(:,2)+B(3).*X(:,3);%»Ø¹éÖµ SST=sum((Y-y1).^2)%×ÜÀë²îƽ·½ºÍ SSE=sum((Y-y2).^2)%²Ð²îƽ·½ºÍ SSR=sum((y2-y1).^2)%»Ø¹éƽ·½ºÍ MSR=SSR/(p-1)%¾ù·½»Ø¹é MSE=SSE/(n-p)%¾ù·½²Ð²î F=MSR/MSE;
R=sqrt(SSR/SST)%¼ÆË㸴Ïà¹ØϵÊý o=SSE/(n-p)%·½²î¹À¼ÆÖµ y3=y2-Y;%²Ð²îÖµ
e=(eye(n)-H)*Y;%²Ð²îÏòÁ¿ plot(Y,e,'*')
未进行Box-Cox变换时的残差分析
进行Box-Cox变换后,残差图为
2.6
load exersice2_6.txt data=exersice2_6; [b,R]=reg(data)
function [b,R]=reg(data) X=data(: , 1 : end-1); Y=data(: , end);
temp=ones(length(data),1); X=[temp,X]; b=inv (X'*X)*X'*Y; YY=X*b;
SST=sum((Y-mean(Y)).^2); SSE=sum((Y-YY).^2); R=sqrt(1-SSE/SST);
2.9
function huigui(X,Y) D=xlsread('2.9.xls'); p=4;n=23; X=D(:,1:end-1); Y=D(:,end); a=ones(n,1); X=[a,X];
B=(inv(X'*X))*X'*Y H=X*inv(X'*X)*X'; y1=mean(Y);%ƽ¾ùÖµ
y2=B(1)+B(2).*X(:,2)+B(3).*X(:,3)+B(4).*X(:,4);%»Ø¹éÖµ SST=sum((Y-y1).^2);%×ÜÀë²îƽ·½ºÍ SSE=sum((Y-y2).^2);%²Ð²îƽ·½ºÍ SSR=sum((y2-y1).^2);%»Ø¹éƽ·½ºÍ MSR=SSR/(p-1);%¾ù·½»Ø¹é MSE=SSE/(n-p);%¾ù·½²Ð²î F=MSR/MSE;
R=sqrt(SSR/SST)%¼ÆË㸴Ïà¹ØϵÊý o=SSE/(n-p);%·½²î¹À¼ÆÖµ y3=y2-Y;%²Ð²îÖµ
e=(eye(n)-H)*Y;%²Ð²îÏòÁ¿ plot(Y,e,'*')
残差分析
《数据分析方法》
课 程 实 验 报 告
1. 模型建立与求解(数据结构与算法描述)
主成分分析:将p个特征变量
通过线性组合综合成尽量小的几个
综合性变量(q
特性,又互不线性相关。
典型相关分析:将两组变量间的相关性凝结为少数几对典型变量间的相关性,通过对相关性较大的少数几对典型变量的研究来了解原来的两组变量的相关性。
主成分分析
设来自于X的一个容量为n的样本观测数据
, i=1,2,…,n
这时,我们便可以用其样本的协方差矩阵S或样本系数相关矩阵R分别作为Σ
或的估计进行主成分分析。 S=
=
R=其中
=
=
, =,j=1,2,…,p
=, j,k=1,2,…,p
设S=为样本协方差矩阵,其特征值为,相
应的正交单位化特征向量为则第k个样本主成分可表示为
这里 (i=1,2,…,p)。
, k=1,2,…,p
其中
x=表示
X=的观测值。当依次代入观测
值
(i=1,2,…,n)时,便得到第k
个样本成分的n
个观测值
(i=1,2,…,n),我们称之为第k个样本主成分的得分。第k个样本主成分的贡
献率为 ,前m个样本主成分的累计贡献率为
典型相关分析
关于X和Y的n组观测数据:
,
,i=1,2,…,n
同主成分分析一样,利用这些观测数据的样本协方差矩阵 S=作为Σ的估计,其中
以S代替Σ所求得的典型变量和典型相关系数分别成为样本典型变 型相关系数。具体地,第k对样本典型相关变量为
=
样本典型相关系数
,k=1,2,…,p x,
=
y,k=1,2,…,p
其中为的特征值
(从而也是
的前p个最大特征值),,,…,和,,…,,分别为和的相应于特征
值的正交单位化特征向量;
x=和
y=分别代表X和Y的观测值,
为的正平方根。同样,我们也
可以求标准化样本的样本典型变量与样本典型相关系数。 典型相关系数的显著性检验 在总体
服从p+q维正态分布
(μ,Σ)条件下,对一般情况的第k
个假设,可用如下的似然比统计量进行检验。令
=-
其中
为各样本典型相关系数的平方,即矩阵
)的特征值,n为样本容量。
(或
采用的是一个在样本容量n较小时有更好逼近精度的渐近服从F分布的统计量 ,即
=
其中
=wt-(p-k+1)(q-k+1)+1
W=n-(p+q+3)
t=
当某个k值使得=2时,取t=1.当为真时,渐近服从
自由度为和的F分布,且大的值意味着应拒绝,故检验p值为
其中为
的观测值。
2. 实验数据与实验结果
习题4.5
相关系数矩阵R为
1.0000 -0.1362 0.0245 -0.2061 0.9445 -0.0655 -0.0535 -0.0466 -0.1362 1.0000 -0.1496 -0.0210 -0.1063 0.9177 -0.0956 -0.0387 0.0245 -0.1496 1.0000 -0.0564 0.0415 -0.0658 0.9398 -0.0044 -0.2061 -0.0210 -0.0564 1.0000 -0.1090 0.0422 0.0200 0.9011 0.9445 -0.1063 0.0415 -0.1090 1.0000 -0.0434 -0.0267 0.0770 -0.0655 0.9177 -0.0658 0.0422 -0.0434 1.0000 -0.0787 0.0652 -0.0535 -0.0956 0.9398 0.0200 -0.0267 -0.0787 1.0000 0.0329 -0.0466 -0.0387 -0.0044 0.9011 0.0770 0.0652 0.0329 1.0000
各主成分的贡献率为
0.2783 0.2531 0.2281 0.2080 0.0138 0.0098 0.0052 0.0037 前两个主成分的累计贡献率为 0.5314 按第一主成分将30个省、区、市排序
0.1 江西 0.0984 西藏 0.0919 湖北 0.0702 黑龙江 0.054 天津 0.0481 贵州 0.0442 海南 0.0436 青海 0.0396 安徽 0.0352 新疆 0.034 北京 0.0332 山东 0.0263 吉林 0.0232 河北 -0.037 山西 -0.04 陕西 -0.0417 福建
-0.0448 江苏 -0.046 云南 -0.0469 辽宁 -0.0484 河南 -0.0528 上海 -0.0551 广西 -0.0569 甘肃 -0.061 湖南 -0.0627 内蒙古 -0.0694 浙江 -0.0699 广东 -0.0702 宁夏 -0.0768 四川
习题4.8
典型相关变量分别为
0.7703 0.3822 -0.6377 0.9241
0.9239 -0.0071 -0.3826 1.0000 样本典型相关系数为
0.1577 0 0 0.0053
显著性检验
k[1]=1 s[1]=0.052772 F[1]=228.011886 d1[1]=4.000000 d2[1]=272.000000
k[2]=2 s[2]=0.004007 F[2]=34054.665184 d1[2]=1.000000
d2[2]=137.000000
所以第一对显著相关 习题4.10
典型变量分别为
-0.5301 0.7009 0.1326 0.6322 0.1540 0.8686 -0.5651 -0.6964 0.4775
0.8057 0.3569 -0.1847 -0.0440 -0.5409 0.9483 -0.5907 0.7616 0.2581 典型相关系数为
557.4797 0 0 0 9.2910 0
0 0 270.9025
显著性检验
k[1]=1 s[1]=3.582420 F[1]=-14.792337 d1[1]=9.000000d2[1]=326.271395
k[2]=2 s[2]=2.563813 F[2]=-25.343880 d1[2]=4.000000d2[2]=270.000000
k[3]=3 s[3]=1.339616 F[3]=-34.478358 d1[3]=1.000000d2[3]=136.000000
所以一,二对典型变量显著相关
4.程序代码清单:
习题4.5
%cwstd.m,用总和标准化法标准化矩阵 function std=cwstd(vector)
cwsum=sum(vector,1); %对列求和
[a,b]=size(vector); %矩阵大小,a为行数,b为列数 for i=1:a for j=1:b
std(i,j)= vector(i,j)/cwsum(j); end end %cwfac.m
function result=cwfac(vector); fprintf('相关系数矩阵:\n')
std=CORRCOEF(vector) %计算相关系数矩阵 fprintf('特征向量(vec)及特征值(val):\n')
[vec,val]=eig(std) %求特征值(val)及特征向量(vec) newval=diag(val) ;
[y,i]=sort(newval) ; %对特征根进行排序,y为排序结果,i为索引 fprintf('特征根排序:\n') for z=1:length(y)
newy(z)=y(length(y)+1-z); end
fprintf('%g\n',newy) rate=y/sum(y);
fprintf('\n贡献率:\n') newrate=newy/sum(newy) sumrate=0; newi=[];
for k=length(y):-1:1
sumrate=sumrate+rate(k); newi(length(y)+1-k)=i(k); if sumrate>0.85 break; end
end %记下累积贡献率大85%的特征值的序号放入newi中 fprintf('主成分数:%g\n\n',length(newi)); fprintf('主成分载荷:\n') for p=1:length(newi) for q=1:length(y)
result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p)); end
end %计算载荷 disp(result) %cwscore.m,计算得分
function score=cwscore(vector1,vector2); sco=vector1*vector2; csum=sum(sco,2);
[newcsum,i]=sort(-1*csum); [newi,j]=sort(i); fprintf('计算得分:\n')
score=[sco,csum,j]
%得分矩阵:sco为各主成分得分;csum为综合得分;j为排序结果 %cwprint.m
function print=cwprint(filename,a,b);
%filename为文本文件文件名,a为矩阵行数(样本数),b为矩阵列数(变量指标数) fid=fopen(filename,'r')
vector=fscanf(fid,'%g',[a b]); fprintf('标准化结果如下:\n') v1=cwstd(vector) result=cwfac(v1); cwscore(v1,result);
习题4.8
function DX n=140;
s11=[1.00 0.63;0.63 1.00]; s22=[1.00 0.42;0.42 1.00]; s12=[0.24 0.06;-0.06 0.07]; s21=[0.24 -0.06;0.06 0.07]; A=inv(s11)*s12*inv(s22)*s21; B=inv(s22)*s21*inv(s11)*s12; [vala,e]=eig(A) [valb,f]=eig(B) p2=diag(vala);
p1=sqrt(vala);p=2;q=2;temp=1; s=ones(1,p);d1=s;d2=s;t=s;F=s; w=140-0.5*(p+q+3); for k=1:p
for i=k:p
temp=temp*(1-p2(k)); end s(k)=temp;
d1(k)=(p-k+1)*(q-k+1);
t(k)=sqrt(((p-k+1)^2*(q-k+1)^2-4)/((p-k+1)^2+(q-k+1)^2-5)); d2(k)=w*t(k)-0.5*(p-k+1)*(q-k+1)+1;
F(k)=d2(k)/d1(k)*((1-s(k)^(1/t(k)))/s(k)^(1/t(k))); end for j=1:p
fprintf('k[%d]=%d\ts[%d]=%f\tF[%d]=%f\td1[%d]=%f\td2[%d]=%f\n',j,j,j,s(j),j,F(j),j,d1(j),j,d2(j)); end
习题4.10
function DX
s=load('exercise4_6.txt'); n=49;X=s(:,2:4);Y=s(:,5:7); s11=cov(X); s22=cov(Y); s12=zeros(3,3); for j=1:3 for k=1:3 sum=0; for i=1:n
sum=sum+(X(i,j)-mean(X(:,j)))*(Y(i,k)-mean(Y(:,k))); end
s12(j,k)=sum; end end s21=s12';
A=inv(s11)*s12*inv(s22)*s21; B=inv(s22)*s21*inv(s11)*s12; [vala,e]=eig(A) [valb,f]=eig(B) p2=diag(vala);
p1=sqrt(vala);p=3;q=3;temp=1; s=ones(1,p);d1=s;d2=s;t=s;F=s; w=140-0.5*(p+q+3); for k=1:p for i=k:p
temp=temp*(1-p2(k)); end
s(k)=temp;
d1(k)=(p-k+1)*(q-k+1);
t(k)=sqrt(((p-k+1)^2*(q-k+1)^2-4)/((p-k+1)^2+(q-k+1)^2-5)); d2(k)=w*t(k)-0.5*(p-k+1)*(q-k+1)+1;
F(k)=d2(k)/d1(k)*((1-s(k)^(1/t(k)))/s(k)^(1/t(k))); end for j=1:p
fprintf('k[%d]=%d\ts[%d]=%f\tF[%d]=%f\td1[%d]=%f\td2[%d]=%f\n',j,j,j,s(j),j,F(j),j,d1(j),j,d2(j));
end
《数据分析方法》
课 程 实 验 报 告
模型建立与求解(数据结构与算法描述)
假定对所研究的对象(总体)在抽样前已有一定的认识,常用先验分布来描述这种认识,然后,基于抽取的样本对先验认识作修正,得到后验分布,而各种统计推断均基于后验分布进行。将Bayes统计的思想用于判别分析,就得到Bayes判别。
设总体别为
,
的协方差矩阵相等且为Σ,这样得到两个正态总体的Bayes判
当
,
及Σ
未知时,它们分别由
,
得训练样本算得的均值
,
及
协方差矩阵的联合估计S(=)来估计,线性判别函数为
误判概率的计算
其中 W(x)=
=
,K=
,d=
3. 实验数据与实验结果
先验概率相等时:
1 2 3 4 5 6 7 8 9
样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体 11 样品属于1总体 12 样品属于1总体 13 样品属于2总体 14 样品属于2总体 15 样品属于2总体 16 样品属于2总体 17 样品属于2总体 18 样品属于2总体 19 样品属于2总体 20 样品属于2总体 21 样品属于2总体 22 样品属于2总体 23 样品属于2总体 24 样品属于2总体 25 样品属于2总体 26 样品属于2总体 27 样品属于2总体 28 样品属于2总体 29 样品属于1总体 30 样品属于2总体 31 样品属于2总体 32 样品属于2总体 33 样品属于2总体 34 样品属于2总体 35 样品属于2总体
误判率的回代估计为:0.057143 1 2 3 4 5 6 7 8 9
样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体 11 样品属于1总体 12 样品属于1总体 13 样品属于2总体 14 样品属于2总体 15 样品属于2总体 16 样品属于2总体 17 样品属于2总体
19 样品属于2总体
20 样品属于2总体
21 样品属于2总体
22 样品属于2总体
23 样品属于2总体
24 样品属于2总体
25 样品属于2总体
26 样品属于2总体
27 样品属于2总体
28 样品属于1总体
29 样品属于1总体
30 样品属于2总体
31 样品属于2总体
32 样品属于2总体
33 样品属于2总体
34 样品属于2总体
35 样品属于1总体
误判率的交叉确认估计为:0.114286
先验概率按比例分配时
1
2
3
4
5
6
7
8
9 样品属于2总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体
11 样品属于2总体
12 样品属于2总体
13 样品属于2总体
14 样品属于2总体
15 样品属于2总体
16 样品属于2总体
17 样品属于2总体
18 样品属于2总体
19 样品属于2总体
20 样品属于2总体
21 样品属于2总体
22 样品属于2总体
23 样品属于2总体
24 样品属于2总体
26 样品属于2总体
27 样品属于2总体
28 样品属于2总体
29 样品属于1总体
30 样品属于2总体
31 样品属于2总体
32 样品属于2总体
33 样品属于2总体
34 样品属于2总体
35 样品属于2总体
误判率的回代估计为:0.171429
1
2
3
4
5
6
7
8
9 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体
11 样品属于1总体
12 样品属于1总体
13 样品属于2总体
14 样品属于2总体
15 样品属于2总体
16 样品属于2总体
17 样品属于2总体
18 样品属于2总体
19 样品属于2总体
20 样品属于2总体
21 样品属于2总体
22 样品属于2总体
23 样品属于2总体
24 样品属于2总体
25 样品属于2总体
26 样品属于2总体
27 样品属于2总体
28 样品属于2总体
29 样品属于1总体
30 样品属于2总体
31 样品属于2总体
32 样品属于2总体
34 样品属于2总体
35 样品属于2总体
误判率的交叉确认估计为:0.057143
4.程序代码清单:
function panbie(X1,X2)
X=load('exercise5_4.txt');
X1=X(1:12,:);
X2=X(13:35,:);
x1=ones(1,7);x2=ones(1,7);i=0;
for i=1:7
x1(i)=mean(X1(:,i));
x2(i)=mean(X2(:,i));
end
S=cov(X);
a1=inv(S)*x1';a2=inv(S)*x2';
% b1=-0.5*x1*inv(S)*x1'+log(12/35);b2=-0.5*x2*inv(S)*x2'+log(23/35);i=0;cou
% nt1=0;count2=0;%按比例分配
b1=-0.5*x1*inv(S)*x1';b2=-0.5*x2*inv(S)*x2';i=0;count1=0;count2=0;%等概率分配
for j=1:12
W1=a1'*X(j,:)'+b1;
W2=a2'*X(j,:)'+b2;
if(W1>=W2)
i=1;
fprintf('%d\t样品属于%d总体\n',j,i);
else
i=2;
fprintf('%d\t样品属于%d总体\n',j,i);
count1=count1+1;
end
end
for k=1:23
W1=a1'*X2(k,:)'+b1;
W2=a2'*X2(k,:)'+b2;
if(W1>=W2)
i=1;
fprintf('%d\t样品属于%d总体\n',k+12,i);
count2=count2+1;
else
i=2;
fprintf('%d\t样品属于%d总体\n',k+12,i);
end
end
p=(count1+count2)/35;
fprintf('误判率的回代估计为:%f\n',p);
function panbie1()
X=load('exercise5_4.txt');
i=0;count1=0;count2=0;
for m=1:12 %对总体一的剔除
if m==1
X1=X(2:12,:);X5=X1;
end
if m==12
X1=X(1:11,:);X5=X1;
else
X1=X(1:m-1,:);X3=X(m+1:12,:);
X5=[X1;X3];
end
X2=X(13:35,:);
X4=[X5;X2];
x1=ones(1,7);x2=ones(1,7);i=0;
for i=1:7
x1(i)=mean(X5(:,i));
x2(i)=mean(X2(:,i));
end
S=cov(X4);
a1=inv(S)*x1';a2=inv(S)*x2';
% b1=-0.5*x1*inv(S)*x1'+log(12/35);b2=-0.5*x2*inv(S)*x2'+log(23/35);
b1=-0.5*x1*inv(S)*x1';b2=-0.5*x2*inv(S)*x2';
W1=a1'*X(m,:)'+b1;
W2=a2'*X(m,:)'+b2;
if(W1>=W2) %对训练样本进行分类
i=1;
fprintf('%d\t样品属于%d总体\n',m,i);
else
i=2;
fprintf('%d\t样品属于%d总体\n',m,i);
count1=count1+1;
end
end
for m=13:35 %对总体二的剔除
if m==13
X2=X(14:35,:);X5=X2;
end
if m==35
X2=X(13:34,:);X5=X2;
else
X2=X(13:m-1,:);X3=X(m+1:35,:);
X5=[X2;X3];
end
X1=X(1:12,:);
X4=[X1;X5];
x1=ones(1,7);x2=ones(1,7);
for i=1:7
x1(i)=mean(X1(:,i));
x2(i)=mean(X5(:,i));
end
S=cov(X4);
a1=inv(S)*x1';a2=inv(S)*x2';
% b1=-0.5*x1*inv(S)*x1'+log(12/35);b2=-0.5*x2*inv(S)*x2'+log(23/35);%按比例分配
b1=-0.5*x1*inv(S)*x1';b2=-0.5*x2*inv(S)*x2';%等概率分配
W1=a1'*X(m,:)'+b1;
W2=a2'*X(m,:)'+b2;
if(W1>=W2) %对训练样本进行分类
i=1;
fprintf('%d\t样品属于%d总体\n',m,i);
count2=count2+1;
else
i=2;
fprintf('%d\t样品属于%d总体\n',m,i);
end
end
p=(count1+count2)/35;
《数据分析方法》
课 程 实 验 报 告
1 实验内容
(1)掌握回归分析的思想和计算步骤;
(2)编写程序完成回归分析的计算,包括后续的显著性检验、残差分析、Box-Cox变换等内容。
2模型建立与求解(数据结构与算法描述)
Y与X1,X2,…,
的关系可表示为Y=f(X1,X2,…,)+ε,回归分析既是
利用Y与X1,X2,…
,的观测数据,并在误差项的某些假定下确定
f(X1,X2,…,)。利用统计推断方法对所确定的函数的合理性一级由此关系
所揭示的Y与X1,X2,…,的关系作分析。
用最小二乘法计算回归参数
总离差平方和
残差平方和
回归平方和
SST=SSE+SSR R为复相关系数
Box-Cox变换
对Y做如下变换:
通过最大似然方法确定问题转化为选择,使
达到最小,其中
3.实验数据为:
习题
2.4
误差方差估计为:4.9445e+004 复相关系数:0.9595 exersice2_6
得 b=-54.9877 4.7082 0.3393 R=0.9376 因此回归方程为:y=-54.9877+4.7082X1+0.3393X2 习题2.9
复相关系数:0.8202
4.程序代码清单:
2.4
function huigui(X,Y) D=xlsread('2.4.xls'); p=3;n=15; X=D(:,1:end-1); Y=D(:,end); a=ones(n,1); X=[a,X];
B=(inv(X'*X))*X'*Y H=X*inv(X'*X)*X'; y1=mean(Y);%ƽ¾ùÖµ
y2=B(1)+B(2).*X(:,2)+B(3).*X(:,3);%»Ø¹éÖµ SST=sum((Y-y1).^2)%×ÜÀë²îƽ·½ºÍ SSE=sum((Y-y2).^2)%²Ð²îƽ·½ºÍ SSR=sum((y2-y1).^2)%»Ø¹éƽ·½ºÍ MSR=SSR/(p-1)%¾ù·½»Ø¹é MSE=SSE/(n-p)%¾ù·½²Ð²î F=MSR/MSE;
R=sqrt(SSR/SST)%¼ÆË㸴Ïà¹ØϵÊý o=SSE/(n-p)%·½²î¹À¼ÆÖµ y3=y2-Y;%²Ð²îÖµ
e=(eye(n)-H)*Y;%²Ð²îÏòÁ¿ plot(Y,e,'*')
未进行Box-Cox变换时的残差分析
进行Box-Cox变换后,残差图为
2.6
load exersice2_6.txt data=exersice2_6; [b,R]=reg(data)
function [b,R]=reg(data) X=data(: , 1 : end-1); Y=data(: , end);
temp=ones(length(data),1); X=[temp,X]; b=inv (X'*X)*X'*Y; YY=X*b;
SST=sum((Y-mean(Y)).^2); SSE=sum((Y-YY).^2); R=sqrt(1-SSE/SST);
2.9
function huigui(X,Y) D=xlsread('2.9.xls'); p=4;n=23; X=D(:,1:end-1); Y=D(:,end); a=ones(n,1); X=[a,X];
B=(inv(X'*X))*X'*Y H=X*inv(X'*X)*X'; y1=mean(Y);%ƽ¾ùÖµ
y2=B(1)+B(2).*X(:,2)+B(3).*X(:,3)+B(4).*X(:,4);%»Ø¹éÖµ SST=sum((Y-y1).^2);%×ÜÀë²îƽ·½ºÍ SSE=sum((Y-y2).^2);%²Ð²îƽ·½ºÍ SSR=sum((y2-y1).^2);%»Ø¹éƽ·½ºÍ MSR=SSR/(p-1);%¾ù·½»Ø¹é MSE=SSE/(n-p);%¾ù·½²Ð²î F=MSR/MSE;
R=sqrt(SSR/SST)%¼ÆË㸴Ïà¹ØϵÊý o=SSE/(n-p);%·½²î¹À¼ÆÖµ y3=y2-Y;%²Ð²îÖµ
e=(eye(n)-H)*Y;%²Ð²îÏòÁ¿ plot(Y,e,'*')
残差分析
《数据分析方法》
课 程 实 验 报 告
1. 模型建立与求解(数据结构与算法描述)
主成分分析:将p个特征变量
通过线性组合综合成尽量小的几个
综合性变量(q
特性,又互不线性相关。
典型相关分析:将两组变量间的相关性凝结为少数几对典型变量间的相关性,通过对相关性较大的少数几对典型变量的研究来了解原来的两组变量的相关性。
主成分分析
设来自于X的一个容量为n的样本观测数据
, i=1,2,…,n
这时,我们便可以用其样本的协方差矩阵S或样本系数相关矩阵R分别作为Σ
或的估计进行主成分分析。 S=
=
R=其中
=
=
, =,j=1,2,…,p
=, j,k=1,2,…,p
设S=为样本协方差矩阵,其特征值为,相
应的正交单位化特征向量为则第k个样本主成分可表示为
这里 (i=1,2,…,p)。
, k=1,2,…,p
其中
x=表示
X=的观测值。当依次代入观测
值
(i=1,2,…,n)时,便得到第k
个样本成分的n
个观测值
(i=1,2,…,n),我们称之为第k个样本主成分的得分。第k个样本主成分的贡
献率为 ,前m个样本主成分的累计贡献率为
典型相关分析
关于X和Y的n组观测数据:
,
,i=1,2,…,n
同主成分分析一样,利用这些观测数据的样本协方差矩阵 S=作为Σ的估计,其中
以S代替Σ所求得的典型变量和典型相关系数分别成为样本典型变 型相关系数。具体地,第k对样本典型相关变量为
=
样本典型相关系数
,k=1,2,…,p x,
=
y,k=1,2,…,p
其中为的特征值
(从而也是
的前p个最大特征值),,,…,和,,…,,分别为和的相应于特征
值的正交单位化特征向量;
x=和
y=分别代表X和Y的观测值,
为的正平方根。同样,我们也
可以求标准化样本的样本典型变量与样本典型相关系数。 典型相关系数的显著性检验 在总体
服从p+q维正态分布
(μ,Σ)条件下,对一般情况的第k
个假设,可用如下的似然比统计量进行检验。令
=-
其中
为各样本典型相关系数的平方,即矩阵
)的特征值,n为样本容量。
(或
采用的是一个在样本容量n较小时有更好逼近精度的渐近服从F分布的统计量 ,即
=
其中
=wt-(p-k+1)(q-k+1)+1
W=n-(p+q+3)
t=
当某个k值使得=2时,取t=1.当为真时,渐近服从
自由度为和的F分布,且大的值意味着应拒绝,故检验p值为
其中为
的观测值。
2. 实验数据与实验结果
习题4.5
相关系数矩阵R为
1.0000 -0.1362 0.0245 -0.2061 0.9445 -0.0655 -0.0535 -0.0466 -0.1362 1.0000 -0.1496 -0.0210 -0.1063 0.9177 -0.0956 -0.0387 0.0245 -0.1496 1.0000 -0.0564 0.0415 -0.0658 0.9398 -0.0044 -0.2061 -0.0210 -0.0564 1.0000 -0.1090 0.0422 0.0200 0.9011 0.9445 -0.1063 0.0415 -0.1090 1.0000 -0.0434 -0.0267 0.0770 -0.0655 0.9177 -0.0658 0.0422 -0.0434 1.0000 -0.0787 0.0652 -0.0535 -0.0956 0.9398 0.0200 -0.0267 -0.0787 1.0000 0.0329 -0.0466 -0.0387 -0.0044 0.9011 0.0770 0.0652 0.0329 1.0000
各主成分的贡献率为
0.2783 0.2531 0.2281 0.2080 0.0138 0.0098 0.0052 0.0037 前两个主成分的累计贡献率为 0.5314 按第一主成分将30个省、区、市排序
0.1 江西 0.0984 西藏 0.0919 湖北 0.0702 黑龙江 0.054 天津 0.0481 贵州 0.0442 海南 0.0436 青海 0.0396 安徽 0.0352 新疆 0.034 北京 0.0332 山东 0.0263 吉林 0.0232 河北 -0.037 山西 -0.04 陕西 -0.0417 福建
-0.0448 江苏 -0.046 云南 -0.0469 辽宁 -0.0484 河南 -0.0528 上海 -0.0551 广西 -0.0569 甘肃 -0.061 湖南 -0.0627 内蒙古 -0.0694 浙江 -0.0699 广东 -0.0702 宁夏 -0.0768 四川
习题4.8
典型相关变量分别为
0.7703 0.3822 -0.6377 0.9241
0.9239 -0.0071 -0.3826 1.0000 样本典型相关系数为
0.1577 0 0 0.0053
显著性检验
k[1]=1 s[1]=0.052772 F[1]=228.011886 d1[1]=4.000000 d2[1]=272.000000
k[2]=2 s[2]=0.004007 F[2]=34054.665184 d1[2]=1.000000
d2[2]=137.000000
所以第一对显著相关 习题4.10
典型变量分别为
-0.5301 0.7009 0.1326 0.6322 0.1540 0.8686 -0.5651 -0.6964 0.4775
0.8057 0.3569 -0.1847 -0.0440 -0.5409 0.9483 -0.5907 0.7616 0.2581 典型相关系数为
557.4797 0 0 0 9.2910 0
0 0 270.9025
显著性检验
k[1]=1 s[1]=3.582420 F[1]=-14.792337 d1[1]=9.000000d2[1]=326.271395
k[2]=2 s[2]=2.563813 F[2]=-25.343880 d1[2]=4.000000d2[2]=270.000000
k[3]=3 s[3]=1.339616 F[3]=-34.478358 d1[3]=1.000000d2[3]=136.000000
所以一,二对典型变量显著相关
4.程序代码清单:
习题4.5
%cwstd.m,用总和标准化法标准化矩阵 function std=cwstd(vector)
cwsum=sum(vector,1); %对列求和
[a,b]=size(vector); %矩阵大小,a为行数,b为列数 for i=1:a for j=1:b
std(i,j)= vector(i,j)/cwsum(j); end end %cwfac.m
function result=cwfac(vector); fprintf('相关系数矩阵:\n')
std=CORRCOEF(vector) %计算相关系数矩阵 fprintf('特征向量(vec)及特征值(val):\n')
[vec,val]=eig(std) %求特征值(val)及特征向量(vec) newval=diag(val) ;
[y,i]=sort(newval) ; %对特征根进行排序,y为排序结果,i为索引 fprintf('特征根排序:\n') for z=1:length(y)
newy(z)=y(length(y)+1-z); end
fprintf('%g\n',newy) rate=y/sum(y);
fprintf('\n贡献率:\n') newrate=newy/sum(newy) sumrate=0; newi=[];
for k=length(y):-1:1
sumrate=sumrate+rate(k); newi(length(y)+1-k)=i(k); if sumrate>0.85 break; end
end %记下累积贡献率大85%的特征值的序号放入newi中 fprintf('主成分数:%g\n\n',length(newi)); fprintf('主成分载荷:\n') for p=1:length(newi) for q=1:length(y)
result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p)); end
end %计算载荷 disp(result) %cwscore.m,计算得分
function score=cwscore(vector1,vector2); sco=vector1*vector2; csum=sum(sco,2);
[newcsum,i]=sort(-1*csum); [newi,j]=sort(i); fprintf('计算得分:\n')
score=[sco,csum,j]
%得分矩阵:sco为各主成分得分;csum为综合得分;j为排序结果 %cwprint.m
function print=cwprint(filename,a,b);
%filename为文本文件文件名,a为矩阵行数(样本数),b为矩阵列数(变量指标数) fid=fopen(filename,'r')
vector=fscanf(fid,'%g',[a b]); fprintf('标准化结果如下:\n') v1=cwstd(vector) result=cwfac(v1); cwscore(v1,result);
习题4.8
function DX n=140;
s11=[1.00 0.63;0.63 1.00]; s22=[1.00 0.42;0.42 1.00]; s12=[0.24 0.06;-0.06 0.07]; s21=[0.24 -0.06;0.06 0.07]; A=inv(s11)*s12*inv(s22)*s21; B=inv(s22)*s21*inv(s11)*s12; [vala,e]=eig(A) [valb,f]=eig(B) p2=diag(vala);
p1=sqrt(vala);p=2;q=2;temp=1; s=ones(1,p);d1=s;d2=s;t=s;F=s; w=140-0.5*(p+q+3); for k=1:p
for i=k:p
temp=temp*(1-p2(k)); end s(k)=temp;
d1(k)=(p-k+1)*(q-k+1);
t(k)=sqrt(((p-k+1)^2*(q-k+1)^2-4)/((p-k+1)^2+(q-k+1)^2-5)); d2(k)=w*t(k)-0.5*(p-k+1)*(q-k+1)+1;
F(k)=d2(k)/d1(k)*((1-s(k)^(1/t(k)))/s(k)^(1/t(k))); end for j=1:p
fprintf('k[%d]=%d\ts[%d]=%f\tF[%d]=%f\td1[%d]=%f\td2[%d]=%f\n',j,j,j,s(j),j,F(j),j,d1(j),j,d2(j)); end
习题4.10
function DX
s=load('exercise4_6.txt'); n=49;X=s(:,2:4);Y=s(:,5:7); s11=cov(X); s22=cov(Y); s12=zeros(3,3); for j=1:3 for k=1:3 sum=0; for i=1:n
sum=sum+(X(i,j)-mean(X(:,j)))*(Y(i,k)-mean(Y(:,k))); end
s12(j,k)=sum; end end s21=s12';
A=inv(s11)*s12*inv(s22)*s21; B=inv(s22)*s21*inv(s11)*s12; [vala,e]=eig(A) [valb,f]=eig(B) p2=diag(vala);
p1=sqrt(vala);p=3;q=3;temp=1; s=ones(1,p);d1=s;d2=s;t=s;F=s; w=140-0.5*(p+q+3); for k=1:p for i=k:p
temp=temp*(1-p2(k)); end
s(k)=temp;
d1(k)=(p-k+1)*(q-k+1);
t(k)=sqrt(((p-k+1)^2*(q-k+1)^2-4)/((p-k+1)^2+(q-k+1)^2-5)); d2(k)=w*t(k)-0.5*(p-k+1)*(q-k+1)+1;
F(k)=d2(k)/d1(k)*((1-s(k)^(1/t(k)))/s(k)^(1/t(k))); end for j=1:p
fprintf('k[%d]=%d\ts[%d]=%f\tF[%d]=%f\td1[%d]=%f\td2[%d]=%f\n',j,j,j,s(j),j,F(j),j,d1(j),j,d2(j));
end
《数据分析方法》
课 程 实 验 报 告
模型建立与求解(数据结构与算法描述)
假定对所研究的对象(总体)在抽样前已有一定的认识,常用先验分布来描述这种认识,然后,基于抽取的样本对先验认识作修正,得到后验分布,而各种统计推断均基于后验分布进行。将Bayes统计的思想用于判别分析,就得到Bayes判别。
设总体别为
,
的协方差矩阵相等且为Σ,这样得到两个正态总体的Bayes判
当
,
及Σ
未知时,它们分别由
,
得训练样本算得的均值
,
及
协方差矩阵的联合估计S(=)来估计,线性判别函数为
误判概率的计算
其中 W(x)=
=
,K=
,d=
3. 实验数据与实验结果
先验概率相等时:
1 2 3 4 5 6 7 8 9
样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体 11 样品属于1总体 12 样品属于1总体 13 样品属于2总体 14 样品属于2总体 15 样品属于2总体 16 样品属于2总体 17 样品属于2总体 18 样品属于2总体 19 样品属于2总体 20 样品属于2总体 21 样品属于2总体 22 样品属于2总体 23 样品属于2总体 24 样品属于2总体 25 样品属于2总体 26 样品属于2总体 27 样品属于2总体 28 样品属于2总体 29 样品属于1总体 30 样品属于2总体 31 样品属于2总体 32 样品属于2总体 33 样品属于2总体 34 样品属于2总体 35 样品属于2总体
误判率的回代估计为:0.057143 1 2 3 4 5 6 7 8 9
样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体 11 样品属于1总体 12 样品属于1总体 13 样品属于2总体 14 样品属于2总体 15 样品属于2总体 16 样品属于2总体 17 样品属于2总体
19 样品属于2总体
20 样品属于2总体
21 样品属于2总体
22 样品属于2总体
23 样品属于2总体
24 样品属于2总体
25 样品属于2总体
26 样品属于2总体
27 样品属于2总体
28 样品属于1总体
29 样品属于1总体
30 样品属于2总体
31 样品属于2总体
32 样品属于2总体
33 样品属于2总体
34 样品属于2总体
35 样品属于1总体
误判率的交叉确认估计为:0.114286
先验概率按比例分配时
1
2
3
4
5
6
7
8
9 样品属于2总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体
11 样品属于2总体
12 样品属于2总体
13 样品属于2总体
14 样品属于2总体
15 样品属于2总体
16 样品属于2总体
17 样品属于2总体
18 样品属于2总体
19 样品属于2总体
20 样品属于2总体
21 样品属于2总体
22 样品属于2总体
23 样品属于2总体
24 样品属于2总体
26 样品属于2总体
27 样品属于2总体
28 样品属于2总体
29 样品属于1总体
30 样品属于2总体
31 样品属于2总体
32 样品属于2总体
33 样品属于2总体
34 样品属于2总体
35 样品属于2总体
误判率的回代估计为:0.171429
1
2
3
4
5
6
7
8
9 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于1总体 样品属于2总体
10 样品属于1总体
11 样品属于1总体
12 样品属于1总体
13 样品属于2总体
14 样品属于2总体
15 样品属于2总体
16 样品属于2总体
17 样品属于2总体
18 样品属于2总体
19 样品属于2总体
20 样品属于2总体
21 样品属于2总体
22 样品属于2总体
23 样品属于2总体
24 样品属于2总体
25 样品属于2总体
26 样品属于2总体
27 样品属于2总体
28 样品属于2总体
29 样品属于1总体
30 样品属于2总体
31 样品属于2总体
32 样品属于2总体
34 样品属于2总体
35 样品属于2总体
误判率的交叉确认估计为:0.057143
4.程序代码清单:
function panbie(X1,X2)
X=load('exercise5_4.txt');
X1=X(1:12,:);
X2=X(13:35,:);
x1=ones(1,7);x2=ones(1,7);i=0;
for i=1:7
x1(i)=mean(X1(:,i));
x2(i)=mean(X2(:,i));
end
S=cov(X);
a1=inv(S)*x1';a2=inv(S)*x2';
% b1=-0.5*x1*inv(S)*x1'+log(12/35);b2=-0.5*x2*inv(S)*x2'+log(23/35);i=0;cou
% nt1=0;count2=0;%按比例分配
b1=-0.5*x1*inv(S)*x1';b2=-0.5*x2*inv(S)*x2';i=0;count1=0;count2=0;%等概率分配
for j=1:12
W1=a1'*X(j,:)'+b1;
W2=a2'*X(j,:)'+b2;
if(W1>=W2)
i=1;
fprintf('%d\t样品属于%d总体\n',j,i);
else
i=2;
fprintf('%d\t样品属于%d总体\n',j,i);
count1=count1+1;
end
end
for k=1:23
W1=a1'*X2(k,:)'+b1;
W2=a2'*X2(k,:)'+b2;
if(W1>=W2)
i=1;
fprintf('%d\t样品属于%d总体\n',k+12,i);
count2=count2+1;
else
i=2;
fprintf('%d\t样品属于%d总体\n',k+12,i);
end
end
p=(count1+count2)/35;
fprintf('误判率的回代估计为:%f\n',p);
function panbie1()
X=load('exercise5_4.txt');
i=0;count1=0;count2=0;
for m=1:12 %对总体一的剔除
if m==1
X1=X(2:12,:);X5=X1;
end
if m==12
X1=X(1:11,:);X5=X1;
else
X1=X(1:m-1,:);X3=X(m+1:12,:);
X5=[X1;X3];
end
X2=X(13:35,:);
X4=[X5;X2];
x1=ones(1,7);x2=ones(1,7);i=0;
for i=1:7
x1(i)=mean(X5(:,i));
x2(i)=mean(X2(:,i));
end
S=cov(X4);
a1=inv(S)*x1';a2=inv(S)*x2';
% b1=-0.5*x1*inv(S)*x1'+log(12/35);b2=-0.5*x2*inv(S)*x2'+log(23/35);
b1=-0.5*x1*inv(S)*x1';b2=-0.5*x2*inv(S)*x2';
W1=a1'*X(m,:)'+b1;
W2=a2'*X(m,:)'+b2;
if(W1>=W2) %对训练样本进行分类
i=1;
fprintf('%d\t样品属于%d总体\n',m,i);
else
i=2;
fprintf('%d\t样品属于%d总体\n',m,i);
count1=count1+1;
end
end
for m=13:35 %对总体二的剔除
if m==13
X2=X(14:35,:);X5=X2;
end
if m==35
X2=X(13:34,:);X5=X2;
else
X2=X(13:m-1,:);X3=X(m+1:35,:);
X5=[X2;X3];
end
X1=X(1:12,:);
X4=[X1;X5];
x1=ones(1,7);x2=ones(1,7);
for i=1:7
x1(i)=mean(X1(:,i));
x2(i)=mean(X5(:,i));
end
S=cov(X4);
a1=inv(S)*x1';a2=inv(S)*x2';
% b1=-0.5*x1*inv(S)*x1'+log(12/35);b2=-0.5*x2*inv(S)*x2'+log(23/35);%按比例分配
b1=-0.5*x1*inv(S)*x1';b2=-0.5*x2*inv(S)*x2';%等概率分配
W1=a1'*X(m,:)'+b1;
W2=a2'*X(m,:)'+b2;
if(W1>=W2) %对训练样本进行分类
i=1;
fprintf('%d\t样品属于%d总体\n',m,i);
count2=count2+1;
else
i=2;
fprintf('%d\t样品属于%d总体\n',m,i);
end
end
p=(count1+count2)/35;