matlab-數(shù)據(jù)分析和統(tǒng)計

數(shù)據(jù)分析和統(tǒng)計

面向列的數(shù)據(jù)集

這年頭似乎十分風(fēng)行"面向”這個詞，這兒故也套用，其英文為“Column-OrientedDataSets”,可理

解為MatLab按列的存儲方式來分析數(shù)據(jù)，下面是一個例子：

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

以上數(shù)據(jù)被保存在一個稱為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

下面，我們調(diào)入此文件，并看看文件的一些參數(shù)

loadcount.dat

[n,p]=size(count)

n=

24

P=

3

創(chuàng)建一個時間軸后，我們可以把圖畫出來：

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次觀測.

基本數(shù)據(jù)分析函數(shù)

(一定注意是面向列的)

繼續(xù)用上面的數(shù)據(jù)，其每列最大值.均值.及偏差分別為:

mx=max(count)

mu=mean(count)

sigma=std(count)

mx=

114145257

mu=

32.000046.541765.5833

sigma=

25.370341.405768.0281

重載函數(shù)，還可以定位出最大.最小值的位置

[mxjndx]=min(count)

mx=

797

indx=

22324

試試看，你能看懂下面的命令是干什么的嗎？

[n,p]=size(count)

e=ones(n,1)

x=count-e*mu

點這看看答案!

下面這句命令則找出了整個矩陣的最小值：

min(count(:))

ans=

7

協(xié)方差及相關(guān)系數(shù)

下面,我們來看看第一列的方差：

cov(count(:,1))

ans=

643.6522

cov()函數(shù)作用于矩陣,則會計算其協(xié)方差矩陣.

corrcoef()用于計算相關(guān)系數(shù)，如:

corrcoef(count)

ans=

1.00000.93310.9599

0.93311.00000.9553

0.95990.95531.0000

數(shù)據(jù)的預(yù)處理

未知數(shù)據(jù)

NaN(NotaNumber—不是一個數(shù))被定義為未經(jīng)定義的算式的結(jié)果，如0/0.在處理數(shù)據(jù)中,NaN常用來表示

未知數(shù)據(jù)或未能獲得的數(shù)據(jù).所有與NaN有關(guān)的運算其結(jié)果都是NaN.

a=magic(3);

a(2,2)=NaN

a=

816

3NaN7

492

sum(a)

ans=

15NaN15

在做統(tǒng)計時，常需要將NaN轉(zhuǎn)化為可計算的數(shù)字或去掉，以下是幾種方法：

注:判斷一個值是否為NaN,只能用isnan(),而不可用x—NaN;

i=find(~isnan(x));

先找出值不是NaN的項的下標(biāo),將這些元素保留

x=x(i)

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

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

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

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

用此法可以從數(shù)據(jù)中去掉不相關(guān)的數(shù)據(jù)，看看下面的命令是干什么用的:

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'),:)=[];

點這看看答案

回歸與曲線擬合

我們經(jīng)常需要把觀測到的數(shù)據(jù)表達(dá)為函數(shù)，假如有如下的對時間的觀測：

t=[0.3.81.11.62.3],;

y=[0.50.821.141.251.351.40],;

plot(t,y「o)

gridon

多項式回歸

由圖可以看出應(yīng)該可以用多項式來表達(dá)：y=a0+al*t+a2*t"2

系數(shù)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即為待求的系數(shù)，畫圖比較可得

T=(0:0.1:2.5),;

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

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

結(jié)果令人失望，但我們可以增加階數(shù)來提高精確度，但更明智的選擇是用別的方法.

線性參數(shù)回歸

形如：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

看起來是不是好多了！

例子研究：曲線擬合

下面我們以美國人口普杳的數(shù)據(jù)來研究一下有關(guān)曲線擬合的問題(MatLab是別人的,教學(xué)文檔是別人

的，例子也是別人的，我只是一個翻譯而已...)

loadcensus

這樣我們得到了兩個變量，cdate是1790至1990年的時間列向量(10年一次)，pop是相應(yīng)人口數(shù)列向量.

上一小節(jié)所講的多項式擬合可以用函數(shù)polyfit()來完成,數(shù)字指明了階數(shù)

p=polyfit(cdate,pop,4)

Warning:Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=5.429790e-20

P=

1.0e+05*

0.0000-0.00000.0000-0.01266.0020

產(chǎn)生警告的原因是計算中的cdata值太大，在計算中的Vandermonde行列式使變換產(chǎn)生了問題，解決的方

法之一是使數(shù)據(jù)標(biāo)準(zhǔn)化.

預(yù)處理:標(biāo)準(zhǔn)化數(shù)據(jù)

數(shù)據(jù)的標(biāo)準(zhǔn)化是對數(shù)據(jù)進(jìn)行縮放，以使以后的計算能更加精確，一種方法是使之成為0均值：

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

現(xiàn)在再進(jìn)行曲線擬合就沒事了！

p=polyfit(sdate,pop,4)

P=

0.70470.921023.470673.859862.2285

pop4=polyval(p,sdate);

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

在上面的數(shù)據(jù)標(biāo)準(zhǔn)化中，也可以有別的算法，如令1790年的人口數(shù)為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,'+')

可以看出，多項式擬合即使提高了階次也無法達(dá)到令人滿意的結(jié)果

指數(shù)擬合

從人口增長圖可以發(fā)現(xiàn)人數(shù)的增長基本是呈指數(shù)增加的，因此我們可以用年份的對數(shù)來進(jìn)行擬合，這

兒，年數(shù)是標(biāo)準(zhǔn)化后的！

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

上面的圖不令人滿意，下面，我們用二階的對數(shù)分析:

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

這種余量分析比多項式擬合的余量分析圖案要隨機(jī)的多(沒有很強(qiáng)的規(guī)律性)，可以預(yù)見，隨著人數(shù)的增

加，余糧所反映的不確定度也在增加，但總的說來，這種擬合方式要強(qiáng)好多！

誤差邊界

誤差邊界常用來反映你所用的擬合方式是否適用于數(shù)據(jù)，為得到誤差邊界，只需在polyfitO中傳遞第二

個參數(shù)，并將其送入polyval().

下面是一個二階多項式擬合模型,年份已被標(biāo)準(zhǔn)化,下面的代碼用了2。，對應(yīng)于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中的差分和濾波基本都是對向量而言的，向量則是存儲取樣信號或序列的.

函數(shù)y=filter(b,a,x)將用a,b描述的濾波器處理向量x,然后將其存儲在向量y中,

filter()函數(shù)可看為是一差分方程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];

我們載入數(shù)據(jù),取其第一列，并計算有：

loadcount.dat

x=count(:,1);

y=filter(b,a,x);

t=1:length(x);

gridon

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

實現(xiàn)所表示的就是濾波后的數(shù)據(jù),它代表了4小時的平均車流量

MatLab的信號處理工具箱中提供了很多用來濾波的函數(shù),可用來處理實際問題!

快速傅立葉變換(FFT)

傅立葉變換能把信號按正弦展開成不同的頻率值，對于取樣信號，用的是離散傅立葉變換.

FET是離散傅立葉變換的一種高速算法,在信號和圖像處理中有極大的用處！

fft離散傅立葉變換

fft2二維離散傅立葉變換

fftnn維離散傅立葉變換

ifft離散傅立葉反變換

ifft2二維離散傅立葉反變換

ifftnn維離散傅立葉反變換

abs幅度

angle相角

unwrap相位按弧度展開,大于n的變換為2n的補(bǔ)角

fftshift把零隊列移至功率譜中央

cplxpair把數(shù)據(jù)排成復(fù)數(shù)對

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雖然是實數(shù),但y是復(fù)數(shù),其中,第一個是因為它是常數(shù)相加的結(jié)果,第五個則對應(yīng)于奈奎斯特頻率,后三個

數(shù)是由于負(fù)頻率的影響,它們是前面三個數(shù)的共能值！

下面，讓我們來驗證一下太陽黑子活動周期是11年!Wolfer數(shù)記錄了300年太陽黑子的數(shù)量及大?。?/p>

loadsunspot.dat

year=sunspot(:,1);

woIfer=sunspot(:,2);

plot(year,wolfer)

title('SunspotData')

SunspotData

200

現(xiàn)在來看看其FFT:

Y=fft(wolfer);

Y的幅度是功率譜,畫出功率譜和頻率的對應(yīng)關(guān)系就得出了周期圖,去掉第一點,因為他只是所有數(shù)據(jù)的和,

畫圖有：

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。是用來計算幅度和相位的

先創(chuàng)建一信號,再進(jìn)行分析,unwarp()把相位大于n的變換為2”的補(bǔ)角:

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