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

5 ?2 H4 ]) s% A: y: _% Z

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

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

1 }3 N6 B' |0 U+ s# f4 U" h

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

+ b5 R3 g- R1 ?$ ]8 e2 h; A* v

load count.dat

mx=max(count)

mx = 114 145 257

1 A( S. Y) {* e5 [- O/ X

mu=mean(count)

mu = 32.0000 46.5417 65.5833

sigma=std(count)

& n7 U9 p+ v' g; n# b' L2 O

sigma = 25.3703 41.4057 68.0281

5 N4 M! B1 r: [1 r5 Z" _

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

[mx,indx]=min(count)

mx = 7 9 7

indx = 2 23 24

4 i; r: b, v( U6 Y6 J- Z) Q2 ~

1、协方差和相关系数

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

, E& H2 @: ?' T

cv=cov(count)

* o5 G, A a6 j& J/ Z& i

cv = 1.0e+003 *

0.6437 0.9802 1.6567

% U3 k; q4 S% k7 m$ @# i. |

0.9802 1.7144 2.6908

& B3 ]6 M& v# M

1.6567 2.6908 4.6278

cr=corrcoef(count)

cr =

6 x# Q4 f3 e( X2 L% t; p( A2 W

1.0000 0.9331 0.9599

2 T# s# Z- `7 H! c" c

0.9331 1.0000 0.9553

0.9599 0.9553 1.0000

2、数据预处理

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

1 n; Z. n, n4 ?5 T4 @3 ?

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

sum(a)

ans = 13 NaN 13

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

& _! m* S: A- }, i U

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

2 B( _1 S" b, g0 k2 _* s: N( G9 M

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

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

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

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

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

- Z; r4 S6 @( M9 V x

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

4 I( c( |8 M& n8 T4 `

3、回归和曲线拟合

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

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

& g5 L8 `) {+ `7 R% `

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

y=[0.5 0.82 1.14 1.25 1.35 1.40];

1 h- i# X' \% e/ f6 q( N0 z

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

' A/ O) j$ r( v: {! {* R) H0 m! e

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

8 }2 c& y4 ~# w1 I7 F

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

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

* ~6 R$ m0 p$ u% w3 W7 j

a=X1\y;

$ J: J2 \9 d2 W9 G5 \

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

: J B, u9 t4 w& K

b=X2\y;

T=[0:.1:2.5];

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

title(多项式回归)

subplot(1,2,2)

# \ a/ n# J# L* ~: c

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

% `. G# L. y) q* x

title(指数函数回归)

1 `/ J [( d$ g; v/ @6 |; _; ?

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

0 |3 ^" E- H g

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

1 U5 t8 n2 i8 C, Z( K$ [/ `

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

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

采用来拟合，则有

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

/ k- f, z) P: w L7 u) M9 v* P

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

: v1 c' O# s4 i, q" \

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

+ {2 L2 `) I% s( d

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

, F- f7 H+ |" Z, c# F) |9 o9 w) s

a=X\y

- d* D6 a9 ^) w9 ?; y; Z

a = 0.1018 0.4844 −0.2847

因此数据的拟合模型为

. I$ S! y, U8 k! r8 C

y=0.1018+0.4844x1−0.2487x2

4、傅里叶分析与FFT

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

: ^' O0 @' } k/ I% u6 u4 S5 [

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

t=0:1/119:1;

. l, p! _0 _1 g9 h

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

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(多频率混合信号)

ylabel(Input \itx),xlabel(Time )

% w" i/ v# A7 ^/ P" d0 g, v- p

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

- \) \/ t) |' V; [: E4 M* |+ S

ylabel(Abs. Magnitude),grid on

xlabel(Frequency (Hertz))

1 ]6 r' Y/ M# @% p5 N. P7 I! F

5 w9 ]3 A$ X- E

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

t=0:1/199:1;

8 c. f5 f8 t; o3 S3 V

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

y=fft(x);

m=abs(y);

1 j8 \8 M9 p3 d! W" Z! q* k

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

' K' e% q$ M/ r

figure(1)

+ ^8 N' h" \' @. D' l

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

( H6 [3 F4 p6 q: @4 z p

title(信号检测)

ylabel(Input \itx),xlabel(Time )

) E \3 ?& A7 X5 H# R2 n, @. M

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

$ Q, ?* G+ |. |0 D {& \6 n

ylabel(Abs. Magnitude),grid on

! X6 R% U9 k6 M

xlabel(Frequency (Hertz))

1 K8 V/ S& I$ f% i8 h

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

) T; ^, L, R1 a

load sunspot.dat

) r& y8 D* k8 E6 [ b4 v

year=sunspot(:,1);

% o |, a& g8 X

wolfer=sunspot(:,2);

3 m T1 y+ q# d" I4 A( E3 G

figure(1)

0 y/ x' \9 E* b/ C) }1 ]! p

subplot(2,1,1)

' i0 F/ a U7 b

plot(year,wolfer)

9 V e2 A$ r! F3 E( @2 w

title(原始数据)

1 \7 G- Z* _ N' c

Y=fft(wolfer);

N=length(Y);

Y(1)=[];

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

( E* ^' L4 t+ u

nyquist=1/2;

+ F: M7 w! E# `

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

period=1./freq;

" z: G' P4 f9 o1 c

subplot(2,1,2)

plot(period,power)

title(功率谱), grid on

* \( Z# M2 i( U( ~6 i8 `! U1 E

axis([0 40 0 2e7])

# g/ m2 e' v1 U6 h

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

2 A4 q' I* E2 q& x * v' u. b J* }) e9 R! ` 0 I; @' v3 J% h5 R; z0 @