https://sklearn.apachecn.org/
1.离群值(点)识别
离群值(outlier)是指在一组数据中出现的与大部分数值相比差异较大的个别值。
关键问题：差异有多大、如何判断？
( 1)利用直方图或盒状图直接判断
置信水平， 显著偏离直方图主体频数区域的数值；或者距离盒状图的箱体(第 1 个四分位数和第 3 四分位数)过远，如超出3倍的IQR距离。
( 2)若一组数据符合某种分布( scipy.stats)，则可用置信区间或显著水平判断
#p-value 和 z-score 关系
P 值不是给定样本结果时原假设为真的概率，而是给定原假设为真时样本结果出现的概率。 计算出 P 值后，将给定的显著性水平 α 与 P 值比较，就可作出检验的结论:如果 α > P值，则在显著性水平 α 下拒绝原假设。如果 α ≤ P 值，则在显著性水平 α 下接受原假设。
# pdf - Probability density function
# cdf - Cumulative distribution function
# ppf - Percent point function (inverse of cdf — percentiles)
#z-score = ppf(p-value), p-value = cdf(z-score)
1 #(2)若一组数据符合某种分布(scipy.stats),则可用置信区间或显著水平判断
2 #3 importnumpy as np4 importmatplotlib.mlab as mlab5 importmatplotlib.pyplot as plt6 from scipy.stats importnorm,gengamma,gamma7
8 #模拟一组带噪音的数据
9 mu, sigma = 10, 0.25
10 x = np.random.normal(mu, sigma, 500) #samples
11 t = np.random.normal(8, 0.5, 10) #outliers
12 x[:10] =t13
14 #假设已知该组数据是正态(高斯)分布， 参数估计获得其参数
15 #mu_,sigma_ = rv_continuous.fit(norm, x)
16 mu_,sigma_ =norm.fit(x)17 rv =norm(mu_,sigma_)18
19 #Pr(c1<=μ<=c2)=1-α
20 #alpha=0.05; confidence level, 1-α; confidence interval [c1, c2]
21 #显著性水平 alpha=0.05，置信水平( 1-0.05=95%)， 置信区间[c1, c2]
22 #c1, c2= norm.interval(0.95, loc=mu_, scale=sigma_)
23 #c1, c2= rv.interval(0.95) #Percent point function (inverse of cdf — percentiles)
24
25 c1= rv.ppf(0.025) #置信极限值,左右两边的z值，小于9.2，大于10.69都是离群值
26 c2= rv.ppf(0.975)27 z_score =c128 p_value =rv.cdf(z_score) #Cumulative distribution function
29
30 #离群值
31 xmin,xmax =c1,c232 inx = np.where( (xxmax) > 0 )#这句代码是什么意思？
33 t_ =x[inx]34
35 #图示
36 n, bins, patches = plt.hist(x,50,density=1)37 y =norm.pdf(bins, mu, sigma)38 y_ =norm.pdf(bins, mu_, sigma_)39 plt.plot(bins, y, 'g--')40 plt.plot(bins, y_, 'r--')41 plt.plot([xmin,xmin],[0,1.5], 'r--')42 plt.plot([xmax,xmax],[0,1.5], 'r--')43 plt.show()
( 3) 给定一组样本数据符合一特定数学模型， 计算并画出其 confidence interval (bands)
非线性拟合方程
1 #( 3) 给定一组样本数据符合一特定数学模型， 计算并画出其 confidence interval (bands)
2 #3 from scipy.optimize importcurve_fit4 from scipy importstats5 importmatplotlib.pyplot as plt6 importnumpy as np7
8 #y = a*ebx + c
9 deffunc(x, a, b, c):10 '''Exponential 3-param function.'''
11 return a * np.exp(b * x) +c12
13 #Generating a dataset
14 n = 500
15 a,b,c = 50,0.2,1000
16 x = 11 + np.random.random(n)*12
17 #y= func(x,a,b,c)
18 y = a * np.exp(b * x) + c + (np.random.random(n)-0.2)*10*(x-11)**2 +\19 (np.random.random(n)-0.5)*10*(x-11)**2#加噪音
20
21 #Define confidence level
22 ci = 0.95
23 #Convert to percentile point of the normal distribution.
24 #See: https://en.wikipedia.org/wiki/Standard_score
25 pp = (1. + ci) / 2.26 #Convert to number of standard deviations.
27 nstd = stats.norm.ppf(pp)#对应n倍的lambda，一会要加在a,b,c的系数上
28
29 #Find best fit.
30 popt, pcov = curve_fit(func, x, y, p0=(20,1,1))31 #Standard deviation errors on the parameters.
32 perr = np.sqrt(np.diag(pcov))#把主对角线的拿过来开方
33
34 print(popt,pcov,perr,nstd)35 '''
36 [7.09214145e+01 1.88412538e-01 8.69468817e+02] 原来的是2037
38 [[ 1.32757358e+02 -7.76643168e-02 -8.36068437e+02] 协方差矩阵39 [-7.76643168e-02 4.55765204e-05 4.83351066e-01]40 [-8.36068437e+02 4.83351066e-01 5.64447127e+03]]41
42 [1.15220379e+01 6.75103847e-03 7.51296963e+01] 将主对角线开方的值43
44 1.95996398454005445 '''
46 #Add nstd standard deviations to parameters to obtain the upper confidence
47 #interval.
48 popt_up = popt + nstd *perr49 popt_dw = popt - nstd *perr50
51 #Plot data and best fit curve.
52 x_ = np.linspace(11, 23, 100)53 plt.scatter(x, y)54 plt.plot(x_, func(x_, *popt), c='g', lw=2.)#100个等间隔的x值放进去，估出来的a,b,c系数
55 plt.plot(x_, func(x_, *popt_up), '--', c='r', lw=2.)#等间隔的100个x
56 plt.plot(x_, func(x_, *popt_dw), '--', c='r', lw=2.)57 #x=12，y=8000位置上沿水平方向上写注记
58 plt.text(12, 8000, '{}% confidence level'.format(ci * 100.))59
60 plt.show()
(4) 点模式分析的离群值发现( MinCovDet - Minimum Covariance Determinant (MCD))
SCIKIT-LEARN中，非高斯分布可采用 OneClassSVM； 高斯分布数据可采用 MinCovDet和 EmpiricalCovariance，效果较好。其中， MinCovDet 鲁棒性好于 EmpiricalCovariance( 最大似然估计， MLE)。 二者都是用来估计出样本集的中心和形状参数，用于决定一个同样本集拟合程度最好的一个理论密度分布曲面。
EmpiricalCovariance( 经验协方差估计， 采用样本的最大似然估计法 MLE)，
， 就是要用似然函数取到最大值时的参数值作为估计值。 为何要使 L 最大呢，因为 X 已经发生了( 基于某一个参数发生的)， 通常它是属于那种最大概率的情况。通
常， EmpiricalCovariance 受噪声影响比较大。
MinCovDet( 最小协方差决定因子算法?)，适用于高斯分布数据，也可以用于单峰分布的数据(多峰不可以！)。 可以处理至多含有
离群值的样本集， 主要思想是寻找
个样本，使得协方差矩阵 determinant 行列式值最小。
确定离群值： 针对一个多维正态分布( 假设为高斯分布， 如： 2 维)，其分布曲面是一个密度曲面。一个样本点到该分布中心的马氏距离( Mahalanobis)， 其数值大小决定了该点
是否离群值。以下代码是用理论分布1生成黑点，用理论分布2生成红点(较少)，希望估算的时候基于黑点估算出它的分布，排除红点噪音。
1 importnumpy as np2 importmatplotlib.pyplot as plt3
4 from sklearn.covariance importEmpiricalCovariance, MinCovDet5
6 n_samples = 125
7 n_outliers = 25
8 n_features = 2
9
10 #generate data
11 gen_cov =np.eye(n_features)12 gen_cov[0, 0] = 2.#主对角线00位置上赋值为2
13 #生成了125个属性点，每个点有2个属性，后面又乘了2行2列的主矩阵
14 X =np.dot(np.random.randn(n_samples, n_features), gen_cov)15 #add some outliers
16 outliers_cov =np.eye(n_features)17 #生成了1个1，相当于将(1,1)位置上换成了7
18 outliers_cov[np.arange(1, n_features), np.arange(1, n_features)] = 7.19 #将倒数的25个样本给换了，从-多少开始，往下数，省略了一个冒号
20 X[-n_outliers:] =np.dot(np.random.randn(n_outliers, n_features), outliers_cov)21
22 #参数估计
23 #fit a Minimum Covariance Determinant (MCD) robust estimator to data
24 #MinCovDet()相当于调构造函数
25 robust_cov =MinCovDet().fit(X)26
27 #全样本得到n个马氏距离，给评价函数，迭代
28 #compare estimators learnt from the full data set with true parameters
29 emp_cov =EmpiricalCovariance().fit(X)30
31 #打一个q就退出来了
32
33 #结果
34 #Estimated robust location (n_features,), 均值
35 robust_cov.location_36 #array([-0.00320609, -0.07327767])
37
38
39 #Estimated robust covariance matrix (n_features, n_features)，协方差矩阵
40 robust_cov.covariance_41 #array([[2.25467783, 0.22281682],
42 #[0.22281682, 1.03483927]])
43
44
45 #规则网格 100x100, numpy.ravel - return a contiguous flattened array
46 xmin,ymin = list(np.min(X,axis=0))#打格子，样本值的第一列
47 xmax,ymax = list(np.max(X,axis=0))48 #xx.ravel().shape 维数1维，里面含有10000个数
49 xx, yy = np.meshgrid(np.linspace(xmin,xmax, 100), np.linspace(ymin,ymax, 100))50 xy = np.concatenate([xx.ravel()[:,np.newaxis], yy.ravel()[:,np.newaxis]], axis=1)51
52 #马氏距离表示点与一个分布之间的距离。 对于一个均值为 μ， 协方差矩阵为 Σ 的多变量
53 #向量，其马氏距离为 sqrt( (x-μ)'Σ^(-1)(x-μ) )
54 mahal_robust = robust_cov.mahalanobis(xy) #表示数据的协方差距离
55 mahal_robust =mahal_robust.reshape(xx.shape)56
57 #mahal_emp = emp_cov.mahalanobis(xy)
58 #mahal_emp = mahal_emp.reshape(xx.shape)
59
60 mahal_X =robust_cov.mahalanobis(X)61 X_ = X[np.where(mahal_X<=3)] #信号数据 mahal马氏距离，
62 X__ = X[np.where(mahal_X>3)] #噪音数据
63
64 ###############################################################################
65 #Display results
66 fig = plt.figure()#定义figure，相当于一张纸
67 #plt.scatter(X[:, 0], X[:, 1], color='black', label='inliers')
68 #从倒数第25个取，倒数125个值
69 plt.scatter(X[:, 0][-n_outliers:], X[:, 1][-n_outliers:], color='red', label='outliers')70 plt.scatter(X_[:, 0], X_[:, 1], color='black', label='selected')71 levels = np.arange(1, 16, 1)72 #等高线
73 robust_cs = plt.contour(xx, yy, np.sqrt(mahal_robust),levels=levels, cmap=plt.cm.YlOrBr_r,linestyles='dotted')74 plt.clabel(robust_cs, inline=1, fmt='%d', fontsize=10)75 #emp_cs = plt.contour(xx, yy, np.sqrt(mahal_emp), cmap=plt.cm.YlOrBr_r, linestyles='dotted')
76 #plt.clabel(emp_cs, inline=1, fmt='%d', fontsize=10)
77 plt.show()
(5)RANSAC (RANdom SAmple Consensus)
RANSAC 是一个迭代模型，通过从包含噪音观测样本中随机取样，估计纯样本(信号) 的理论分布模型。
主要分为两大步骤：
1)首选选择一个样本子集，数目必须满足一个基数，估算出模型参数。
2)然后检查整个数据集的样本是否和该模型一致，给定一个阈值或评估函数，偏离该模型的样本为离群值，符合者构成一个新样本集。
重复上述两个步骤，直到获得一个含有足够多样本的数据子集(有效数据)。其中， 步骤( 2)可能获得一个样本很少的子集，模型被拒绝；也可能获得一个样本多的子集，模型参数可以
被进一步细化。
RANSAC 中用于建模的分布模型，必须具有 fit(X, y)， score(X, y)两个函数。 前者用样本拟合模型，后者返回样本数据精度，以便判断是否结束迭代。
例如： 拟合线性模型，剔除离群值
1 #RANSAC
2 #避开异常点做迭代回归，参数估计，正向评价
3
4 """
5 ===========================================6 Robust linear model estimation using RANSAC7 ===========================================8 """
9 importnumpy as np10 from matplotlib importpyplot as plt11 from sklearn importlinear_model, datasets12
13 n_samples = 1000
14 n_outliers = 50
15
16 #生成1000个样本，feature只有1个，噪音含10个
17 X, y, coef = datasets.make_regression(n_samples=n_samples, n_features=1, \18 n_informative=1, noise=10, \19 coef=True, random_state=0)20
21 #Add outlier data
22 np.random.seed(0)23 #生成50个样本，放在1000个样本前面
24 X[:n_outliers] = 3 + 0.5 * np.random.normal(size=(n_outliers, 1))25 y[:n_outliers] = -3 + 10 * np.random.normal(size=n_outliers)26
27 #Fit line using all data
28 model =linear_model.LinearRegression()29 model.fit(X, y)30
31 #Robustly fit linear model with RANSAC algorithm
32 #RANSAC迭代的每一个步骤还是线性回归
33 #估算还是靠的linear_model
34 #inlier,outlier信号，噪音
35 #mask，把true的取出来
36 model_ransac =linear_model.RANSACRegressor(base_estimator=linear_model.LinearRegression())37 model_ransac.fit(X, y)38 inlier_mask =model_ransac.inlier_mask_39 outlier_mask = np.logical_not(inlier_mask)#取反，就是噪音
40
41 #Predict data of estimated models
42 line_X = np.arange(-5, 5)#生成了10个x值，等间隔为1
43 line_y = model.predict(line_X[:, np.newaxis])#预测出它的y值出来，画出黑色的直线
44 line_y_ransac =model_ransac.predict(line_X[:, np.newaxis])45 print ('y = ax +b, a=? b=?')46 #将RANSAC估计的a,b值打出来
47 print(model_ransac.estimator_.coef_, model_ransac.estimator_.intercept_)48
49 #Compare estimated coefficients
50 print ("Estimated coefficients (true, normal, RANSAC):")51 print(coef, model.coef_, model_ransac.estimator_.coef_)52
53 plt.plot(X[inlier_mask], y[inlier_mask], '.g', label='Inliers')54 plt.plot(X[outlier_mask], y[outlier_mask], '.r', label='Outliers')55 plt.plot(line_X, line_y, '-k', label='Linear regressor')56 plt.plot(line_X, line_y_ransac, '-b', label='RANSAC regressor')57 plt.legend(loc='lower right')58 plt.show()
2. 分类问题( 一个遥感图像分类的例子 – 贝叶斯分类器，决策树，知识向量机等)
分类是机器学习的一大类问题。解决分类需要理解( 1) 数据预处理、( 2) 特征提取、( 3)分类方法、以及( 4) 结果评估等关键环节。
#数据预处理
数据的标准化(Standardization)/归一化(Normalization)， 具体定义？
将数据按比例缩放，使之落入一个小的特定区间。在某些比较和评价的指标处理中经常会用到，去除数据的单位限制，将其转化为无量纲的纯数值，便于不同单位或量级的指标能够进行比较和加权。 主要目的:
1 把数变为特定区间的数值， 如： [0,1]内的小数为了数据处理方便提出来的，把数据映射到特定范围之内， 提高处理便捷性。
2 把有量纲表达式变为无量纲表达式
简化计算方式， 将有量纲的表达式， 变换为无量纲的表达式，成为纯量。
常见方法:
(1) min-max 标准化(Min-max standardization) --->[0,1]
(x-min)/(max-min)
(2) z-score 标准化(zero-mean normalization),
(x-mean)/std
(3) log 函数转换
(4)atan 函数转换
用反正切函数也可以实现数据的归一化。
… …
#通常，机器学习算法中的梯度方法对于数据的缩放和尺度都是很敏感的， scikit-learn 预处理使得特征数据缩放到 0-1 范围中
https://www.cnblogs.com/pejsidney/p/8031250.html
实际上下面这段代码几乎都不懂
1 #通常，机器学习算法中的梯度方法对于数据的缩放和尺度都是很敏感的， scikit-learn 预处
2 理使得特征数据缩放到 0-1范围中3 #https://www.cnblogs.com/pejsidney/p/8031250.htm
4 #5 from sklearn importpreprocessing6
7 #scale, (X - X_mean)/X_std #standard scaler对数据做归一化的公式
8 x = np.arange(9).reshape(3,3)#原始数据，3行3列
9 scaler =preprocessing.StandardScaler()10 scaler.fit(x) #算scaler归一化的参数
11 x_scale = scaler.transform(x)#将x集合里的每一个样本减均值除以标准差
12 np.min(x_scale), np.max(x_scale), np.mean(x_scale), np.std(x_scale)13 #14 x15
16 #scale, X_std=(X - min)/(max-min)做拉伸
17 x = np.arange(9).reshape(3,3)18 scaler =preprocessing.MinMaxScaler()19 scaler.fit(x)20 x_scale =scaler.transform(x)21 np.min(x_scale), np.max(x_scale), np.mean(x_scale), np.std(x_scale)22
23
24 #binarize二值化，给阈值，0,1
25 x = np.arange(9).reshape(3,3)26 binarizer = preprocessing.Binarizer(threshold= 4)27 x_bin =binarizer.transform(x)28 np.min(x_bin), np.max(x_bin), np.mean(x_bin), np.std(x_bin)29 #30 x_bin31
32 #normalize #第一个是按列来算的，后面的是按所有的来算的
33 x = np.arange(9).reshape(3,3)34 norm =preprocessing.Normalizer()35 norm.fit(x)36 x_norm =norm.transform(x)37 np.min(x_norm), np.max(x_norm), np.mean(x_norm), np.std(x_norm)38
39 #遥感图像读取与分类
40 #41 importnumpy as np42 importsys43 sys.path.append( r'E:\选课复习\研一上\Python\week7\chapter_06')44 from gwp_image importIMAGE45
46 #图像和训练区读取
47 drv =IMAGE()48 im_proj,im_geotrans,im_data = drv.read_img('taihu_gf2.tif')49 ao_proj,ao_geotrans,ao_data = drv.read_img('taihu_gf2_aoi.tif')50 band,row,col =im_data.shape51
52 im_data = np.swapaxes(np.swapaxes(im_data,0,1),1,2)#0轴1轴交换，维数变为2400,2000，第一次交换的结果1,2又做一次交换
53 print(im_data.shape, ao_data.shape)54 #(2000, 2000, 4) (2000, 2000)
55
56 #57 im_data[0,0]58 #array([442, 380, 215, 74], dtype=uint16)
59
60 #reshape into m_sample*n_feaure
61 im_data = im_data.reshape(row*col, band)62 ao_data = ao_data.reshape(row*col)63 print(im_data.shape, ao_data.shape)64 #(4000000,4)(4000000,)
65
66 im_data.shape67 #(4000000, 4)
68
69 #im_data[0,:]
70 #array([442, 380, 215, 74], dtype=uint16)#变换的结果要对
71
72 ao_data.shape73 #(4000000,)
74
75 np.uint(ao_data)76
77 #样本数据
78 from sklearn importpreprocessing79 from sklearn.model_selection importKFold80
81 #sampling data
82 x_sample = im_data[np.where(ao_data>0)]#ao_data像元属于哪一个类
83 y_sample = ao_data[np.where(ao_data>0)]84
85 samples = np.concatenate([x_sample, y_sample [:,np.newaxis]],axis=1)#沿着列合，沿着1轴
86 #samples.shape
87 #samples[0,:]
88 #像元的波段值，属于哪一类
89
90 #n_splits=3, 2/3 training, 1/3 verfying， 交叉验证 3 次
91 k_fold = KFold(n_splits=3, shuffle=False) #shuffle=True， 随机数种子启用， 则每次分组不同
92 cp=093 for train, test ink_fold.split(samples):94 #print('Train: %s | test: %s' % (train_indices, test_indices))
95 cp+=1
96 print(samples[ train].shape, samples[test].shape)97
98 #train set
99 x_train = samples [train][:,:-1]100 y_train = samples [train][:,-1]101
102 #test set
103 x_test = samples [test][:,:-1]104 y_test = samples [test][:,-1]105
106 #特征选择(Feature Selection)
107 '''
108 样本具有多个特征(或属性)，对于分类或回归等建模问题，这些特征起到了不同的作用，109 有的相关性强则贡献大些，而有的则比较弱。特征选择是困难的，经验和专业知识很重要，110 但是 scikt-learn 有现成的算法支持特征选择。111 '''
112 #树算法(Tree algorithms)是用于计算特征的信息量的一种常见方法
113 #114 from sklearn importmetrics115 from sklearn.ensemble importExtraTreesClassifier116 model =ExtraTreesClassifier()117
118 scaler =preprocessing.StandardScaler ()119 scaler.fit(x_train)120 x_train = scaler.transform(x_train)#原始数据0~255个值来训练结果
121 x_test =scaler.transform(x_test)122
123 model.fit(x_train, y_train)124 print(model.feature_importances_)125 #[0.30527885 0.32615108 0.18510445 0.18346563]
126
127 #朴素贝叶斯
128 #著名的机器学习算法，还原训练样本数据的分布密度，其在多类别分类中有很好的效果
129 #最大似然分类是朴素贝叶斯分类的特例，假设各类别的先验概率相等。
130 from sklearn importmetrics131 from sklearn.naive_bayes import GaussianNB#高斯分类器
132 model =GaussianNB()133
134 model.fit(x_train, y_train)135
136 y_pred = model.predict(x_train)#预测
137 len(np.where((y_pred-y_train)==0)[0])138 #33153
139 y_pred =model.predict(x_test)140 len(np.where((y_pred-y_test)==0)[0])#3365应该接近y_test.shape
141 #3365
142 #y_test.shape
143 #(16947,)实际上并不接近
144
145 #146 im_labe = model.predict( scaler.transform(im_data) )#把原始数据全部做一遍分类
147 #im_data[0,]
148 #array([442,380,215,74])看一看im_data是说明数据
149 im_labe =im_labe.reshape(row,col)150 print (np.unique(np.array(im_labe), return_counts=True) ) #??
151 im_labe = ((im_labe/10).astype(np.uint8))*10#将2000行2000列里的数据都除以10，按大类来，10,20,30
152
153 drv.write_img('taihu_gf2_lab01.tif',im_proj,im_geotrans,im_labe)#分类效果感觉并不好
154
155 #np.unique(im_labe)
156 #array([10,20,30,40])
157
158 #2/4=0.5 3.0以后算出小数来了，和过去不一样
#决策树
#分类与回归树(Classification and Regression Trees ,CART)算法常用于特征含有类别信息的分类或者回归问题，适用于多分类情况
#计算每类在一个属性上的信息熵，根据信息熵之间的差异，选定信息熵最大的属性最为最好的分类依据先进行分类。
唯一不同是构造函数不同，fit训练模型，predict用训练好的模型对整幅图像进行预测。
1 from sklearn importmetrics2 from sklearn.tree importDecisionTreeClassifier3 model =DecisionTreeClassifier()4 model.fit(x_train, y_train)5 y_pred =model.predict(x_train)6 len(np.where((y_pred-y_train)==0)[0])7 #33896
8 y_pred =model.predict(x_test)9 len(np.where((y_pred-y_test)==0)[0])10 #3388
11 #12 im_labe =model.predict( scaler.transform(im_data) )13 im_labe =im_labe.reshape(row,col)14 print (np.unique(np.array(im_labe), return_counts=True) ) #15 im_labe = (im_labe/10)*10
16
17 drv.write_img('taihu_gf2_lab02.tif',im_proj,im_geotrans,im_labe)
#支持向量机
#SVM 是非常流行的机器学习算法，主要用于分类问题，适用于一对多或多分类
1 #2 from sklearn importmetrics3 from sklearn.svm importSVC4
5 model.fit(x_train, y_train)6
7 y_pred =model.predict(x_train)8 len(np.where((y_pred-y_train)==0)[0])9 #33896
10 y_pred =model.predict(x_test)11 len(np.where((y_pred-y_test)==0)[0])12 #3388
13 #14 im_labe =model.predict( scaler.transform(im_data) )15 im_labe =im_labe.reshape(row,col)16 print (np.unique(np.array(im_labe), return_counts=True) ) #17 im_labe = (im_labe/10)*10
18
19 drv.write_img('taihu_gf2_lab03.tif',im_proj,im_geotrans,im_labe)
#结果评估或模型选择，针对三个模型依次执行，查看混淆矩阵的结果
1 print(model)2 '''
3 DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=None,4 max_features=None, max_leaf_nodes=None,5 min_impurity_decrease=0.0, min_impurity_split=None,6 min_samples_leaf=1, min_samples_split=2,7 min_weight_fraction_leaf=0.0, presort=False,8 random_state=None, splitter='best')9
10 '''
11 #make predictions
12 measured =y_test13 predicted =model.predict(x_test)14 #summarize the fit of the model
15 print(metrics.classification_report(measured, predicted))16 #非常低，只有0.2，...都是3000多
17 print(metrics.confusion_matrix(measured, predicted))18
19 #50000多个样本，Kfold分成两部分，给成false，关闭了采样的随机性，导致可能一些样本都没有
3.再论回归问题， 网格搜索与交叉验证的用途
给定一组样本(含因变量值和自变量值)，数学建模时模型参数的求解是一个关键问题。模型参数分两大类：一类是利用样本数据，通过数学模型反演和机器学习可以获得的参数，一类是凭借经验值或者统计与数据集特征相关的参数。
网格搜索(随机搜索) 与交叉验证主要用于在一个给定的参数分布空间内，搜索那些数学模型不能通过机器学习估算的参数，通过交叉验证，选出得分最高的参数作为模型最优输入。
例如：知识向量机的罚函数 C 与 gamma ( RBF)， Ridge 模型的 alpha 参数， Lasso 模型的 alpha参数等。
网格搜索 – 类似穷举法， 获得每个枚举参数组合的模型误差， 使得模型误差最小的参数对为优化解。 但是可以通过粗网格初选和细网格精选提高参数搜索效率；随机搜索 – 参数符合一个数学分布，多个参数搜索仍能保持较高搜索效率。
(1)以一个房屋价格回归模型为例。
https://archive.ics.uci.edu/ml/datasets.html
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names
1 #3.再论回归问题，网格搜索与交叉验证的用途
2
3 #数据读取
4 importnumpy as np5 fh= open("housing.data")6 #506行的大字符串
7 #对这一行先去了首尾的空白字符，包括空格，table，回车键\n,unix的\f
8 #成了506行，14列
9 raw_data = [[float(v) for v in ln.strip().split()] for ln infh.readlines()]10 fh.close()11 #强制转为32位，减少内存占用
12 dataset =np.array(raw_data).astype(np.float32)13 #要了13行以前的，第14行是价格
14 X = dataset[:,0:13]15 y = dataset[:,13]16
17 #多元回归+惩罚函数
18 importnumpy as np19 from sklearn.linear_model import Ridge#类名，惩罚函数的系数，只能是咱们给
20 #from sklearn.grid_search import GridSearchCV 由于版本问题，这一句不能执行
21 from sklearn.model_selection importGridSearchCV22 #prepare a range of alpha values to test
23 alphas = np.array([1,0.1,0.01,0.001,0.0001,0])24 alphas = np.array([0.01,0.001,0.0001,0])25 #create and fit a ridge regression model, testing each alpha
26 model = Ridge(fit_intercept=True, solver='svd')#ridge回归模型初始化，每个自变量前面加个待定系数，后边再加个Bata
27 #solver : {'auto', 'svd', 'cholesky', 'lsqr', 'sparse_cg', 'sag'}#多种办法帮助求出待定系数的w
28 #给定α的情况下求待定系数cv，用ridge模型做回归,α用刚才给的那个数，param_grid是写死的
29 #若是以简单参数必填的，直接给参数值，否则就用dict的格式
30 #出的结果都会放在grid里
31 grid = GridSearchCV(estimator=model, param_grid=dict(alpha=alphas))32 grid.fit(X, y)33 print(grid)34 '''
35 GridSearchCV(cv='warn', error_score='raise-deprecating',36 estimator=Ridge(alpha=1.0, copy_X=True, fit_intercept=True,37 max_iter=None, normalize=False, random_state=None,38 solver='svd', tol=0.001),39 iid='warn', n_jobs=None,40 param_grid={'alpha': array([0.01 , 0.001 , 0.0001, 0. ])},41 pre_dispatch='2*n_jobs', refit=True, return_train_score=False,42 scoring=None, verbose=0)43
44 '''
45
46 #summarize the results of the grid search
47 print(grid.best_score_)48 #-1.5507356822683775
49 print(grid.best_estimator_.alpha) #alpha
50 #0.01
51 print(grid.best_estimator_.coef_) #斜率
52 '''
53 [-1.0795517e-01 4.6437558e-02 2.0076692e-02 2.6850116e+0054 -1.7652218e+01 3.8107793e+00 5.9260055e-04 -1.4738853e+0055 3.0578226e-01 -1.2343894e-02 -9.5147777e-01 9.3177296e-0356 -5.2488232e-01]57
58 '''
59 print(grid.best_estimator_.intercept_) #截距
60 #36.378204
61
62 #给个新区间，多级网格搜索
63
64 #平均误差 STD，平均房价( 1000$)
65 #delta是预测价格和实际价格的差异
66 delta = np.dot(X,grid.best_estimator_.coef_) + grid.best_estimator_.intercept_ -y67 #求平方取均值开方
68 print(np.sqrt(np.mean(np.square(delta))), np.mean(y))69 #4.6791954 22.532806
70
71 #特征选取， SelectFromModel， #from version 0.17
72 #13个变量，想挑出其中特别重要的
73 from sklearn.feature_selection importSelectFromModel74 model = SelectFromModel(grid.best_estimator_, prefit=True)75 X_new =model.transform(X)76 X_new.shape77 #(506, 3)
78 X.shape79 #(506, 13)
80
81 X_new[2,:]82 #从13列里挑出了3列
83
84 #支持向量机回归(知识向量机?)
85
86 ##
87 #Grid-search and Cross-validated estimators
88 #These estimators are called similarly to their counterparts, with ‘CV’ appended to their name
89
90 from sklearn.model_selection importGridSearchCV91 from sklearn.model_selection importKFold, cross_val_score92 from sklearn.svm importSVR93 importnumpy as np94 importmatplotlib.pyplot as plt95
96 ###############################################################################
97 #Generate sample data
98 X = np.sort(5 * np.random.rand(40, 1), axis=0) #40个随机数，乘以5,0-5之间
99 y = np.sin(X).ravel()#求了sin值
100
101 ###############################################################################
102 #Add noise to targets
103 y[::5] += 3 * (0.5 - np.random.rand(8)) #间隔 5，取数据，0.5减去8个随机数，
104 #范围在-1.5 - +1.5原来样本里的8个数据，加了噪音
105 #40个数每隔5个数取出来，取了8个数出来
106 y107 y[::5]108
109 ###############################################################################
110 #Fit regression model
111 svr_rbf = SVR(kernel='rbf', C=1e3, gamma=0.1)112 svr_lin = SVR(kernel='linear', C=1e3)#线性函数作为支持向量机的核函数
113 svr_poly = SVR(kernel='poly', C=1e3, degree=2)114 y_rbf =svr_rbf.fit(X, y).predict(X)115 y_lin =svr_lin.fit(X, y).predict(X)116 y_poly =svr_poly.fit(X, y).predict(X)117 ###############################################################################
118 #look at the results
119 plt.scatter(X, y, c='k', label='data')120 plt.hold('on') #保持当前绘图和所有轴属性(包括当前颜色和线型)，以便随后绘图命令不重#
121 #上一句代码有点问题，会报错，无法执行
122 #AttributeError: module 'matplotlib.pyplot' has no attribute 'hold'
123 #置颜色和线型； off 返回默认模式， PLOT 命令借此擦出前面的绘图并重置新图所有坐标属
124 性。125 plt.plot(X, y_rbf, c='g', label='RBF model')126 plt.plot(X, y_lin, c='r', label='Linear model')127 plt.plot(X, y_poly, c='b', label='Polynomial model')128 plt.xlabel('data')129 plt.ylabel('target')130 plt.title('Support Vector Regression')131 plt.legend()132 plt.show()133 ###############################################################################
134 #grid search of svm for perfect parameters
135 C = np.array([2000,1000,500,250,125])136 g = np.array([1,0.5,0.25,0.125])137 svr_rbf = SVR(kernel='rbf')138 grid = GridSearchCV( estimator= svr_rbf, param_grid=dict(C=C,gamma=g) )#给区间网格搜索
139 grid.fit(X, y)140 print(grid)141 '''
142 GridSearchCV(cv='warn', error_score='raise-deprecating',143 estimator=SVR(C=1.0, cache_size=200, coef0=0.0, degree=3,144 epsilon=0.1, gamma='auto_deprecated', kernel='rbf',145 max_iter=-1, shrinking=True, tol=0.001,146 verbose=False),147 iid='warn', n_jobs=None,148 param_grid={'C': array([2000, 1000, 500, 250, 125]),149 'gamma': array([1. , 0.5 , 0.25 , 0.125])},150 pre_dispatch='2*n_jobs', refit=True, return_train_score=False,151 scoring=None, verbose=0)152
153 '''
154 print(grid.best_estimator_)155 '''
156 SVR(C=500, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma=0.25,157 kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False)158
159 '''
4.聚类分析
4.1 Scikt-learn支持的聚类方法
http://scikit-learn.org/stable/modules/clustering.html#clustering
常见聚类方法的 API 调用方法：
不同的聚类方法需要输入参数不同，如：类别数目、阈值、每类最少样本数等。有的自动化程度高些，有的需要用户不断测试参数直到获得一个较好质量的结果。
下面这部分代码能够跑一遍，但是对于具体的意义并没有理解。
1 #4 聚类分析
2 #4.1 Scikit-learn支持的聚类方法
3 #让稀疏地方的点向密集地方的点位移，第二圈再迭代，重新计算密度，所有点都向最大密度地方做位移
4 #无限循环下会形成几个聚类中心
5 #最后一点点位移到聚类中心
6 #迭代过程每一步样本点所有移的位置串成一根线，像河流形成聚类中心
7 #生成样本
8 #属性加距离形成权重
9 importnumpy as np10 from sklearn.datasets.samples_generator importmake_blobs11 #we create separable points
12 X, Y = make_blobs(n_samples=100, centers=2, random_state=0, cluster_std=0.60)13 X, Y = make_blobs(n_samples=1500, cluster_std=[1.0, 2.5, 0.5], random_state=170)#3个聚类中心
14 X.shape15 #(1500,2)
16 Y.shape17 #(1500,)
18 #Y代表在坐标上面会有属性值
19
20
21 #聚类 – MeanShift
22 cluster_name = "MeanShift"
23 from sklearn.cluster importMeanShift, estimate_bandwidth24 #ms = MeanShift(bandwidth=estimate_bandwidth(X, quantile=0.2, n_samples=1500), bin_seeding=True)
25 #带宽指针对某一列的值求标准差，相当于求波动范围
26 ms = MeanShift(bandwidth=estimate_bandwidth(X, n_samples=1500), bin_seeding=True)27 #fit函数放进去
28 ms.fit(X)29 #ms.labels_分了若干个类
30 labels =ms.labels_31 #np.unique(labels)
32 #array([0, 1, 2], dtype=int64),0是未分类的，1是一类，2是一类
33 #1500个点，分两簇
34 cluster_centers =ms.cluster_centers_35 cluster_centers.shape36 #(3,2)每个样本有两个属性
37 cluster_centers38 '''
39 array([[ 1.85568743, 0.44433028],40 [-8.92178147, -5.45806097],41 [-4.68205739, -0.28572227]])42
43 '''
44 labels.shape45 #(1500,)
46 X.shape47 #(1500,2)
48 np.where(labels==1)49 #打算取数据
50 X[np.where(labels==2),:].shape51 #(1,423,2)多了1，有错误,拿出了432行不是一个数组，而是数组外又套了个数组，应该要的是行的下标
52 X[np.where(labels==2)].shape53 #(423, 2)
54 X[np.where(labels==2)[0],:].shape#结果同上
55 Y.shape56 #(1500,)
57 np.unique(Y)58 #array([0,1,2])
59 Y[np.where(labels==2)[0]].shape60 #(423,)
61 #预测出的0有可能对应Y里的2
62
63 n_clusters_ =len(np.unique(labels))64 print("number of estimated clusters : %d" %n_clusters_)65 #number of estimated clusters : 3
66
67
68 #聚类 – DBSCAN
69 '''
70 ( 参数： 半径 eps， MinPts；关键词： 邻居、 核心点/边缘点/离群点、 密度(直接) 可达， 密度(直接)71 可达的若干核心对象及其邻域内的边缘点形成一个聚类簇)72 '''
73
74 cluster_name = "DBSCAN"
75 from sklearn.cluster importDBSCAN76 #eps : The maximum distance between two samples for them to be considered as in the sameneighborhood.
77 #min_samples : The number of samples (or total weight) in a neighborhood for a point to be
78 #considered as a core point. This includes the point itself.
79 #80 ms= DBSCAN (eps=0.75,min_samples=10)#eps找邻居找半径，看密度
81 ms.fit (X)82 '''
83 DBSCAN(algorithm='auto', eps=0.75, leaf_size=30, metric='euclidean',84 metric_params=None, min_samples=10, n_jobs=None, p=None)85
86 '''
87 labels =ms.labels_88 #cluster_centers = ms.cluster_centers_
89 n_clusters_ =len(np.unique(labels))90 print("number of estimated clusters : %d" %n_clusters_)91 #number of estimated clusters : 5
92
93 #聚类 - KMeans#M类给了M个初始化的中心，每一轮迭代都对聚类中心有一定的校正，当聚类中心
94 #不再偏移
95 cluster_name = "KMeans"
96 from sklearn.cluster importKMeans97 ms=KMeans(n_clusters=3, random_state=170)98 ms.fit(X)99 labels =ms.labels_100 cluster_centers =ms.cluster_centers_101 n_clusters_ =len(np.unique(labels))102 print("number of estimated clusters : %d" %n_clusters_)103 #number of estimated clusters : 3
104
105 #聚类 – AgglomerativeClustering (Hierarchical clustering 层次聚类)
106 #层次聚类能够不用给定分多少类，可以给定阈值自动决定类数
107 #n个样本，有n个类，每个样本聚类中心就是样本所在的位置
108 #任何两个点之间算距离，算最短距离，挑邻居合成一类，n类->n-1类
109 #新的聚类中心，两个加起来求平均
110 #3,4个点有12中连接方式
111 cluster_name = "AgglomerativeClustering"
112 from sklearn.cluster importAgglomerativeClustering113 #用于层次聚类的“ 距离” 定义， 用于判断是否合并
114 #complete - 元素对之间的最大距离
115 #single - 元素对之间的最小距离
116 #group average – 组间样本距离的平均值
117 #ward – “距离”指两簇 ESS 的和， 同合并后新簇 ESS 的差值( 增量)， ESS=Σ(xi-x_mean)**2
118 ms= AgglomerativeClustering (linkage="ward", n_clusters=3) #'ward', 'average', 'complete'
119 ms.fit(X)120 '''
121 AgglomerativeClustering(affinity='euclidean', compute_full_tree='auto',122 connectivity=None, distance_threshold=None,123 linkage='ward', memory=None, n_clusters=3,124 pooling_func='deprecated')125 '''
126 labels =ms.labels_127 #cluster_centers = ms.cluster_centers_
128 n_clusters_ =len(np.unique(labels))129 print("number of estimated clusters : %d" %n_clusters_)130 '''
131 number of estimated clusters : 3132 '''
133 #聚类 - spectral_clustering
134 ##谱聚类
135 #从图像分割来的
136 #将每一个像元看成一个点，每个点上有一组属性值
137 #任何一个像元都与周围像元连了8根线，任何一个点能到达任一个点
138 #连线需要计算权重，直线连接，对角线连接，波谱值相近程度，波谱值接近
139 #这两个点属于同一类，找出权限最弱的那一个打断，最后形成孤岛
140 #141 cluster_name = "spectral_clustering"
142 from sklearn.neighbors importkneighbors_graph143 from sklearn.cluster importspectral_clustering144
145 #没法标准化，计算速度慢，前面fit一下就可以，这里不行，它会把样本转为一个图
146 graph = kneighbors_graph(X, 25)147 labels= spectral_clustering(graph, n_clusters=3, eigen_solver='arpack')148 n_clusters_ =len(np.unique(labels))149 print("number of estimated clusters : %d" %n_clusters_)150 #number of estimated clusters : 3
151
152 #绘图显示聚类结果
153 importmatplotlib.pyplot as plt154 X0 = X[np.where(labels==0)]155 X1 = X[np.where(labels==1)]156 X2 = X[np.where(labels==2)]157 plt.plot(X[:,0],X[:,1],'ko')158 plt.plot(X0[:,0],X0[:,1],'ro')159 plt.plot(X1[:,0],X1[:,1],'bo')160 plt.plot(X2[:,0],X2[:,1],'go')161 plt.title(cluster_name)162 plt.show()163
164 #4.2非常非常慢，图像分割会分成非常少，50x50个，2000个像素的样本
5867
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
