研究生数学建模比赛感触

研究生数学建模比赛感触

一.题目

本次研究生建模比赛，我深有感触，我们队友一起选的B题，题目是关于空气质量的二次建模问题，此题明显是考察数据分析的，我做两问，以及第三问借鉴一下别人的想法。总体感觉难度有点大，因为是第一次参加，对模型不够熟练，只做一部分，不过也收获了不少，比如数据预处理，表的合并，写入文件，以及缺失值处理，还有统计方法，除此外还用到模型。主要是python的pandas,numpy以及matplotlib库的使用，还有聚类（KMeans），逻辑回归，支持向量机的应用.

题目

二.建模过程以及代码、图片、数据等

问题一 这个应该就是一个公式计算问题

代码如下（使用 jupyter notebook）:

代码语言：javascript复制

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import pandas as pd 
import numpy as np
import math


# In[2]:


data = pd.read_excel("附件1 监测点A空气质量预报基础数据.xlsx",sheet_name=2)
data1 =data[497:501]
data1


# In[3]:


#得到各个时间段的浓度值
so2 = data1.iloc[:,2].tolist()
no2 = data1.iloc[:,3].tolist()
pm10 = data1.iloc[:,4].tolist()
pm25 = data1.iloc[:,5].tolist()
o3 = data1.iloc[:,6].tolist()
co = data1.iloc[:,7].tolist()


# In[4]:


#co的AOI计算函数
def co_js(co):
    if co >= 0 and co <2:
        BPhi=2
        BPlo=0
        IAQIhi=50
        IAQIlo=0
        Cp=co
    elif co >= 2 and co <4:
        BPhi=4
        BPlo=2
        IAQIhi=100
        IAQIlo=50
        Cp=co
    elif co >= 4 and co < 14:
        BPhi=14
        BPlo=4
        IAQIhi=100
        IAQIlo=150
        Cp=co
    elif co >= 14 and co < 24:
        BPhi=24
        BPlo=14
        IAQIhi=200
        IAQIlo=150
        Cp=co
    elif co >= 24 and co < 36:
        BPhi=36
        BPlo=24
        IAQIhi=300
        IAQIlo=200
        Cp=co
    elif co >=36 and co < 48:
        BPhi=48
        BPlo=36
        IAQIhi=400
        IAQIlo=300
        Cp=co
    elif co >= 48 and co < 60:
        BPhi=60
        BPlo=48
        IAQIhi=500
        IAQIlo=400
        Cp=co
    else:
        IAQIp=0
    IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo) IAQIlo))
    return IAQIp


# In[5]:


#co的AOI计算结果
co_IAQI=[]
for i in co:
    co_IAQI.append(co_js(i))
co_IAQI


# In[6]:


def so2_js(so2):
    if so2 >= 0 and so2 < 50:
        BPhi=50
        BPlo=0
        IAQIhi=50
        IAQIlo=0
        Cp=so2
    elif so2 >= 50 and so2 < 150:
        BPhi=150
        BPlo=50
        IAQIhi=100
        IAQIlo=50
        Cp=so2
    elif so2 >= 150 and so2 < 475:
        BPhi=475
        BPlo=150
        IAQIhi=150
        IAQIlo=100
        Cp=so2
    elif so2 >= 475 and so2 < 800:
        BPhi=800
        BPlo=475
        IAQIhi=200
        IAQIlo=150
        Cp=so2
    elif so2 >= 800 and so2 < 1600:
        BPhi=1600
        BPlo=800
        IAQIhi=300
        IAQIlo=200
        Cp=so2
    elif so2 >= 1600 and so2 < 2100:
        BPhi=2100
        BPlo=1600
        IAQIhi=400
        IAQIlo=300
        Cp=so2
    elif co >= 2100 and co < 2620:
        BPhi=2620
        BPlo=2100
        IAQIhi=500
        IAQIlo=400
        Cp=so2
    else:
        IAQIp=0
    IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo) IAQIlo))
    return IAQIp


# In[7]:


so2_IAQI=[]
for i in so2:
    so2_IAQI.append(so2_js(i))
so2_IAQI


# In[8]:


def no2_js(no2):
    if no2 >= 0 and no2 < 40:
        BPhi=40
        BPlo=0
        IAQIhi=50
        IAQIlo=0
        Cp=no2
    elif no2 >= 40 and no2 < 80:
        BPhi=80
        BPlo=40
        IAQIhi=100
        IAQIlo=50
        Cp=no2
    elif no2 >= 80 and no2 < 180:
        BPhi=180
        BPlo=80
        IAQIhi=150
        IAQIlo=100
        Cp=no2
    elif no2 >= 180 and so2 < 280:
        BPhi=280
        BPlo=180
        IAQIhi=200
        IAQIlo=150
        Cp=no2
    elif no2 >= 280 and no2 < 565:
        BPhi=565
        BPlo=280
        IAQIhi=300
        IAQIlo=200
        Cp=no2
    elif no2 >= 565 and no2 < 750:
        BPhi=750
        BPlo=565
        IAQIhi=400
        IAQIlo=300
        Cp=no2
    elif co >= 750 and co < 940:
        BPhi=940
        BPlo=750
        IAQIhi=500
        IAQIlo=400
        Cp=no2
    else:
        IAQIp=0
    IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo) IAQIlo))
    return IAQIp


# In[9]:


no2_IAQI=[]
for i in so2:
    no2_IAQI.append(no2_js(i))
no2_IAQI


# In[10]:


def pm10_js(pm10):
    if pm10>= 0 and pm10 < 50:
        BPhi=50
        BPlo=0
        IAQIhi=50
        IAQIlo=0
        Cp=pm10
    elif pm10 >= 50 and pm10 < 150:
        BPhi=150
        BPlo=50
        IAQIhi=100
        IAQIlo=50
        Cp=pm10
    elif pm10 >= 150 and pm10 < 250:
        BPhi=250
        BPlo=150
        IAQIhi=150
        IAQIlo=100
        Cp=pm10
    elif pm10 >= 250 and pm10 < 350:
        BPhi=350
        BPlo=250
        IAQIhi=200
        IAQIlo=150
        Cp=pm10
    elif pm10 >= 350 and pm10 < 420:
        BPhi=420
        BPlo=360
        IAQIhi=300
        IAQIlo=200
        Cp=pm10
    elif pm10 >= 420 and pm10< 500:
        BPhi=500
        BPlo=420
        IAQIhi=400
        IAQIlo=300
        Cp=pm25
    elif pm10 >= 500 and pm10 < 600:
        BPhi=600
        BPlo=500
        IAQIhi=500
        IAQIlo=400
        Cp=pm10
    else:
        IAQIp=0
    IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo) IAQIlo))
    return IAQIp


# In[11]:


pm10_IAQI=[]
for i in pm10:
    pm10_IAQI.append(pm10_js(i))
pm10_IAQI


# In[12]:


def pm25_js(pm25):
    if pm25 >= 0 and pm25 < 35:
        BPhi=35
        BPlo=0
        IAQIhi=50
        IAQIlo=0
        Cp=pm25
    elif pm25 >= 35 and pm25 < 75:
        BPhi=75
        BPlo=35
        IAQIhi=100
        IAQIlo=50
        Cp=pm25
    elif pm25 >= 75 and pm25 < 115:
        BPhi=115
        BPlo=75
        IAQIhi=150
        IAQIlo=100
        Cp=pm25
    elif pm25 >= 115 and pm25 < 150:
        BPhi=350
        BPlo=250
        IAQIhi=200
        IAQIlo=150
        Cp=pm25
    elif no2 >= 150 and no2 < 250:
        BPhi=250
        BPlo=150
        IAQIhi=300
        IAQIlo=200
        Cp=pm25
    elif no2 >= 250 and no2 < 350:
        BPhi=350
        BPlo=250
        IAQIhi=400
        IAQIlo=300
        Cp=pm25
    elif co >= 350 and co < 500:
        BPhi=500
        BPlo=350
        IAQIhi=500
        IAQIlo=400
        Cp=pm25
    else:
        IAQIp=0
    IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo) IAQIlo))
    return IAQIp


# In[13]:


pm25_IAQI=[]
for i in pm25:
    pm25_IAQI.append(pm25_js(i))
pm25_IAQI


# In[14]:


def o3_js(o3):
    if o3 >= 0 and o3 < 100:
        BPhi=100
        BPlo=0
        IAQIhi=50
        IAQIlo=0
        Cp=o3
    elif o3 >= 100 and o3 < 160:
        BPhi=160
        BPlo=100
        IAQIhi=100
        IAQIlo=50
        Cp=o3
    elif o3 >= 160 and o3 < 215:
        BPhi=115
        BPlo=75
        IAQIhi=150
        IAQIlo=100
        Cp=o3
    elif o3 >= 215 and o3 < 265:
        BPhi=350
        BPlo=250
        IAQIhi=200
        IAQIlo=150
        Cp=o3
    elif no2 >= 265 and no2 < 800:
        BPhi=800
        BPlo=265
        IAQIhi=300
        IAQIlo=200
        Cp=o3
    else:
        IAQIp=0
    IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo) IAQIlo))
    return IAQIp


# In[15]:


o3_IAQI=[]
for i in o3:
   o3_IAQI.append(o3_js(i))
o3_IAQI


# In[16]:


#整合数据
AQIA = np.transpose(np.array([co_IAQI,so2_IAQI,no2_IAQI,o3_IAQI,pm10_IAQI,pm25_IAQI]))
AQIA


# In[17]:


#统计浓度最大值
AQIA_end = AQIA.max(axis=1)
AQIA_end


# In[18]:


AQIA_index = np.argmax(AQIA,axis=1)
AQIA_index


# In[19]:


#首要污染物名称
AQIA_columns = []
columns = ["co","so2","no2","o3","pm10","pm25"]
for i in AQIA_index:
    AQIA_columns .append(columns[i])
AQIA_columns


# In[20]:


#首要污染物
wrw = []
for i in range(4):
    if AQIA_end[i] <= 50:
        wrw.append("无")
    else:
        wrw.append(AQIA_columns[i])
wrw


# In[21]:


sj = data1.iloc[:,0]
sj = np.array(sj,dtype="datetime64[D]")
dd = data1.iloc[:,1]


# In[22]:


sj


# In[23]:


#写入文件
columns =["监测日期","地点","AQI","首要污染物"]
DATA = pd.DataFrame(data=np.array([sj,dd,AQIA_end,wrw]).T,columns=columns)
DATA.to_excel("AQI.xlsx")

运行结果如图：

	监测日期	地点	AQI	首要污染物
0	2020-08-25	监测点A	60
1	2020-08-26	监测点A	46	无
2	2020-08-27	监测点A	218
3	2020-08-28	监测点A	258

问题二：做的分类，不知道是否正确，分成六类

代码语言：javascript复制

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import pandas as pd 
import matplotlib.pyplot as plt
import numpy as np 
from sklearn.linear_model import LinearRegression


# In[2]:


data = pd.read_excel("附件1 监测点A空气质量预报基础数据.xlsx",sheet_name=1)
data = data.replace("—",np.nan)
data1 = data.dropna(axis=0,how="any")


# In[3]:


data1


# In[4]:


so2 = np.array(data1.iloc[:,2])
no2 = np.array(data1.iloc[:,3])
pm10 = np.array(data1.iloc[:,4])
pm25 = np.array(data1.iloc[:,5])
o3 = np.array(data1.iloc[:,6])
co = np.array(data1.iloc[:,7])
wd = np.array(data1.iloc[:,8])
sd = np.array(data1.iloc[:,9])
qy = np.array(data1.iloc[:,10])
fs = np.array(data1.iloc[:,11])
fx = np.array(data1.iloc[:,12])


# In[5]:


y = [so2,no2,pm10,pm25,o3,co]
colors =["red","black","orange","blue","green","purple"]
list = ["so2","no2","pm10","pm25","o3","co"]


# In[6]:


for i in range(6):
    plt.figure(figsize=(20,8))
    plt.subplot(3,2,i 1)
    plt.scatter(wd,y[i],c=colors[i])
    plt.xlabel("wd")
    plt.ylabel(list[i])
    plt.savefig("./wd/" list[i] "_wd" ".jpg")
    plt.show()


# In[7]:


for i in range(6):
    plt.figure(figsize=(20,8))
    plt.subplot(3,2,i 1)
    plt.scatter(sd,y[i],c=colors[i])
    plt.xlabel("sd")
    plt.ylabel(list[i])
    plt.savefig("./sd/" list[i] "_sd" ".jpg")
    plt.show()


# In[8]:


for i in range(6):
    plt.figure(figsize=(20,8))
    plt.subplot(3,2,i 1)
    plt.scatter(qy,y[i],c=colors[i])
    plt.xlabel("qy")
    plt.ylabel(list[i])
    plt.savefig("./qy/" list[i] "_qy" ".jpg")
    plt.show()


# In[9]:


for i in range(6):
    plt.figure(figsize=(20,8))
    plt.subplot(3,2,i 1)
    plt.scatter(fx,y[i],c=colors[i])
    plt.xlabel("fx")
    plt.ylabel(list[i])
    plt.savefig("./fx/" list[i] "_fx" ".jpg")
    plt.show()


# In[10]:


for i in range(6):
    plt.figure(figsize=(20,8))
    plt.subplot(3,2,i 1)
    plt.scatter(fs,y[i],c=colors[i])
    plt.xlabel("fs")
    plt.ylabel(list[i])
    plt.savefig("./fs/" list[i] "_fs" ".jpg")
    plt.show()


# In[11]:


from sklearn.decomposition import PCA


# In[12]:


columns = ["wd","sd","qy","fs","fx"]
tzdata = np.array([wd,sd,qy,fs,fx]).T
type(tzdata)


# In[13]:


pac_sum = []
for i in range(1,6):
    pca = PCA(i).fit(tzdata)
    pac_sum.append(pca.explained_variance_ratio_.sum())
pac_sum


# In[14]:


wd1 = pac_sum[0]
sd1 = pac_sum[1]-pac_sum[0]
qy1 = pac_sum[2]-pac_sum[1]
fs1 = pac_sum[3]-pac_sum[2]
fx1 = pac_sum[4]-pac_sum[3]
list = [wd1,sd1,qy1,fs1,fx1]
list


# In[15]:


import numpy as np
plt.rcParams['font.family']=['sans-serif']
plt.rcParams['font.sans-serif']=['SimHei']
plt.figure(figsize=(20,8))
plt.plot(range(1,6),pac_sum)
plt.xticks([1,2,3,4,5])
plt.ylim(0.9,1.01)
plt.xlabel("选择特征个数")
plt.ylabel("方差贡献率误差总和")
plt.savefig("tz_" "var_sum" ".jpg")
plt.show()


# In[16]:


from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score,silhouette_samples


# In[17]:


#查看样本标签
x = KMeans(5).fit(co.reshape(-1,1))
x.labels_


# In[18]:


#查看质心
x.inertia_


# In[19]:


x.cluster_centers_


# In[20]:


#测试打分情况
silhouette_score(co.reshape(-1,1),x.labels_)


# In[21]:


co.reshape(-1,1)


# In[22]:


silhouette_samples(co.reshape(-1,1),x.labels_)


# In[23]:


#确定分类的个数，为5-6类最好
list_inertia = []
for i in range(2,10):
    x = KMeans(n_clusters=i).fit(co.reshape(-1,1))
    list_inertia.append(x.inertia_)
plt.figure(figsize=(25,8))
plt.plot(range(2,10),list_inertia)
plt.xlabel("n")
plt.ylabel("inertia")
plt.savefig("km_" "co" ".jpg")
plt.show()


# In[24]:


def cs_ct(sp):
     x = KMeans(6).fit(sp.reshape(-1,1)).labels_
     return x


# In[25]:


cs_ct(co)


# In[26]:


cl = []
for i in cs_ct(co):
    cl.append(str(i))
cl


# In[27]:


m = cs_ct(co)
m.shape


# In[28]:


co


# In[29]:


#测试区间最大值以及最小值
data = pd.Series(data=co,index=cs_ct(co))
data[1].max()
data[1].min()


# In[30]:


np.array(data[0])


# In[31]:


type(data)


# In[32]:


#获取分类区间
def qj_lb(sp):
    x = KMeans(6).fit(sp.reshape(-1,1)).labels_
    data = pd.Series(data=sp,index=x)
    list1 = []
    list2 = []
    for i in range(6):
        list1.append(data[i].min())
        list2.append(data[i].max())
    list3 = []
    for i in range(6):
        print("第{}个类别的范围为：{}到{}".format(i 1,list1[i],list2[i]))
        #list3.append("第{}个类别的范围为：{}到{}".format(i 1,list1[i],list2[i]))
# for i in range(6):
# plt.figure(figsize=(20,8))
# plt.subplot(3,2,i 1)
# plt.plot(np_sj,sp)
# plt.xlabel("lb")
# plt.ylabel("fw")
# plt.show()
    print("-"*50)


# In[33]:


#遍历循环，打印分类数据
sj =[so2,no2,pm10,pm25,o3,co,wd,sd,qy,fs,fx]
bq =["so2","no2","pm10","pm25","o3","co","wd","sd","qy","fs","fx"]
for i in range(len(bq)):
    print(bq[i] "浓度影响")
    qj_lb(sj[i])
    #写入文件格式不对，未写入，故注释
    #f = open(r"./nd.txt","w")
    #f.write(str(qj_lb(sj[i])))


# In[112]:


def zx_wz(sp):
    x = KMeans(6).fit(sp.reshape(-1,1)).labels_
    data = pd.Series(data=sp,index=x)
    list1 = []
    list2 = []
    list3=[]
    for i in range(6):
        list1.append(data[i].min())
        list2.append(data[i].max())
    return (list1,list2)


# In[ ]:





# In[84]:


zx_wz(co)


# In[79]:


def zx_d(sj):
    list=[]
    for i in range(6):
        list.append(zx_wz(sj)[1][i]-zx_lb(co)[0][i])
    return list


# In[87]:


zx_d(co)


# In[111]:


colors = ["red","blue","green","yellow","black","purple","orange","pink","coral","gold"]
for i in range(10):
    plt.figure(figsize=(60,25))
    plt.subplot(2,5,i 1)
    plt.scatter(range(1,7),zx_d(sj[i]),c=colors[i])
    plt.xlabel(bq[i])
    plt.ylabel("分类质心位置")
    plt.title(bq[i] "分类图")
    plt.savefig("./分类/" bq[i] "分类.jpg")
    plt.grid()
    plt.show()

图片展示：

方差

特征

代码语言：javascript复制

so2浓度
第1个类别的范围为：7.0到9.0
第2个类别的范围为：-2.0到4.0
第3个类别的范围为：18.0到47.0
第4个类别的范围为：13.0到17.0
第5个类别的范围为：5.0到6.0
第6个类别的范围为：10.0到12.0
--------------------------
no2浓度
第1个类别的范围为：33.0到51.0
第2个类别的范围为：77.0到112.0
第3个类别的范围为：2.0到18.0
第4个类别的范围为：52.0到76.0
第5个类别的范围为：113.0到211.0
第6个类别的范围为：19.0到32.0
--------------------------
pm10浓度
第1个类别的范围为：82.0到116.0
第2个类别的范围为：25.0到40.0
第3个类别的范围为：-4.0到24.0
第4个类别的范围为：41.0到58.0
第5个类别的范围为：59.0到81.0
第6个类别的范围为：117.0到217.0
--------------------------
pm25浓度
第1个类别的范围为：25.0到37.0
第2个类别的范围为：-5.0到12.0
第3个类别的范围为：54.0到79.0
第4个类别的范围为：80.0到163.0
第5个类别的范围为：13.0到24.0
第6个类别的范围为：38.0到53.0
--------------------------
o3浓度
第1个类别的范围为：188.0到405.0
第2个类别的范围为：-2.0到22.0
第3个类别的范围为：84.0到128.0
第4个类别的范围为：51.0到83.0
第5个类别的范围为：129.0到187.0
第6个类别的范围为：23.0到50.0
--------------------------
co浓度
第1个类别的范围为：0.7到0.8
第2个类别的范围为：1.1到1.4
第3个类别的范围为：0.1到0.5
第4个类别的范围为：0.9到1.0
第5个类别的范围为：1.5到2.5
第6个类别的范围为：0.6到0.6
--------------------------
wd浓度
第1个类别的范围为：23.9到27.6
第2个类别的范围为：14.8到19.6
第3个类别的范围为：31.2到38.2
第4个类别的范围为：19.7到23.8
第5个类别的范围为：5.8到14.7
第6个类别的范围为：27.7到31.1
--------------------------
sd浓度
第1个类别的范围为：75.0到84.0
第2个类别的范围为：42.0到56.0
第3个类别的范围为：85.0到99.0
第4个类别的范围为：14.0到41.0
第5个类别的范围为：57.0到66.0
第6个类别的范围为：67.0到74.0
--------------------------
qy浓度
第1个类别的范围为：1003.0到1006.9
第2个类别的范围为：1011.3到1015.5
第3个类别的范围为：1007.0到1011.2
第4个类别的范围为：1020.1到1029.2
第5个类别的范围为：1015.6到1020.0
第6个类别的范围为：993.5到1002.9
--------------------------
fs浓度
第1个类别的范围为：0.7到1.0
第2个类别的范围为：2.1到2.6
第3个类别的范围为：1.1到1.5
第4个类别的范围为：2.7到5.8
第5个类别的范围为：0.1到0.6
第6个类别的范围为：1.6到2.0
--------------------------
fx浓度
第1个类别的范围为：40.4到84.0
第2个类别的范围为：261.7到311.7
第3个类别的范围为：174.1到261.6
第4个类别的范围为：84.1到172.8
第5个类别的范围为：0.0到40.3
第6个类别的范围为：311.8到360.0
--------------------------

后面的没做，就不展示。

不管咋样，这次竞赛收获还是很多的，有付出才有收获，加油。

python numpy jupyter matplotlib pandas

0 人点赞