最新MATLAB在数据统计中的应用.docx
《最新MATLAB在数据统计中的应用.docx》由会员分享，可在线阅读，更多相关《最新MATLAB在数据统计中的应用.docx（15页珍藏版）》请在冰豆网上搜索。
最新MATLAB在数据统计中的应用
MATLAB在数据统计中的应用
MATLAB在数据统计中的应用
______________________________________________
目录：
1、一元线性回归的matlab实现（含检验）【更新】
2、一维数据滑动平均的matlab实现
3、多元线性回归的matlab实现
4、K阶自回归拟合及二阶自回归预测的Matlab实现
5、一次指数平滑预测的matlab实现
6、n次指数平滑及其预测
7、一维数据移动平滑的matlab实现
8、K阶自相关系数的matlab实现（含置信度检验）
说明：
1.正文中命令部分可以直接在Matlab中运行，作者（Yangfd09）在MATLABR2009a（7.8.0.347）中运行通过。
2.限于作者水平问题，文中难免疏漏和错误，如蒙赐教，不胜感激！
3.原创作品，仅供学习交流之用，会有不定期更新。
一元线性回归的matlab实现（含检验）【更新】
%求一元线性回归方程
%数据要求：
两行。
第一行存放x的观察值，第二行存放y的观察值
%数据文件名：
data_yyhg.mat；变量名：
test
%loaddata_yyhg.mat
N=length（test（1,:
））;%注：
也可以用[M,N]=size（test）
%但不能用N=size（test（1,:
））
sx=0;sx2=0;sy=0;sy2=0;sxy=0;Lxy=0;Lyy=0;
fori=1:
N
sx=sx+test（1,i）;
sx2=sx2+test（1,i）^2;
sy=sy+test（2,i）;
sy2=sy2+test（2,i）^2;
sxy=sxy+test（1,i）*test（2,i）;
Lxy=Lxy+（test（1,i）-sum（test（1,:
））/N）*（test（2,i）-sum（test（2,:
）/N））;
Lyy=Lyy+（test（2,i）-sum（test（2,:
））/N）^2;
end
r=[N,sx;sx,sx2]\[sy;sxy];
a=r
（1）;b=r
（2）;
%F分布检验
U=b*Lxy;
Q=Lyy-U;
F=（N-2）*U/Q;
%拟合优度检验
x=test（1,:
）;y=a+b*x;eq=sum（test（2,:
））/N;
ssd=0;ssr=0;
fori=1:
N
ssd=ssd+（test（2,i）-y（i））^2;
ssr=ssr+（y（i）-eq）^2;
end
sst=ssd+ssr;
RR=ssr/sst;
%命令窗口中显示回归方程
str=[blanks（5）,'y=','（',num2str（a）,'）','+','（',num2str（b）,'）','*x'];
disp（''）
disp（'回归方程为:
'）
disp（str）
disp（'R^2拟合优度检验：
'）
strin=['R^2=',num2str（RR）];
disp（strin）
disp（'F-分布显著性检验：
'）
stri=['F计算值：
',num2str（F）,blanks（4）,'自由度：
f1=1,f2=',num2str（N-2）];
disp（stri）
disp（'注：
请对照F-分布表找到所需置信水平下的F临界值Fa，若F>Fa，则通过检验。
'）
%绘制x-y散点图和回归直线
yy=a+b*test（1,:
）;
plot（test（1,:
）,test（2,:
）,'r.'）,holdon
plot（test（1,:
）,yy,'b-'）,holdoff
title（str）
附（可以直接粘贴到.mat文件中）：
3.8
4
5.8
8
11.3
14.4
16.5
16.2
13.8
10.8
6.7
4.7
77.7
51.2
60.1
54.1
55.4
56.8
45
55.3
67.5
73.3
76.6
79.6
一维数据滑动平均的matlab实现
%滑动平均
%数据格式：
单行（按时间序列排序）
%数据文件名：
data_ph.mat,变量名：
test
loaddata_ph.mat
M=length（test）;
disp（'请输入单侧平滑点数（时距）'）
k=input（'（输入1对应于三点平滑，2对应五点平滑）：
'）;
y=zeros（1,M）;
if2*k+1<=M
fori=1:
M-2*k
forj=i:
i+2*k
y（i+k）=y（i+k）+test（j）;
end
y（i+k）=y（i+k）/（2*k+1）;
end
y（[1:
k,M-k+1:
M]）=NaN;
str=[int2str（k）,'点滑动平均结果如下：
'];
disp（str）
formatcompact
data=test,result=y
format
else
disp（'Error:
数据个数不足！
'）
end
附（直接复制到.mat文件中即可）：
（某城市1999-2004年用水量数据）
211.3
260.18
209.1
248.79
241
250
多元线性回归的matlab实现
%多元线性回归
%数据要求：
M行N列。
第一行为y的观察值，其余行分别为x1,x2,...,x（M-1）
%数据文件名：
data_dyhg.mat;变量名：
test
loaddata_dyhg.mat
[M,N]=size（test）;
Y=test（1,:
）';
X=cat（2,ones（N,1）,test（2:
end,:
）'）;
A=X'*X;
B=X'*Y;
b=A\B;
formatlongg
str=['系数按升幂排列为：
',num2str（b'）];
disp（str）
%F-分布检验
eqx=sum（test（2:
M,:
）,2）/N;eqy=sum（Y）/N;
Liy=zeros（1,M-1）;
Lyy=0;
fori=1:
N
Lyy=Lyy+（Y（i）-eqy）^2;
end
%求回归平方和（方法一）：
用y的拟合值
y=zeros（1,N）;
UU=0;
fori=1:
N
forj=1:
M
y（i）=y（i）+b（j）*X（i,j）;
end
UU=UU+（y（i）-eqy）^2;
end
Q0=Lyy-UU;
F0=（UU/（M-1））/（Q0/（N-M））;
%求回归平方和（方法二）：
用偏回归系数
fori=1:
M-1
forj=1:
N
Liy（i）=Liy（i）+（test（i+1,j）-eqx（i））*（test（1,j）-eqy）;
end
end
U=0;
fori=1:
M-1
U=U+b（i+1）*Liy（i）;
end
Q1=Lyy-U;
F1=（U/（M-1））/（Q1/（N-M））;
%命令窗口中显示回归方程
disp（''）
disp（'F-分布检验：
'）
stri='F计算值：
';
strin=['方法一：
',num2str（F0）,'方法二：
',num2str（F1）];
string=['自由度：
','f1=',num2str（M-1）,',f2=',num2str（N-M）];
st=cat（2,stri,strin）;
disp（string）
disp（st）
disp（'注：
1.请对照F-分布表找到所需置信水平下的F临界值Fa，若F>Fa，则通过检验。
'）
disp（'2.方法一求回归平方和时用到了因变量的拟合值，方法二用的是偏回归系数。
'）
附（数据较多，显示不便，直接复制到.mat文件中即可）：
48.25
193.72
413.94
358.6
615.04
752.42
435.43
238.55
87.85
316
503.73
554.04
502.07
611.78
603.66
501.67
540.16
264.15
427.11
513.09
478.21
395.25
650.14
480.24
85.79
144.38
39.17
65.32
71.88
58.57
54.33
106.33
257.21
114.53
127.49
194.9
331.09
110.57
194.42
163.89
389.9
541.5
573.03
521.31
645.21
466.28
558.83
621.02
515.02
545.72
786.75
584.89
574
40.5
36.6
35.53
37.48
35.43
33.82
35.63
36.57
39.77
36.05
34.2
35.38
35.62
34
34.38
34.73
34.58
37.2
35.12
35.42
34.73
35.85
33.75
33.4
41.8
41.58
40.15
40.26
40.72
40
40.3
39.37
38.83
39.15
38.93
38.78
38.45
38.63
38.23
37.92
37.2
36.58
35.73
35.73
35.15
35.52
32.95
34.03
34.7
35
34.21
35.43
36.14
1170.8
1707.2
1908.8
2072.4
2136.4
930.8
2025.1
1397.8
1477.2
1517.2
1410
1886.6
1917
3471.4
2314.6
1250
1131.7
2726.7
1765
2450
1495
1873.7
970
1079
1770
2159
1138.7
1526
1591
1270.2
1177.4
1332.2
2311.8
1453.7
1482.7
1764.6
2271
1367
1976.1
1530.8
3045.1
1255.6
1421.9
1346.6
1360
1650
1014.3
1753.2
2810.2
2915.7
3362.7
1221.2
1111.7
K阶自回归拟合及二阶自回归预测的Matlab实现
%k阶自回归及二阶自回归预测
%数据格式：
单行
%数据文件名：
data_zhg.mat,变量名：
test
clear,clc
loaddata_zhg.mat
M=length（test）;
k=input（'请输入自回归阶数：
'）;
y=zeros（k+1,M-k）;
fori=1:
k+1
y（i,:
）=test（i:
i+M-k-1）;
end
y=flipud（y）;
[m,N]=size（y）;
Y=y（1,:
）';
X=cat（2,ones（N,1）,y（2:
end,:
）'）;
A=X'*X;
B=X'*Y;
b=A\B;
formatlongg
str=['首项为常数项，其余系数按距离预测值远近排列为：
',num2str（b'）];
disp（str）
%二元自回归预测
ifk==2
Yuce2=b
（1）+b
（2）*test（M）+b（3）*test（M-1）;
str=['二阶自回归预测下一个值为：
',num2str（Yuce2）];
disp（str）
end
附（直接粘贴到.mat文件中即可）：
3149.44
3303.66
3010.3
3109.61
3639.21
3253.8
3466.5
3839.9
3894.66
4009.61
4253.25
4101.5
4119.88
4258.65
4401.79
一次指数平滑预测的matlab实现
%一次指数平滑及其预测
%数据格式：
单行（按时间序列排序）
%数据文件名：
data_ph.mat,变量名：
test
loaddata_ph.mat
M=length（test）;
disp（'参考：
时间序列稳定、数据波动较小时，a取（0.05,0.3）；否则取（0.7,0.95）'）
a=input（'请输入平滑系数a：
'）;
y=zeros（1,M）;
fori=1:
M
forj=0:
i-1
y（i）=y（i）+test（i-j）*a*（1-a）^j;
end
end
yy=a*test（M）+（1-a）*y（M）;
formatcompact
formatshort
data=test,result=y
format
str=['下一时段数值预测：
',num2str（yy）];
disp（str）
附（直接复制到.mat文件中即可）：
（某城市1999-2004年用水量数据）
211.3
260.18
209.1
248.79
241
250
n次指数平滑及其预测
%数据格式：
单行（按时间序列排序）
%数据文件名：
data_ph.mat,变量名：
test
loaddata_ph.mat
M=length（test）;
disp（'参考：
时间序列稳定、数据波动较小时，a取（0.05,0.3）；否则取（0.7,0.95）'）
a=input（'请输入平滑系数a：
'）;
k=input（'请输入平滑次数：
'）;
y0=test;
y=zeros（1,M）;
fort=1:
k
fori=1:
M
forj=0:
i-1
y（i）=y（i）+y0（i-j）*a*（1-a）^j;
end
end
ift~=k
y0=y;
end
end
yy=a*y0（M）+（1-a）*y（M）;
formatcompact
formatshort
data=test,result=y
format
附（直接复制到.mat文件中即可）：
（某城市1999-2004年用水量数据）
211.3
260.18
209.1
248.79
241
250
一维数据移动平滑的matlab实现
%平滑是为了消除偶然因素对地理变量的影响
%以使地理变量随时间发展变化的趋势和方向明显化
%常用的有三种：
移动平均、滑动平均、指数平滑
%数据格式：
单行（按时间序列排序）
%数据文件名：
data_ph.mat,变量名：
test
loaddata_ph.mat
M=length（test）;
k=input（'请输入移动点数（时距）：
'）;
y=zeros（1,M）;
ifkfori=1:
M-k
forj=1:
k
y（i+k）=y（i+k）+test（i+k-j）;
end
y（i+k）=y（i+k）/k;
end
y（1:
k）=NaN;yy=y（M）+（y（M）-y（M-k））/k;
str=[int2str（k）,'点移动平均结果如下：
'];
disp（str）
formatcompact
formatshort
data=test,result=y
if~isnan（y（M-k））
stri=['预测下一时段数值为：
',num2str（yy）];
disp（stri）
else
disp（'移动点数太大，无法预测下一时段数值！
'）
end
format
else
disp（'Error:
数据个数不足！
'）
end
附（直接复制到.mat文件中即可）：
（某城市1999-2004年用水量数据）
211.3
260.18
209.1
248.79
241
250
K阶自相关系数的matlab实现（含置信度检验）
%k阶自相关系数的计算
%数据格式：
单行
%数据文件名：
data_zxgk.mat，变量名：
test
clear,clc
loaddata_zxgk.mat
M=length（test）;
k=input（'请输入自相关阶数K：
'）;
ifk<=M-1
eqt=sum（test（1:
M-k））/（M-k）;
eqtk=sum（test（1+k:
M））/（M-k）;
A=0;B=0;C=0;
fori=1:
M-k
A=A+（test（i）-eqt）*（test（i+k）-eqtk）;
B=B+（test（i）-eqt）^2;
C=C+（test（i+k）-eqtk）^2;
end
rk=A/sqrt（B*C）;
str=['##',int2str（k）,'阶自相关系数为：
r（',int2str（k）,'）=',num2str（rk）];
disp（str）
stri=['##自由度（f=n-K）：
',int2str（M-k）];
disp（stri）
disp（'##欲检验置信度，请查阅相关系数检验临界值ra.若r>ra则通过检验'）
else
disp（'Error:
K值太大！
'）
end
附（直接复制到.mat文件中即可）：
3149.44
3303.66
3010.3
3109.61
3639.21
3253.8
3466.5
3839.9
3894.66
4009.61
4253.25
4101.5
4119.88
4258.65
4401.79
虽然matlab里面有这些函数，但是要求自己编写，计算机视觉上有这个实验，是别人编写的。
别人到网上找了半天才零散的找到一些碎片，整理以后发上来的！
MatLab自编的均值滤波、中值滤波、高斯滤波图像处理函数。
%自编的均值滤波函数。
x是需要滤波的图像,n是模板大小（即n×n）
functiond=avefilt（x,n）
a（1:
n,1:
n）=1; %a即n×n模板,元素全是1
p=size（x）; %输入图像是p×q的,且p>n,q>n
x1=double（x）;
x2=x1;
%A（a:
b,c:
d）表示A矩阵的第a到b行,第c到d列的所有元素
fori=1:
p
（1）-n+1
forj=1:
p
（2）-n+1
c=x1（i:
i+（n-1）,j:
j+（n-1））.*a; %取出x1中从（i,j）开始的n行n列元素与模板相乘
s=sum（sum（c））; %求c矩阵（即模板）中各元素之和
x2（i+（n-1）/2,j+（n-1）/2）=s/（n*n）; %将模板各元素的均值赋给模板中心位置的元素
end
end
%未被赋值的元素取原值
d=uint8（x2）;
%自编的中值滤波函数。
x是需要滤波的图像,n是模板大小（即n×n）
functiond=midfilt（x,n）
p=size（x）; %输入图像是p×q的,且p>n,q>n
x1=double（x）;
x2=x1;
fori=1:
p
（1）-n+1
forj=1:
p
（2）-n+1
c=x1（i:
i+（n-1）,j:
j+（n-1））;%取出x1中从（i,j）开始的n行n列元素,即模板（n×n的）
e=c（1,:
）; %是c矩阵的第一行
foru=2:
n
e=[e,c（u,:
）]; %将c矩阵变为一个行矩阵
end
mm=median（e）; %mm是中值
x2（i+（n-1）/2,j+（n-1）/2）=mm; %将模板各元素的中值赋给模板中心位置的元素
end
end
%未被赋值的元素取原值
d=uint8（x2）;
%自编的高斯滤波函数，S是需要滤波的图象，n是均值，k是方差
functiond=gaussfilt（k,n,s）
Img=double（s）;
n1=floor（（n+1）/2）;%计算图象中心
fori=1:
n
forj=1:
n
b（i,j）=exp（-（（i-n1）^2+（j-n1）^2）/（4*k））/（4*pi*k）;
end
end
%生成高斯序列b。
Img1=conv2（Img,b,'same'）;%用生成的高斯序列卷积运算，进行高斯滤波
d=uint8（Img1）;
%此为程序主文件，包含主要功能单元，以及对子函数进行调用
try
%实验步骤一：
彩色、灰度变换
h=imread（'photo.jpg'）; %读入彩色图片
c=rgb2gray（h）; %把彩色图片转化成灰度图片，256级
figure,imshow（c）,title（'原始图象'）;%显示原始图象
g=imnoise（c,'gaussian',0.1,0.002）;%加入高斯噪声
figure,imshow（g）,title（'加入高斯噪声之后的图象'）; %显示加入高斯噪声之后的图象
%实验步骤二：
用系统预定义滤波器进行均值滤波
n=input（'请输入均值滤波器模板大小\n'）;
A=fspecial（'average',n）;%生成系统预定义的3X3滤波器
Y=filter2（A,g）/255; %用生成的滤波器进行滤波,并归一化
figure,imshow（Y）,title（'用系统函数进行均值滤波后的结果'）;%显示滤波后的图象
%实验步骤三:
用自己的编写的函数进行均值滤波
Y2=avefilt（g,n）; %调用自编函数进行均值滤波，n为模板大小
figure,imshow（Y2）,title（'用自己的编写的函数进行均值滤波之后的结果'）;%显示滤波后的图象
%实验步骤四:
用Matlab系统函数进行中值滤波
n2=input（'请输入中值滤波的模板的大小\n'）;
Y3=medfilt2（g,
升级会员 