matlab-数据分析和统计

数据分析和统计

面向列的数据集

这年头似乎十分风行"面向”这个词，这儿故也套用，其英文为“Column-OrientedDataSets”,可理

解为MatLab按列的存储方式来分析数据，下面是一个例子：

TimeLocation1Location2Location3

OlhOO11119

02h0071311

03h00141720

04h0011139

05h00435169

06h00384676

07h0061132186

08h0075135180

09h003888115

lOhOO283655

HhOO121214

12h00182730

13h00181929

14h00171518

15h00193648

16h00324710

17h00426592

18h005766151

19h00445590

20h00114145257

21h00355868

22h00111215

23h0013915

24h001097

以上数据被保存在一个称为count,dat的文件中.

11119

71311

141720

11139

435169

384676

61132186

75135180

3888115

283655

121214

182730

181929

171518

193648

324710

426592

5766151

445590

114145257

355868

111215

13915

1097

下面，我们调入此文件，并看看文件的一些参数

loadcount.dat

[n,p]=size(count)

n=

24

P=

3

创建一个时间轴后，我们可以把图画出来：

t=1:n;

set(O,'defau代axeslinestyleorder','■卜・卜.')

se^O/defaultaxescolororde^JO00])

plot(t,count),legend('Location1VLocation2'；Location3',0)

xlabel('Time'),ylabel('VehicleCount'),gridon

300

—Location1

—Location2

250•一--•-Location3

200

-

S

o

o

9

P

s

>

100

50

足以证明，以上是对3个对象的24次观测.

基本数据分析函数

(一定注意是面向列的)

继续用上面的数据，其每列最大值.均值.及偏差分别为:

mx=max(count)

mu=mean(count)

sigma=std(count)

mx=

114145257

mu=

32.000046.541765.5833

sigma=

25.370341.405768.0281

重载函数，还可以定位出最大.最小值的位置

[mxjndx]=min(count)

mx=

797

indx=

22324

试试看，你能看懂下面的命令是干什么的吗？

[n,p]=size(count)

e=ones(n,1)

x=count-e*mu

点这看看答案!

下面这句命令则找出了整个矩阵的最小值：

min(count(:))

ans=

7

协方差及相关系数

下面,我们来看看第一列的方差：

cov(count(:,1))

ans=

643.6522

cov()函数作用于矩阵,则会计算其协方差矩阵.

corrcoef()用于计算相关系数，如:

corrcoef(count)

ans=

1.00000.93310.9599

0.93311.00000.9553

0.95990.95531.0000

数据的预处理

未知数据

NaN(NotaNumber—不是一个数)被定义为未经定义的算式的结果，如0/0.在处理数据中,NaN常用来表示

未知数据或未能获得的数据.所有与NaN有关的运算其结果都是NaN.

a=magic(3);

a(2,2)=NaN

a=

816

3NaN7

492

sum(a)

ans=

15NaN15

在做统计时，常需要将NaN转化为可计算的数字或去掉，以下是几种方法：

注:判断一个值是否为NaN,只能用isnan(),而不可用x—NaN;

i=find(~isnan(x));

先找出值不是NaN的项的下标,将这些元素保留

x=x(i)

x=x(find(~isnan(x)))同上,去掉NaN

x二x(、isnan(x));更快的做法

x(isnan(x))=口；消掉NaN

X(any(isnan(X)'),:)=口;把含有NaN的行都去掉

用此法可以从数据中去掉不相关的数据，看看下面的命令是干什么用的:

mu=mean(count);

sigma=std(count);

[n,p]=size(count)

outliers=abs(count—mu(ones(n,1),:))>3*sigma(ones(n,1),:);

nout=sum(outliers)

nout=

100

count(any(outliers'),:)=[];

点这看看答案

回归与曲线拟合

我们经常需要把观测到的数据表达为函数，假如有如下的对时间的观测：

t=[0.3.81.11.62.3],;

y=[0.50.821.141.251.351.40],;

plot(t,y「o)

gridon

多项式回归

由图可以看出应该可以用多项式来表达：y=a0+al*t+a2*t"2

系数aO,al,a2可以由最小平方拟合来确定，这一步可由反除号”\”来完成

解下面的三元方程组可得：

X=[ones(size(t))tt.A2]

X=

1.000000

1.00000.30000.0900

1.00000.80000.6400

1.00001.10001.2100

1.00001.60002.5600

1.00002.30005.2900

a=X\y

a=

0.53180.9191-0.2387

a即为待求的系数，画图比较可得

T=(0:0.1:2.5),;

Y=[ones(size(T))TT.A2]*a;

plot(T,Y,Tt,y,'o[),gridon

结果令人失望，但我们可以增加阶数来提高精确度，但更明智的选择是用别的方法.

线性参数回归

形如：y=aO+al*exp(-t)+a2*t*exp(­t)

计算方法同上：

X=[ones(size(t))exp(-1)t.*exp(-1)];

a=X\y

a=

1.3974-0.89880.4097

T=(0:0.1:2.5),;

Y=[ones(size(T))exp(-T)T.exp(-T)]*a;

plot(T,Y,'-',t,y,'o'),gridon

1.5

看起来是不是好多了！

例子研究：曲线拟合

下面我们以美国人口普杳的数据来研究一下有关曲线拟合的问题(MatLab是别人的,教学文档是别人

的，例子也是别人的，我只是一个翻译而已...)

loadcensus

这样我们得到了两个变量，cdate是1790至1990年的时间列向量(10年一次)，pop是相应人口数列向量.

上一小节所讲的多项式拟合可以用函数polyfit()来完成,数字指明了阶数

p=polyfit(cdate,pop,4)

Warning:Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=5.429790e-20

P=

1.0e+05*

0.0000-0.00000.0000-0.01266.0020

产生警告的原因是计算中的cdata值太大，在计算中的Vandermonde行列式使变换产生了问题，解决的方

法之一是使数据标准化.

预处理:标准化数据

数据的标准化是对数据进行缩放，以使以后的计算能更加精确，一种方法是使之成为0均值：

sdate=(cdate-mean(cdate))./std(cdate)

现在再进行曲线拟合就没事了！

p=polyfit(sdate,pop,4)

P=

0.70470.921023.470673.859862.2285

pop4=polyval(p,sdate);

plot(cdate,pop4,'-'.cdate,pop,gridon

在上面的数据标准化中，也可以有别的算法，如令1790年的人口数为0.

余量分析

p1=polyfit(sdate,pop,1);

pop1=polyval(p1,sdate);

res1=pop-pop1;

figure,plot(cdate,res1

p=polyfit(sdate,pop,2);

pop2=polyval(p,sdate);

res2=pop-pop2;

figure,plot(cdate,res2,,+,)

p=polyfit(sdate,pop,4)；

pop4=polyval(p,sdate);

plot(cdate,pop4,'-,,cdate,pop,'+')

res4=pop-pop4;

figure,plot(cdate,res4,'+')

可以看出，多项式拟合即使提高了阶次也无法达到令人满意的结果

指数拟合

从人口增长图可以发现人数的增长基本是呈指数增加的，因此我们可以用年份的对数来进行拟合，这

儿，年数是标准化后的！

logpl=polyfit(sdate,log10(pop),1);

logpredl=10.Apolyval(logp1,sdate);

semilogy(cdate,logpred1,'-',cdate,pop,'+');

gridon

Iogres2=Iog10(pop)-polyval(logp2,sdate);

plot(cdate,logres2,'+')

0.03

0.02

0.01

0

-0.01

-0.02

-0.03

-0.04

175018001850190019502000

上面的图不令人满意，下面，我们用二阶的对数分析:

Iogp2=polyfit(sdate,log10(pop),2);

Iogpred2=10.Apolyval(logp2,sdate);

,

semilogy(cdateJogpred25'-,cdateJpopJ'4-');

gridon

r=pop-10.A(polyval(logp2,sdate));

plot(cdate,r,'+')

10

+

5

十

+

+T-

+

+

+

0+

++++'十

+

-5

-10

十：

-15」

175018001850190019502000

这种余量分析比多项式拟合的余量分析图案要随机的多(没有很强的规律性)，可以预见，随着人数的增

加，余粮所反映的不确定度也在增加，但总的说来，这种拟合方式要强好多！

误差边界

误差边界常用来反映你所用的拟合方式是否适用于数据，为得到误差边界，只需在polyfitO中传递第二

个参数，并将其送入polyval().

下面是一个二阶多项式拟合模型,年份已被标准化,下面的代码用了2。，对应于95%的可置信度：

[p2,S2]=polyfit(sdate,pop,2);

[pop2,del2]=polyval(p2,sdate,S2);

plot(cdate,pop,'+',cdate,pop2,'g-',cdate,pop2+2*del2,'匚

cdate,pop2-2*del2,r)

gridon

300

差分方程和滤波

MatLab中的差分和滤波基本都是对向量而言的，向量则是存储取样信号或序列的.

函数y=filter(b,a,x)将用a,b描述的滤波器处理向量x,然后将其存储在向量y中,

filter()函数可看为是一差分方程aly(n)=bl*x(1)+b2*x(2)+...-a2*y(2)-…

如有以下差分方程：y(n)=l/4*x(n)+l/4*x(n-l)+l/4*x(n-2)+l/4*x(n-3),则

a=1;

b=[1/41/41/41/4];

我们载入数据,取其第一列，并计算有：

loadcount.dat

x=count(:,1);

y=filter(b,a,x);

t=1:length(x);

gridon

legend('OriginalData','SmoothedData',2)

实现所表示的就是滤波后的数据,它代表了4小时的平均车流量

MatLab的信号处理工具箱中提供了很多用来滤波的函数,可用来处理实际问题!

快速傅立叶变换(FFT)

傅立叶变换能把信号按正弦展开成不同的频率值，对于取样信号，用的是离散傅立叶变换.

FET是离散傅立叶变换的一种高速算法,在信号和图像处理中有极大的用处！

fft离散傅立叶变换

fft2二维离散傅立叶变换

fftnn维离散傅立叶变换

ifft离散傅立叶反变换

ifft2二维离散傅立叶反变换

ifftnn维离散傅立叶反变换

abs幅度

angle相角

unwrap相位按弧度展开,大于n的变换为2n的补角

fftshift把零队列移至功率谱中央

cplxpair把数据排成复数对

nextpow2下两个更高的功率

向量x的FFT可以这样求：

x=[437-91000],

y=fft(x)

y=

6.0000

11.4853-2.7574i

-2.0000-12.0000i

-5.4853+11.2426i

18.0000

-5.4853-11.2426i

-2.0000+12.0000i

11.4853+2.7574i

x虽然是实数,但y是复数,其中,第一个是因为它是常数相加的结果,第五个则对应于奈奎斯特频率,后三个

数是由于负频率的影响,它们是前面三个数的共能值！

下面，让我们来验证一下太阳黑子活动周期是11年!Wolfer数记录了300年太阳黑子的数量及大小：

loadsunspot.dat

year=sunspot(:,1);

woIfer=sunspot(:,2);

plot(year,wolfer)

title('SunspotData')

SunspotData

200

现在来看看其FFT:

Y=fft(wolfer);

Y的幅度是功率谱,画出功率谱和频率的对应关系就得出了周期图,去掉第一点,因为他只是所有数据的和,

画图有：

N=length(Y);

丫⑴=□；

power=abs(Y(1:N/2)).A2;

nyquist=1/2;

freq=(1:N/2)/(N/2)*nyquist;

plot(freq,power),

gridon

xlabelCcycles/year')

title(,Periodogram,)

上面的图看起来不大方便,下面我们画出频谱-周期图

period=1./freq;

plot(period,power),

axis([04002e7]),

gridon

ylabelCPower')

xlabel(^Period(Years/Cycle),)

为了得出精确一点的解,如下：

[mpindex]=max(power);

period(index)

ans=

11.0769

变换后的幅度和相位

abs()和angle。是用来计算幅度和相位的

先创建一信号,再进行分析,unwarp()把相位大于n的变换为2”的补角:

t=0:1/99:1;

x=sin(2*pi*15*t)+sin(2*pi*40*t);

y=fft(x);

m=abs(y);

p=unwrap(angle(y));

f=(0:length(y)-1)'*99/length(y);

subplot(2,1,1),plot(f,m),

ylabel('Abs.Magnitude'),gridon

subplot(2,1,2),plot(f,p*180/pi)

ylabel('Phase[Degrees]'),gridon

xlabel('Frequency[Hertz]')

s【

①

①

6」

名

】

①

s

e

e

d

-1000

-12001