Matlab数据分析

数学建模 Comments

目录

1 多项式及其函数

1.1 多项式的创建与表达

· 输入多项式 $2x^4-85x^3+15x^2+100$

1
a=ploy2sym([2 -85 15 100])

1.2 多项式的四则运算

1.2.1 加、减法

1
2
3
4
5
a=[1 2 3 4 5];
b=[2 3 4 0 5];
c=a+b;
c
poly2sym(c)        %显示多项式

4-1


1.2.2 乘法

1
2
3
a=[1 2];
b=[3 4];
c=poly2sym(conv(a,b))

4-2


1.2.3 除法

1
2
3
4
5
a=[1 4 9 15];
b=[1 2 -3 20 -16 21];
[q,r]=deconv(b,a)
poly2sym(q)
poly2sym(r)

$[q,r]=deconv(b,a)$:
  q为多项式b除以a的商户式；r为多项式b除以a的余式。返回的q,r仍为多项式系数向量。

4-3


1.3 多项式求导

k = polyder(p) : 求多项式p的导函数多项式。
k = polyder(a,b) : 求多项式a和b 乘积的导函数多项式。
[p,q] =polyder(a,b):求多项式b与多项式a相除的导函数多项式，导函数的分子存在p中，分母存在q中。其中，参数a,b是多项式的系数向量，返回结果p,q也是多项式的系数向量。

1
2
3
4
a=[2 4 6 -8 10]
poly2sym(a)
b=polyder(a)
poly2sym(b)

4-4



· 已知多项式$p(x)=2x^3-x^2+3$和$q(x)=2x=1$，求如下导数：$p’、(p·q)’、(p/q)’$。 

1
2
3
4
5
6
clear all;
k1=polyder([2,-1,0,3]);
poly2sym(k1)
k2=polyder([2,-1,0,3],[2,1]);
poly2sym(k2)
[k3,d]=polyder([2,-1,0,3],[2,1])

4-5



1.4 多项式求值

1.4.1 polyval函数

· 已知多项式 $x^4+6x^3-9$ ，分别取$x=1.3$和一个给定的矩阵2x3个元素为自变量，计算该多项式的值。

1
2
3
4
5
6
7
8
clear all;
A=[1 6 0 0 -9];     %定义的多项式
x=1.3;              %定义的自变量
y1=polyval(A,x);    
y1
B=[-1 1.3 1.5;2 -1.8 1.6];      %定义的矩阵为自变量
y2=polyval(A,B);
y2

· 调用此函数时，有两种形式：
  $y=poluval(A，x)$ —其中x为单独的自变量
  $y=poluval(A，B)$ —其中B为一个自变量矩阵
4-6


1.4.2 polyvalm函数

· 已知多项式 $p(x)=3^3-x^2+5$ ，求矩阵 A=$\begin{pmatrix} -1 & 2 \\ -2 & 1 \end{pmatrix}$ 的多项式的值。

1
2
3
4
5
clear all;
p=[3 -1 0 5];
x=[-1 2; -2 1];
polyval(p,x)
polyvalm(p,x)

4-7


1.5 多项式求根

· 已知多项式 $x^3-6x^2-72x-27$ 的根。

1
2
3
4
clear all;
p=[1 -6 -72 -27];
k=poly2sym(p)
r=roots(p)

4-8

· 若已知多项式的根，可以反求多项式

1
2
3
4
5
r=[2 5 8 7];    %已知根
p=poly(r)       %反求多项式系数
k=poly2sym(p)   %导出多项式
r=roots(p)      %再次求出多项式的根验证正确性


· 注意有根才能求根！
4-9


1.6 部分分式展开

对有理多项式 $\frac{10(s+2)}{(s+1)(s+3)(s+4)}$ 进行展开

1
2
3
4
clear all;
num=10*[1 2];
den=poly([-1;-3;-4]);
[res,poles,k]=residue(num,den)

结果是余数、极点和部分分式展开的常数项，如下：

$\frac{10(s+2)}{(s+1)(s+3)(s+3)}=\frac{-6.66667}{s+4}+\frac{5}{s+3}+\frac{1.6667}{s+1}$

4-10


1.7 求解线性方程组

$\begin{cases} 2x+3y-z=2 \\8x-2y-3z=4 \\ 45x+3y+9z=23 \end{cases}$

方法一、

1
2
3
a=[2 3 -1;8 2 3;45 3 9];
b=[2;4;23];
z=rats(inv(a)*b)

方法二：

1
2
syms x y z
[x y z]=solve(2*x+3*y-z-2,8*x+2*y+3*z-4,45*x+3*y+9*z-23)


1.8 解定积分

$I= \int_0^1 xlin(1+x)d_x$

· 方法一：

1
quad('x.*log(1+x)',0,1)

· 方法二：
1
2
syms x
int(x*log(1+x),0,1)

1.9 曲面拟合

最小二乘拟合

· 某地一煤矿，为估计其储量便于开采，先在该地进行开采，先在该地进行勘探。假设该地是一长方形区域，长为4km，宽为5km.勘探得到如下表所示数据：

编号 坐标/km 煤层厚度/km * 编号 坐标/km 煤层厚度/km
1 (1,1) 13.2 * 11 (3,1) 23.28
2 (1,2) 25.80 * 12 (3,2) 26.48
3 (1,3) 8.47 * 13 (3,3) 29.14
4 (1,4) 25.27 * 14 (3,4) 12.04
5 (1,5) 22.32 * 15 (3,5) 14.58
6 (2,1) 15.47 * 16 (4,1) 19.95
7 (2,2) 21.33 * 17 (4,2) 23.73
8 (2,3) 14.49 * 18 (4,3) 15.35
9 (2,4) 24.83 * 19 (4,4) 18.01
10 (2,5) 26.19 * 20 (4,5) 16.29 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
x=[1 2 3 4 5 6 7 8 9 10]
y=[1 2 3 4 4 4.7 5 5.2 6 7.2 ]


%尝试一次多项式拟合  
p1=polyfit(x,y,1)

%尝试第二次拟合  
p3=polyfit(x,y,3)

%对拟合曲线进行比较
x2=1:0.1:10
y1=polyval(p1,x2)
y3=polyval(p3,x2)
plot(x,y,'*',x2,y1,'b',x2,y3,'r')

拟合效果

4-12


采用拟合的方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

  clc;
 x=1:4;
 y=1:5;
 [X,Y]=meshgrid(x,y);
 Z=-[13.72 25.80 8.47 25.27 22.32;15.47 21.33 14.49 24.83 26.19;23.28 26.48 29.14 12.04 14.58;19.95 23.73 15.35 18.01 16.29]';
 surf(X,Y,Z);
 x=reshape(X,20,1);
 y=reshape(Y,20,1);
 z=reshape(Z,20,1);
 p=4;
 q=5;
 U=leftmatrix(x,p,y,q);
 V=rightmatrix(x,p,y,q,z);
 a_n=U\V;
 aa=zeros();
 for i=1:length(a_n)
     ii=quotient(i-1,q)+1;
     jj=mod(i-1,q)+1;
     aa(ii,jj)=a_n(i,1);
 end
 
 aa;
 
 m=31;n=41;
 [X1,Y1]=meshgrid(linspace(1,4,m),linspace(1,5,n));
 xx=reshape(X1,m*n,1);
 yy=reshape(Y1,m*n,1);
 zz=zeros(m*n,1);
 xy=zeros(m*n,1);
 xt=zeros(m*n,1);
 yt=zeros(m*n,1);
 for i=1:p
     for j=1:q
         xt=xx.^(i-1);
         yt=yy.^(j-1);
         xy=xt.*yt;
         zz=zz+aa(i,j).*xy;
     end
 end
 Z1=reshape(zz,n,m);
 surf(X1,Y1,Z1);
 p=sum(Z);
 q=p(:,[2,3,4]);
 h=sum(q')/15;
 v=2000*4000*h;
 
 function U =leftmatrix(x,p,y,q)
     %U*a=V,a为系数列矩阵，长度为p*q;
     %U为q*p矩阵；
     %x,y为长度一致的列矩阵,给定点的坐标；
     %p,q为拟合的函数的x,y的幂的最高次数；
     m=length(x);
     if(nargin~=4)&(m~=length(y))
         error('错误，检查错误！！');
     end
     U_L=p*q;%U为p*q阶方程
     U=zeros(U_L,U_L);
     for i=1:p*q
         for j=1:p*q
             x_z=quotient(j-1,q)+quotient(i-1,q); %x的幂的次数，quotient为求商
             y_z=mod(j-1,q)+mod(i-1,q);
             U(i,j)=qiuhe(x,x_z,y,y_z);
         end
     end
 end
 

 function V=rightmatrix(x,p,y,q,z)
     %U*a=V
     %V为一个列向量，长为p*q
     %x,y,z为点的坐标
     %p,q分别为x,y的幂的最高次数

     if nargin~=5
         error('错误，检查检查rightmatrix');
     end
     V=zeros(p*q,1);
     for i=1:p*q
         x_z=quotient(i-1,q);
         y_z=mod(i-1,q);
         V(i,1)=qiuhe(x,x_z,y,y_z,z);
     end
 end

 
 
 
function sh =quotient(x,y)
    %sh为x/y的商
    sh=(x-mod(x,y)/y);
end
    
function he =qiuhe(x,p,y,q,z)
    %x,y向量长度相同
    %p,q分别为x,y的幂的次数
    m=length(x);
    if (nargin<4)&&(m~=length(y))
        error('错误，检查!!')
    end
    if nargin ==4
        z=ones(m,1);
    end
    he=0;
    for i=1:m
        he=he+x(i)^p*y(i)^q*z(i);
    end
end

出了问题，暂时搁置，以后解决！！



数值插值

函数的极限

数值积分

多元统计分析

假设检验

7  回归分析

  在客观世界中普遍存在着变量之间的关系，像某个现象的发生或某种结果的得出，往往与其他某个或某些因素有关，但这种关系又不是很确定的，只是从数据上可以看出“有关”的趋势。回归分析就是用来研究区域，有这种特征的变量之间的线性关系，用一定的线性回归模拟来你和因变量和自变量的数据，并通过确定模型参数得到回归方程.


7.1一元线性回归

  当线性回归中的自变量只有一个时，称为一元线性回归假设，对于x(在某个区间内)的每一个值有:

$y - N(a+bx,\sigma^2)$

其中, $a、b、\sigma^2$ 都是不依赖于 $x$ 的未知参数。记对于做这样的正态假设相当于假设

$\begin{cases} Y=a+bx+\epsilon\\\epsilon-N(0,\sigma^2)\end{cases}$

其中，未知参数都不依赖 $x$。
  通常采用最小二乘法来确定。上面两个待定参数 $a$ 和 $b$ 及要求观测值与利用上面回归模型得到的拟合值之间的差值的平方和最小。差值平方和达到最小时的模型参数便作为待定参数的最终值，代入模型，即可得到回归方程。
  在Matlab中，提供了polyfit函数来实现从一次到高次多项式的回归法。polyfit函数在前面章节已经介绍过，在此不再介绍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
%% 输入数据
clc
clear all
x=[23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40];
y=[41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0];	
% 绘制散点图，判断是否具有线性关系
figure
plot(x,y,'r*')%作散点图
xlabel('x（职工工资总额）','fontsize', 12)           %横坐标名
ylabel('y（商品零售总额）', 'fontsize',12)           %纵坐标名
set(gca,'linewidth',2);
%% 采用LinearModel.fit函数进行回归

m2 = LinearModel.fit(x,y)
%% 采用regress函数进行回归
Y=y';
X=[ones(length(y),1),x'];
[b, bint, r, rint, s] = regress(Y, X)
%绘图
b1=b(1);b2=b(2);
xx=min(x):0.01:max(x);
yy=b1+b2*xx;
plot(x,y,'r*',xx,yy,'b')
xlabel('x（职工工资总额）','fontsize', 12)           %横坐标名
ylabel('y（商品零售总额）', 'fontsize',12)           %纵坐标名
set(gca,'linewidth',2);

4-13

  • 采用LinearModel.fit函数进行回归
  • $[b, bint, r, rint, s] = regress(Y, X)$函数用于另外一种回归!

4-14


非线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
%用nlinfit或fitnlm
%% 输入数据
clc, clear all, close all
x=[1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];
y=[7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];  
plot(x,y,'*','linewidth',2);
set(gca,'linewidth',2);
xlabel('销售额x/万元','fontsize', 12)         
ylabel('流通费率y/%', 'fontsize',12)           


%% 对数形式
m1 = @(b,x) b(1) + b(2)*log(x);
nonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])
b=nonlinfit1.Coefficients.Estimate;
Y1=b(1,1)+b(2,1)*log(x);
hold on
plot(x,Y1,'--k','linewidth',2)
%% 指数形式拟合
m2 = 'y ~ b1*x^b2';
nonlinfit2 = fitnlm(x,y,m2,[1;1])
b1=nonlinfit2.Coefficients.Estimate(1,1);
b2=nonlinfit2.Coefficients.Estimate(2,1);
Y2=b1*x.^b2;
hold on
plot(x,Y2,'r','linewidth',2)
legend('原始数据','a+b*lnx','a*x^b')

4-15


逐步回归

4-16

TuTo

天地与我并生，万物与我为一。