MATLAB应用——数据分析与统计

7 \* _$ k2 H3 M7 g8 C0 O d

数学建模是用数学方法解决各种实际问题的桥梁，它已经渗透到各个领域，而且发挥出越来越重要的作用。面对自然科学和工程应用中的难题，大部分人无从入手，而个别人却能短时间内给出切实可行的解决方案，其差别往往在于驾驭数学知识的能力不同。现代计算机技术的应用不仅减少了计算错误，而且加强了数学应用者解决问题的能力。MATLAB是一款常用的数据处理软件，为了更好的应用MATLAB软件，我将整理好的MATLAB函数分享到今日头条上，以利己利人查阅。

MATLAB提供的很多数据分析与统计函数都是面向列的，即矩阵中的每一列代表一个变量的多个观测值，其列数对应于变量数，行数对应于测量点数。

max和min函数可求出数据的最大值和最小值，mean和std函数可求出数据的均值和标准差，sum和prod函数可求出数据元素和与数据元素积。例如，对MATLAB内含的某城市24小时的车流量数据count.dat可作分析：

load count.dat

* `$ v- N- k. {' N! z( y2 F- N

mx=max(count)

mx = 114 145 257

- F- F* W" N; n% G/ Q* A

mu=mean(count)

6 g: ~. T- D: n6 w7 b: M$ P7 @! `: n" ^

mu = 32.0000 46.5417 65.5833

sigma=std(count)

sigma = 25.3703 41.4057 68.0281

对有些函数还可给出位置，例如，在求出最小值的同时，可得到最小值所在的位置(行号)：

* N' H' W9 z+ e6 r" \2 _) t7 W& J h

[mx,indx]=min(count)

+ V: ~# p# E, i) x: }% F

mx = 7 9 7

. u% L, }- y2 f9 x* W+ y0 ?, a m

indx = 2 23 24

1、协方差和相关系数

8 H6 d4 u1 @. s* b/ F! x

cov函数可以求出单个变量的协方差，而corrcoef函数可求出两个变量之间的相关系数，例如：

" m# F( z8 T a

cv=cov(count)

" a0 J& V! Z4 P& S# ]5 C' i

cv = 1.0e+003 *

0.6437 0.9802 1.6567

9 t2 }& p$ R5 W7 I

0.9802 1.7144 2.6908

u$ _0 ]( A7 a( V) i

1.6567 2.6908 4.6278

: Q, t0 L$ @8 e+ I

cr=corrcoef(count)

7 u8 \$ Z/ \, e- Q6 k% J

cr =

1.0000 0.9331 0.9599

0.9331 1.0000 0.9553

. @4 \$ t: u4 f# Q. d' c1 B

0.9599 0.9553 1.0000

~9 J+ |6 A# n1 M/ M

2、数据预处理

在MATLAB中遇到超出范围的数据时均用NaN (非数值) 表示，而且在任何运算中，只要包含NaN，就将它传递到结果中，因此在对数据进行分析前，应对数据中出现的NaN作剔除处理。例如：

4 ]6 L, H1 W! p9 E; q' u# H7 @

a=[1 2 3;5 NaN 8;7 4 2];

1 A9 J+ _6 |0 c/ e5 G' N+ s2 t6 C

sum(a)

ans = 13 NaN 13

* z2 L+ d+ D2 b8 p

在矢量x中删除NaN元素，可有下列四种方法：

(1)  i=find(~isnan(x))；x=x(i)。

8 V& z2 y# ?' }! ~; Y

(2)  x=x(find(~isnan(x)))。

(3)  x=x(~isnan(x))。

# `1 b! I3 y, Y2 @0 C# b

(4)  x(isnan(x))=[ ]。

9 n5 q2 w% E& t" ]8 l9 s6 z

在矩阵X中删除NaN所在的行，可输入

# v& V( M' ~ \0 [

X(any(isnan(X))，：)=[ ]；

+ W) h' X' {( P, \1 J7 m

经过这种预处理后的数据，可进行各种分析和统计操作。

: f1 e4 I2 H$ f7 u" D

3、回归和曲线拟合

+ ^# d: p* K2 G9 f0 T

对给定的数据进行拟合，可采用多项式回归，也可采用其它信号形式的回归，其基本原理是最小二乘法，这一功能实现在MATLAB中显得轻而易举。

例1：设通过测量得到一组时间t与变量y的数据：

t=[0 .3 .8 1.1 1.6 2.3];

' G# f$ P! j! l. m

y=[0.5 0.82 1.14 1.25 1.35 1.40];

L& ]- Y K) q0 ~- ^ m: V

8 Q: C8 I' q+ Y5 s3 X

进行回归，可得到两种不同的结果。MATLAB程序如下：

" G# o! X( W1 l7 }

t=[0 .3 .8 1.1 1.6 2.3];

" q/ r7 `" o5 K4 [! [& n) M

y=[.5 .82 1.14 1.25 1.35 1.40];

8 {2 E3 L U( E; \ K! w

X1=[ones(size(t)) t t.^2];

% x# c2 {: x. G j/ b& d

a=X1\y;

8 ?1 W8 B }1 a1 y) X& D& X

X2=[ones(size(t)) exp(–t) t.*exp(–t)];

7 E: X# p [, A) d1 K0 J: C

b=X2\y;

T=[0:.1:2.5];

- r# g' J- D% r& Z, {% ^7 D

Y1=[ones(size(T)) T T.^2]*a;

Y2=[ones(size(T)) exp(-T) T.*exp(-T)]*b;

figure(1)

subplot(1,2,1)

plot(T,Y1,-,t,y,o),grid on

) z O' Q! `4 Q- V, `7 N

title(多项式回归)

2 m6 E/ ~) q- A9 T/ m6 u3 L9 q

subplot(1,2,2)

% ?/ {. [- R, J8 V4 V! {

plot(T,Y2,-,t,y,o),grid on

6 P$ o2 X2 t) q e

title(指数函数回归)

n0 k7 ]. a0 R2 }( t# S/ n4 ^' q

; h, i- a0 G* I1 E' h& }

例2 已知变量y与x1，x2有关，测得一组数据为

, ]. j$ Z# N- F" g2 ]

　　x1=[.2 .5 .6 .8 1.0 1.1 ];

' ^0 Z9 r; \) g' u: D0 I

　　x2=[.1 .3 .4 .9 1.1 1.4 ];

+ K7 L5 D6 I& k

　　y=[.17 .26 .28 .23 .27 .24];

采用来拟合，则有

% e; i% P9 H. V# T. b& O2 g% F: m

x1=[.2 .5 .6 .8 1.0 1.1];

. e2 v! e" G/ _, o; p

x2=[.1 .3 .4 .9 1.1 1.4];

y=[.17 .26 .28 .23 .27 .24];

X=[ones(size(x1)) x1 x2];

8 H6 E- a( u+ [

a=X\y

a = 0.1018 0.4844 −0.2847

3 o' m3 W- L9 @

因此数据的拟合模型为

5 k& P" [( _. a8 U' b N' p

y=0.1018+0.4844x1−0.2487x2

4、傅里叶分析与FFT

利用MATLAB提供的FFT函数可方便地计算出信号的傅里叶变换，从而在频域上对信号进行分析。

例1 ：混合频率信号成分分析。有一信号x由三种不同频率的正弦信号混合而成，通过得到信号的DFT，确定出信号的频率及其强度关系，程序如下：

t=0:1/119:1;

x=5*sin(2*pi*20*t)+3*sin(2*pi*30*t)+sin(2*pi*45*t);

& B6 l2 y b) u: u. L7 T1 F. ^

y=fft(x);

m=abs(y);

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

figure(1)

subplot(2,1,1),plot(t,x),grid on

title(多频率混合信号)

7 V6 X, F% y( M% }' S

ylabel(Input \itx),xlabel(Time )

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

ylabel(Abs. Magnitude),grid on

xlabel(Frequency (Hertz))

3 U4 f# I9 F z" t; r' r+ h1 k

0 h+ J7 C6 c# ]/ u. @* k) t

例2 ：信号在传输过程中，由于受信道或环境影响，在接收端得到的是噪声环境下的信号。我们利用FFT函数对这一信号进行傅里叶分析，从而确定信号的频率，程序如下：

t=0:1/199:1;

x=sin(2*pi*50*t)+1.2*randn(size(t)); ％噪声中的信号

y=fft(x);

0 P! q2 Y% k( ^- W3 ]& s3 U* }

m=abs(y);

. o* L9 q2 \8 m( Y9 O' A# u

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

' z. f4 {' L# G5 L+ }) R

figure(1)

' p$ c' B( O( r$ A5 E+ G4 ^

subplot(2,1,1),plot(t,x),grid on

title(信号检测)

ylabel(Input \itx),xlabel(Time )

+ _- w9 G+ }$ L

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

ylabel(Abs. Magnitude),grid on

xlabel(Frequency (Hertz))

5 q9 ] w0 q9 Z# i

! `# s* U( u: V

例3 ：天文学家记录了300年来太阳黑子的活动情况，我们对这组数据进行傅里叶分析，从而得出太阳黑子的活动周期。MATLAB程序如下：

# m# C- g# R, F

load sunspot.dat

year=sunspot(:,1);

- n0 Q, z' g# C2 }

wolfer=sunspot(:,2);

# B# t: Y# r6 O% z0 ^7 @3 N

figure(1)

Z6 z! [1 A' D: n% _

subplot(2,1,1)

8 C+ t6 J$ W- f6 m+ |- L; F) h

plot(year,wolfer)

- L3 h" V7 r$ q5 B

title(原始数据)

Y=fft(wolfer);

N=length(Y);

Y(1)=[];

8 I" V1 l% L1 n) n. K9 i. Z3 c. x

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

9 G* b% @7 g/ v: d

nyquist=1/2;

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

period=1./freq;

/ |" l" S* c( z! k+ P8 h6 n

subplot(2,1,2)

- y- ]& l& j, w' H# y7 f) j

plot(period,power)

title(功率谱), grid on

axis([0 40 0 2e7])

. F& {% K5 b% k' ]6 G0 e" {/ T

各位读者朋友，感谢您的阅读，您若对工程应用中的数学问题感兴趣，欢迎关注我，愿我们一起讨论和成长！！！

' f1 n7 ?$ |# \8 } ' v: r; Y) m2 I1 A6 F& v! n* z! g / c2 f. q8 \4 o0 t& ~% B/ v# ? . F% e6 b0 l# m! K