第三章 数学建模

一、数学建模概述
1.数学建模算法内容

2.线性回归

3.KNN最邻近分类

4.PCA主成分分析

5.K-means聚类

6.蒙塔卡罗模拟

二、线性回归
1.线性回归理论概述



2.线性回归的python实现方法
线性回归通常是人们在学习预测模型时首选的技术之一。在这种技术中，因变量是连续的，自变量可以是连续的也可以是离散的，回归线的性质是线性的。
线性回归使用最佳的拟合直线（也就是回归线）在因变量（Y）和一个或多个自变量（X）之间建立一种关系
① 数据示例

from sklearn.linear_model import LinearRegression   # 导入线性回归模块
rng = np.random.RandomState(1)   # np.random.RandomState → 随机数种子，对于一个随机数发生器，只要该种子（seed）相同，产生的随机数序列就是相同的
xtrain = 10 * rng.rand(30)
ytrain = 8 + 4 * xtrain + rng.rand(30)   # 样本关系：y = 8 + 4*x
xtest = np.linspace(0,10,1000)   # 创建测试数据xtest，并根据拟合曲线求出ytest
ytest = model.predict(xtest[:,np.newaxis])   # model.predict → 预测

② 误差

ytest2 = model.predict(xtrain[:,np.newaxis])  # 样本数据x在拟合直线上的y值
plt.plot([xtrain,xtrain],[ytrain,ytest2],color = 'gray')  # 误差线

③ 求解a（斜率）、b（截距）

model = LinearRegression()   # LinearRegression → 线性回归评估器，用于拟合数据得到拟合直线
model.fit(xtrain[:,np.newaxis],ytrain)   # 回归拟合 # model.fit(x,y) → 拟合直线，参数分别为x与y   # x[:,np.newaxis] → 将数组变成(n,1)形状
print('斜率a为：%.4f' % model.coef_[0])
print('截距b为：%.4f' % model.intercept_)
print('线性回归函数为：\ny = %.4fx + %.4f' % (model.coef_[0],model.intercept_))   # 参数输出

④ 多元线性回归

df = pd.DataFrame(xtrain, columns = ['b1','b2','b3','b4'])
df['y'] = ytrain
model = LinearRegression()
model.fit(df[['b1','b2','b3','b4']],df['y'])   # 多元回归拟合
print('斜率a为：' ,model.coef_)
print('截距b为：%.4f' % model.intercept_)
print('线性回归函数为：\ny = %.1fx1 + %.1fx2 + %.1fx3 + %.1fx4 + %.1f' 
      % (model.coef_[0],model.coef_[1],model.coef_[2],model.coef_[3],model.intercept_))

3.线性回归模型评估
通过几个参数验证回归模型
SSE(和方差、误差平方和)：The sum of squares due to error
MSE(均方差、方差)：Mean squared error
RMSE(均方根、标准差)：Root mean squared error
R-square(确定系数) Coefficient of determination

ytest = model.predict(xtrain[:,np.newaxis])  # 求出预测数据
mse = metrics.mean_squared_error(ytrain,ytest)  # 求出均方差
rmse = np.sqrt(mse)  # 求出均方根

ssr = ((ytest - ytrain.mean())**2).sum()  # 求出预测数据与原始数据均值之差的平方和
sst = ((ytrain - ytrain.mean())**2).sum()  # 求出原始数据和均值之差的平方和
r2 = ssr / sst # 求出确定系数

r2 = model.score(xtrain[:,np.newaxis],ytrain)  # 求出确定系数
# 确定系数R-square非常接近于1，线性回归模型拟合较好

三、KNN最邻近分类

from sklearn import neighbors   # 导入KNN分类模块
import warnings
from sklearn import datasets
warnings.filterwarnings('ignore')    # 不发出警告
knn = neighbors.KNeighborsClassifier()   # 取得knn分类器
knn.fit(data[['fight','kiss']], data['type'])
print('预测电影类型为:', knn.predict([18, 90]))
# 加载数据，构建KNN分类模型
# 预测未知数据

四、PCA主成分分析

from sklearn.decomposition import PCA    # 加载主成分分析模块PCA # sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False)  
pca = PCA(n_components = 1)    # n_components = 1 → 降为1维 # n_components:  PCA算法中所要保留的主成分个数n，也即保留下来的特征个数n
pca.fit(df)     # 构建模型 # fit(X,y=None) → 调用fit方法的对象本身。比如pca.fit(X)，表示用X对pca这个对象进行训练
print(pca.explained_variance_)    # 输出特征值 # explained_variance_ratio_：返回 所保留的n个成分各自的方差百分比。
print(pca.components_)   # 输出特征向量 # components_：返回具有最大方差的成分
x_pca = pca.fit_transform(df)   # 数据转换 # fit_transform(X) → 用X来训练PCA模型，同时返回降维后的数据，这里x_pca就是降维后的数据
x_new = pca.inverse_transform(x_pca)   #inverse_transform() → 将降维后的数据转换成原始数据 # 主成分分析，生成新的向量x_pca

from sklearn.datasets import load_digits  # 多维数据降维
digits = load_digits()

五、K_means聚类

# 创建数据
from sklearn.datasets.samples_generator import make_blobs   # make_blobs聚类数据生成器
x,y_true = make_blobs(n_samples = 300,   # 生成300条数据，n_samples → 待生成的样本的总数。
                     centers = 4,        # 四类数据，centers → 类别数
                     cluster_std = 0.5,  # 方差一致，cluster_std → 每个类别的方差，如多类数据不同方差，可设置为[1.0,3.0]（这里针对2类数据）
                     random_state = 0)   # random_state → 随机数种子
# 构建K均值模型
from sklearn.cluster import KMeans   # 导入模块
kmeans = KMeans(n_clusters = 4)  # 这里为4簇
kmeans.fit(x)
y_kmeans = kmeans.predict(x)   # 构建模型，并预测出样本的类别y_kmeans
centroids = kmeans.cluster_centers_   # kmeans.cluster_centers_：得到不同簇的中心点

六、蒙塔卡罗模拟

① π的计算

from matplotlib.patches import Circle
# numpy.random.uniform(low,high,size) → 从一个均匀分布[low,high)中随机采样，均匀分布
d = np.sqrt((x-a)**2 + (y-b)**2)   # 计算点到圆心的距离
res = sum(np.where(d < r, 1, 0))   # 统计落在圆内的点的数目
pi = 4 * res / n   # 计算 pi （π）的近似值 → Monte Carlo方法：用统计值去近似真实值
plt.plot(x,y,'ro',markersize = 1)   #绘制正方形
circle = Circle(xy = (a,b),radius = r, alpha = 0.5 ,color = 'gray')
axes.add_patch(circle)   #绘制圆

② 计算积分 y = x**2

def f(x):
    return x**2   # 创建函数 y = x**2
res = sum(np.where(y < f(x), 1, 0))   # 统计 落在函数 y=x^2图像下方的点的数目
integral = res / n   # 计算定积分的近似值，落在下方的概率  # n:投点次数

③ 排队上厕所问题

startingtime[0] = arrivingtime[0]  # 第一个人之前没有人，所以开始时间 = 到达时间
finishtime[0] = startingtime[0] + workingtime[0]   # 第一个人完成时间 = 开始时间 + “工作”时间
waitingtime[0] = startingtime[0]-arrivingtime[0]   # 第一个人不用等待
for i in range(1,len(arrivingtime)):
    if finishtime[i-1] > arrivingtime[i]:
        startingtime[i] = finishtime[i-1]
    else:
        startingtime[i] = arrivingtime[i]
        emptytime[i] = arrivingtime[i] - finishtime[i-1]
    finishtime[i] = startingtime[i] + workingtime[i]
    waitingtime[i] = startingtime[i] - arrivingtime[i]
    print('第%d个人：到达时间 开始时间 “工作”时间 完成时间 等待时间\n' %i,
          arrivingtime[i],startingtime[i],workingtime[i],finishtime[i],waitingtime[i],
         '\n')
print('arerage waiting time is %f' %np.mean(waitingtime))