數據分析和統計

面向列的數據集

這年頭似乎十分風行"面向”這個詞，這兒故也套用，其英文為“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