数据分析方法实验

《数据分析方法》

课 程 实 验 报 告

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;


相关文章

  • 天然药物化学设计性镐药剂09实验方案
  • 天然药物化学设计性实验方案 生命科学与工程学院药物制剂专业 任课教师:李楠 孔阳 天然药物化学是药学专业的一门专业基础课,它涉及到有机化学.中药药剂学.生药学.波谱学等多个学科的内容,是运用化学.生药学的原理和方法来研究天然植物药的有效化学 ...查看


  • 塞曼效应实验数据分析与处理方法改进
  • 第30卷 2010年5月 第5期 物理实验 V01.30No.5 PHYSICSEXPERlMENTATION May,2010 塞曼效应实验数据分析与处理方法改进 杨 冰,丁 蔻,李丽华,董瑞新,闫循领 (聊城大学物理科学与信息工程学院, ...查看


  • 试剂分析性能评估模板
  • 胆固醇测定试剂盒 分析性能评估资料 山东高密彩虹分析仪器有限公司 目 录 1 概述 2 胆固醇测定试剂盒及相关信息 3 性能评估资料 3.1 检测限评估资料 3.2 线性范围评估资料 3.3 可报告范围评估资料 3.4 准确性(回收实验)评 ...查看


  • 实验数据的统计分析与科研论文的撰写
  • 2009年1月第19卷 第1期 中国比较医学杂志 CHINESE J OURNAL OF COMPAR ATIVE MEDICINE January , 2009Vol . 19 No . 1 知识讲座 实验数据的统计分析与科研论文的撰写 ...查看


  • 药理学实验报告.doc
  • 巢湖职业技术学院 ChaoHu vocational and Technical College 药理学实验报告 (供药学专业使用) 班级____________ 姓名____________ 组别____________ 指导教师____ ...查看


  • 无机及分析化学实验教学大纲
  • <无机及分析化学实验>教学大纲 [课程编号] [学时学分]120 学时: [开课模式]必修 [实验学时]36 学时 [上机学时] [课程类型] 专业课 [考核方式]考试 [先修课程] [开课单位]石油化工学院 [授课对象]石油化 ...查看


  • 教育学研究方法
  • 教育研究方法 一.教育研究概述  教育研究的界说 1.教育研究的含义:以发现或发展科学知识体系为导向,通过对教育现象的解释.预测和控制,以促进一般化原理.原则的发展: 2.教育研究的意义:①探索教育规律,以解决重要的教育理论与实践问题为导 ...查看


  • 教育技术研究方法
  • <教育技术学研究方法> 复习题答案及评分标准 一.名词解释(每个5分,必须简要说明其内含) 1.教育技术研究方法:P2就是人们为深刻认识应用教育技术进行教育活动的过程和现象而采用的途径.手段和工具等(3分).教育技术学研究方法是 ...查看


  • 中学物理实验教学研究
  • 中学物理实验教学研究 第一章 物理实验与物理学发展 第一节.物理实验在物理学发展中的地位. 物理学是一门实验科学,在物理学中,每个概念的建立.每个规律的发现,都有其坚实的实验基础.可以说没有实验就没有物理学. 一.物理实验与古代物理学发展. ...查看


  • 密立根油滴实验的几种数据处理方法的比较
  • 墨理工科研 密立根油滴实验的几种数据处理方法的比较 芦明霞李斌王天会王 (吉林大学珠海学院公共基础与教学研究中心 丹 广东・珠海519041) 中图分类号:0423.9 文献标识码:A 摘要比较了目前常用的几种密立根油滴实验数据处理的方法, ...查看


热门内容