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

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

* G' ?9 H7 z5 e' B% O

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

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

load count.dat

, ~9 W1 g$ C# x' |1 _) A: S6 U8 `

mx=max(count)

mx = 114 145 257

mu=mean(count)

9 q, C: L1 c2 h0 y3 f8 s

mu = 32.0000 46.5417 65.5833

sigma=std(count)

4 R9 }/ f( a! ~, T

sigma = 25.3703 41.4057 68.0281

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

8 T1 e; m( p' \. J- U- l5 L7 b

[mx,indx]=min(count)

mx = 7 9 7

indx = 2 23 24

( k, \) f, b& h8 K

1、协方差和相关系数

8 K% b. r8 z% g$ V* w; t* E: I

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

6 C% P* c3 ]: K d# n' l

cv=cov(count)

! n* F' q& Z5 O

cv = 1.0e+003 *

0.6437 0.9802 1.6567

0.9802 1.7144 2.6908

1.6567 2.6908 4.6278

9 v8 M% j, T" E5 K/ E& g4 D

cr=corrcoef(count)

cr =

1.0000 0.9331 0.9599

: C5 _4 M7 X2 f# T+ T% ^

0.9331 1.0000 0.9553

6 V9 T M' v& y* I; ^9 B5 F' _' E

0.9599 0.9553 1.0000

2、数据预处理

3 E( k" a/ y) o: _3 \ E

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

- i* b% C7 @3 |8 _& S

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

sum(a)

" c/ y3 x3 H1 B

ans = 13 NaN 13

# j/ N) }! b8 t. |6 q7 {0 z

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

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

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

3 E; t9 z2 r+ ^

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

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

; W0 Z2 n2 m E4 a$ Z

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

# r7 B- T9 k: y

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

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

8 h4 N# L' Q7 Y' K

3、回归和曲线拟合

6 h2 m- k( q' X6 j

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

7 @) b9 b! D U; C+ a- f

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

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

: \0 q' L! s& r6 o

y=[0.5 0.82 1.14 1.25 1.35 1.40];

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

3 t/ U5 P4 r3 _+ c

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

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

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

2 G. s4 {& |0 \3 C$ i5 Y' g& e: f

a=X1\y;

" u) ]% x) e8 @& M

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

* A0 h3 J0 U8 G6 q) i

b=X2\y;

9 k+ q( W/ S H) N

T=[0:.1:2.5];

- x9 f7 t! N" R! W* b% p

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

- Q# W: @: j! t% a; l! c1 X7 p( R& g

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

figure(1)

1 _4 _1 [& f9 e3 x) d9 m

subplot(1,2,1)

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

n" H$ h/ P4 X% p% j% p4 {' q g

title(多项式回归)

subplot(1,2,2)

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

* i( s7 l2 L' l+ K7 w3 J

title(指数函数回归)

9 _$ G6 [+ _$ w7 z! D1 d

8 x) |5 O! B* o4 Y& L4 ?( I! Y3 m

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

6 O/ B! j# R) r

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

" D/ T3 x% Z" Z* u# a2 Z, i

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

% ~( l; [/ C7 h, F/ ?6 A

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

采用来拟合，则有

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

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

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

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

a=X\y

/ I; W# x3 p* P

a = 0.1018 0.4844 −0.2847

/ w, c% ]3 x" \- a! z

因此数据的拟合模型为

. A! E9 Y7 g6 \8 f* @6 N

y=0.1018+0.4844x1−0.2487x2

4、傅里叶分析与FFT

: C$ b E! c9 M- @

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

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

t=0:1/119:1;

P7 Q: g* w! K; M$ e6 z/ |; b

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

% Y; m! Q3 V; X5 Q" Y2 d5 r" y9 q

y=fft(x);

- C1 E: @/ c7 l7 O- ?! K

m=abs(y);

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

* U/ n1 H5 ~/ N! t- }

figure(1)

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

% w8 ^7 O7 ~2 g2 q3 B3 J. }& p( V2 _

title(多频率混合信号)

; ?/ L s0 Q6 ]. r% [9 R ]. g

ylabel(Input \itx),xlabel(Time )

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

ylabel(Abs. Magnitude),grid on

' B+ j1 b- J4 ?4 {* w+ n

xlabel(Frequency (Hertz))

7 C3 r* H S# `! C, J$ ~+ r; B, D

; Y7 ^; u* N# h: s! X

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

t=0:1/199:1;

* \2 u7 o9 o& C+ j# m/ c3 y8 D

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

7 F0 t* w) A G& E: o9 d" X0 J

y=fft(x);

m=abs(y);

% N5 e- s' m5 X+ y/ F3 _% w _

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

9 A2 r7 K. B) z9 V; M% \

figure(1)

+ Z3 q: m. S5 [9 A

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

' u6 `' ]* X# [

title(信号检测)

ylabel(Input \itx),xlabel(Time )

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

ylabel(Abs. Magnitude),grid on

xlabel(Frequency (Hertz))

& `% L3 `: l; n

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

0 u- j0 h' Z* H

load sunspot.dat

$ z4 N) p' M4 G% @

year=sunspot(:,1);

wolfer=sunspot(:,2);

figure(1)

subplot(2,1,1)

' v0 q5 z0 C1 M- _. V1 p- d

plot(year,wolfer)

) Q7 f* j: C' @$ D

title(原始数据)

2 J7 L) C, b2 B' b5 J. u

Y=fft(wolfer);

; B7 J2 I1 t6 q) c! p4 E

N=length(Y);

Y(1)=[];

: Y" Q( t% X1 r& H

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

nyquist=1/2;

4 T+ M* b' A/ T# ]) ]# W

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

8 F# E; J2 I0 L b9 m0 e; q6 t

period=1./freq;

) Z7 M) V G* N% G4 L1 \

subplot(2,1,2)

B5 f. K5 E- ^% G0 z* \

plot(period,power)

, H+ U5 T/ p9 N5 f4 u5 D

title(功率谱), grid on

axis([0 40 0 2e7])

* v9 s, v4 [; q' c

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

8 Y" Q- e) \7 y7 Z/ a1 S- _# t % h; @ t. E2 k l' U c