数学建模—大数据统计


阅读数: 次   

EverydayOneCat

卡其脱离太！

知识点

1.大数据的MATLAB导入导出

1.1调用xlsread函数读取数据

常用格式：

  • num = xlsread(filename, sheet, range)
    • sheet可省略，默认是’Sheet1’
    • range是左上角到右下角
    • 切记在全国大学生数学建模赛中不要用绝对路径
  • num = xlsread(filename, -1)
    • 这个是打开excel表让你自己选择区域，也很常用

将数据A.xlsx放入默认路径下，读取文件A.xlsx第1个工作表中单元格A2:H4中的数据

1
num=xlsread('A.xlsx', 'A2:H4')

1.2把数据写入Excel文件

  • xlswrite(filename, M, sheet, range)
    • M是需要插入的数据名称
    • 如果没有该文件，会自动创建一个

把矩阵x写入文件B.xls(放在默认路径下)的第2个工作表中的单元格区域D6:I10，并返回操作信息

1
[s,t] = xlswrite(‘B.xls', x, Sheet2, 'D6:I10‘)

定义一个元胞数组，将它写入Excel文件B.xls的自命名工作表的指定区域；把元胞数组x写入文件D盘的B.xls的指定工作表（Sheet1）中的单元格区域A3:F5

1
2
3
4
5
6
7
x = {1,60101,6010101,'陈亮',63,'';2,60101,6010102,'李旭',73,'';3,60101,...
6010103,'刘鹏飞',0,'缺考'} % 定义一个元胞数组，这种有中文的用大括号
x =
[1] [60101] [6010101] '陈亮' [63] ''
[2] [60101] [6010102] '李旭' [73] ''
[3] [60101] [6010103] '刘鹏飞' [ 0] '缺考'
xlswrite('D:\B.xls',x,'Sheet1','A3:F5')

2.大数据的清洗

2.1缺失值处理：插值

在实际中，常常要处理由实验或测量所得到的一些离散数据。插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。此类问题为插值问题。

MATLAB 实现：实现分段线性插值不需要编制函数程序，它自身提供了内部的功能函数：
interp1 (一维插值) intep2 (二维) interp3 (三维) intern (n维)

image-20200727140428703

例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值。

1
2
3
4
5
6
hours=1:12;
temps=[5 8 9 15 25 29 31 30 22 25 27 24];
h=1:0.1:12;
t=interp1(hours,temps,h,'spline');
plot(hours,temps,'+',h,t,hours,temps,'r:')
xlabel('Hour'),ylabel('Degrees Celsius')

image-20200727141037113

+++

image-20200727141113311

例：测得平板表面3×5网格点处的温度分别为：
82 81 80 82 84
79 63 61 65 81
84 84 82 85 86
试作出平板表面的温度分布曲面z=f(x,y)的图形．

1.先在三维坐标画出原始数据，画出粗糙的温度分布曲线图

1
2
3
4
x=1:5;
y=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
mesh(x,y,temps)

image-20200727141242753

2．以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值，接着画图的图像就更有观测性

1
2
3
4
xi=1:0.2:5;
yi=1:0.2:3;
zi=interp2(x,y,temps,xi,yi','cubic');%这里xi,yi必须有一个要逆
mesh(xi,yi,zi)

image-20200727141446271

+++

image-20200727141502758

例：在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）×（-50，150）里的哪些地方船要避免进入．

image-20200727141610212

做题步骤：

  1. 输入插值基点数据
  2. 在矩形区域(75,200)×(-50,150)进行插值。
  3. 作海底曲面图
  4. 作出水深小于5的海域范围,即z=5的等高线.
1
2
3
4
5
6
7
8
9
10
11
%程序1：插值并作海底曲面图
x =[129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0
162.0 162.0 117.5 ];
y =[ 7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5
84.0 -33.5 ];
z =[ 4 8 6 8 6 8 8 9 9 8 8 9 4 9 ];
x1=75:1:200;
y1=-50:1:150;
[x1,y1]=meshgrid(x1,y1);%平面xy曲线
z1=griddata(x,y,z,x1,y1,'v4');
meshc(x1,y1,z1)

海底曲面图：image-20200727141935022

1
2
3
4
5
6
7
8
9
10
11
12
%程序2：插值并作出水深小于5的海域范围。
x =[129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0
162.0 117.5 ];
y =[ 7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -
33.5 ];
z =[ 4 8 6 8 6 8 8 9 9 8 8 9 4 9 ];
x1=75:1:200;
y1=-50:1:150;
[x1,y1]=meshgrid(x1,y1);
z1=griddata(x,y,z,x1,y1,'v4'); %插值
z1(z1>=5)=nan; %将水深大于5的置为nan，这样绘图就不会显示出来
meshc(x1,y1,z1)

水深小于5的海域范围：image-20200727142044346

2.2异常值处理

异常是在数据集中偏离大部分数据的数据，使人怀疑这些数据的偏离并非由随机因素产生，而是产生于完全不同的机制。

异常挖掘（outlier mining）问题由两个子问题构成：
• (1)如何度量异常；
• (2)如何有效发现异常。

(1)基于统计的异常检测

image-20200727190443750

基于统计方法异常点检测技术的优缺点：

• 优点：
1）异常点检测的统计学方法具有坚实的基础，建立在标准的
统计学技术(如分布参数的估计)之上。
2）当存在充分的数据和所用的检验类型的知识时，这些检验
可能非常有效。
• 缺点：
1）大部分统计方法都是针对单个属性的，对于多元数据技术方法较少。
2）在许多情况下, 数据分布是未知的。
3）对于高维数据, 很难估计真实的分布。

(2)基于聚类的异常检测

物以类聚——相似的对象聚合在一起，基于聚类的异常点检测方法有两个共同特点：
（1）先采用特殊的聚类算法处理输入数据而得到聚类，再在聚类的基础上来检测异常。
（2）只需要扫描数据集若干次，效率较高，适用于大规模数据集。

image-20200727190703883

其中对象p到每个类之间的距离d(p,Ci)有两种计算方法：

  • p与类Ci的重心之间的距离；
  • p与类Ci中每个样本之间的距离的平均值。

以上异常点检测方法也称为两阶段法，也简称为TOD。

基于聚类的方法：具体程序见程序TOD.m

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
function yichang=TOD(A)
N=size(A,1);
D=N;
yichang=zeros(N,1);
K=3;
[u re]=kmeans(A,K);  %K-均值聚类，分为k类——u代表哪一类，re代表重心
tt=sort(u);%排序
%聚类结果
    [t1,t2]=find(u==1);
    class1=A(t1,:);
    [t1,t2]=find(u==2);
    class2=A(t1,:);
    [t1,t2]=find(u==K);
    classK=A(t1,:);
    N1=size(class1,1);
    N2=size(class2,1);
    NK=size(classK,1);
   for i=1:size(A,1)
       D1=0;
       D2=0;
       DK=0;
       for j=1:N1
          D1=sqrt((A(i,:)-class1(j,:))*(A(i,:)-class1(j,:))')+D1; 
       end
       temp1=N1/N*D1;
       for j=1:N2
          D2=sqrt((A(i,:)-class2(j,:))*(A(i,:)-class2(j,:))')+D2; 
       end
       temp2=N2/N*D2;
       for j=1:NK
          DK=sqrt((A(i,:)-classK(j,:))*(A(i,:)-classK(j,:))')+DK; 
       end
       tempK=NK/N*DK;
       op(i)=temp1+temp2+tempK;
   end
xmean=mean(op);
xvar=std(op);
bata=1.28;
for i=1:N
    if op(i)>xmean+bata*xvar
      yichang(i)=1;
    end
end
  yichang

在命令窗口输入

1
A=xlsread('A.xlsx'); yichang=TOD(A)

运行，得到运行结果：yichang显示1的就是异常

作业

1.异常点

用所学的方法找出book_A中的异常点（红色标号）

文件链接：https://pluto-1300780100.cos.ap-nanjing.myqcloud.com/download/Book_A.xlsx

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
A=xlsread('Book_A.xlsx','C2:R790');
N=size(A,1);
D=N;
yichang=zeros(N,1);
K=3;
[u re]=kmeans(A,K);  %K-均值聚类
tt=sort(u);
%聚类结果
    [t1,t2]=find(u==1);
    class1=A(t1,:);
    [t1,t2]=find(u==2);
    class2=A(t1,:);
    [t1,t2]=find(u==K);
    classK=A(t1,:);
    N1=size(class1,1);
    N2=size(class2,1);
    NK=size(classK,1);
   for i=1:size(A,1)
       D1=0;
       D2=0;
       DK=0;
       for j=1:N1
          D1=sqrt((A(i,:)-class1(j,:))*(A(i,:)-class1(j,:))')+D1; 
       end
       temp1=N1/N*D1;
       for j=1:N2
          D2=sqrt((A(i,:)-class2(j,:))*(A(i,:)-class2(j,:))')+D2; 
       end
       temp2=N2/N*D2;
       for j=1:NK
          DK=sqrt((A(i,:)-classK(j,:))*(A(i,:)-classK(j,:))')+DK; 
       end
       tempK=NK/N*DK;
       op(i)=temp1+temp2+tempK;
   end
xmean=mean(op);
xvar=std(op);
bata=1.28;
for i=1:N
    if op(i)>xmean+bata*xvar
      yichang(i)=1;
    end
end
  disp(yichang);

2.城市表层土壤重金属污染分析

随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。

现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。

附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。

对2011年全国建模A的数据，进行插值，画出该地区的地形图，同时绘出各金属污染物的密度分布图

地形图：

1
2
3
4
5
6
7
8
x=xlsread('cumcm2011A附件_数据.xls','附件1','B4:B322');
y=xlsread('cumcm2011A附件_数据.xls','附件1','C4:C322');
z=xlsread('cumcm2011A附件_数据.xls','附件1','D4:D322');
x1=min(x):200:max(x);
y1=min(y):200:max(y);
[x1,y1]=meshgrid(x1,y1);
z1=griddata(x,y,z,x1,y1,'v4');
meshc(x1,y1,z1)

image-20200728104826032

金属污染物密度分布图：

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
x=xlsread('cumcm2011A附件_数据.xls','附件1','B4:B322');
y=xlsread('cumcm2011A附件_数据.xls','附件1','C4:C322');
z=xlsread('cumcm2011A附件_数据.xls','附件1','D4:D322');
x1=min(x):200:max(x);
y1=min(y):200:max(y);
[x1,y1]=meshgrid(x1,y1);
z1=griddata(x,y,z,x1,y1,'v4');
figure(1);meshc(x1,y1,z1);title('地形图')
As=xlsread('cumcm2011A附件_数据.xls','附件2','B4:B322');
Cd=xlsread('cumcm2011A附件_数据.xls','附件2','C4:C322');
Cr=xlsread('cumcm2011A附件_数据.xls','附件2','D4:D322');
Cu=xlsread('cumcm2011A附件_数据.xls','附件2','E4:E322');
Hg=xlsread('cumcm2011A附件_数据.xls','附件2','F4:F322');
Ni=xlsread('cumcm2011A附件_数据.xls','附件2','G4:G322');
Pb=xlsread('cumcm2011A附件_数据.xls','附件2','H4:H322');
Zn=xlsread('cumcm2011A附件_数据.xls','附件2','I4:I322');
As1=griddata(x,y,As,x1,y1,'v4');
figure(2);meshc(x1,y1,As1);title('As分布图')
Cd1=griddata(x,y,Cd,x1,y1,'v4');
figure(3);meshc(x1,y1,Cd1);title('Cd分布图')
Cr1=griddata(x,y,Cr,x1,y1,'v4');
figure(4);meshc(x1,y1,Cr1);title('Cr分布图')
Cu1=griddata(x,y,Cu,x1,y1,'v4');
figure(5);meshc(x1,y1,Cu1);title('Cu分布图')
Hg1=griddata(x,y,Hg,x1,y1,'v4');
figure(6);meshc(x1,y1,Hg1);title('Hg分布图')
Ni1=griddata(x,y,Ni,x1,y1,'v4');
figure(7);meshc(x1,y1,Ni1);title('Ni分布图')
Pb1=griddata(x,y,Pb,x1,y1,'v4');
figure(8);meshc(x1,y1,Pb1);title('Pb分布图')
Zn1=griddata(x,y,Zn,x1,y1,'v4');
figure(9);meshc(x1,y1,Zn1);title('Zn分布图'

As密度分布图image-20200728110945243

Cd密度分布图image-20200728111015903

Cr密度分布图image-20200728111144092

Cu密度分布图image-20200728111155668

Hg密度分布图image-20200728111214321

Ni密度分布图image-20200728111225990

Pb密度分布图image-20200728111243182

Zn密度分布图image-20200728111301019

结语

Let life be beautiful like summer flowers and death like autumn leaves

要为六级做准备了😥

赏

谢谢你请我吃糖果

支付宝
微信

扫一扫，分享到微信

微信分享二维码

tag:

    缺失模块。
    1、请确保node版本大于6.2
    2、在博客根目录（注意不是yilia根目录）执行以下命令：
    npm i hexo-generator-json-content --save

    3、在根目录_config.yml里添加配置： 

      jsonContent:
    meta: false
    pages: false
    posts:
      title: true
      date: true
      path: true
      text: false
      raw: false
      content: false
      slug: false
      updated: false
      comments: false
      link: false
      permalink: false
      excerpt: false
      categories: false
      tags: true

很惭愧<br><br>只做了一点微小的工作<br>谢谢大家