半部秘籍–分类、回归、无监督与集成学习

作者: MoBaiGeneral 分类: 传统机器学习,分类,回归,无监督,集成学习 发布时间: 2021-10-20 00:48 访问量:13,558
FavoriteLoading收藏

原文

1 分类器

1.1 常用分类器

from sklearn.ensemble import VotingClassifier
from sklearn.linear_model import LogisticRegression
#pip install xgboost -i https://pypi.tuna.tsinghua.edu.cn/simple
from xgboost.sklearn import XGBClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.ensemble import  AdaBoostClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.svm import SVC
from sklearn.svm import LinearSVC
from sklearn import decomposition
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.ensemble import BaggingClassifier
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import f_regression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV
from catboost import CatBoostClassifier
from mlxtend.classifier import StackingClassifier
from sklearn.neural_network import MLPClassifier
from sklearn import cluster, mixture
from sklearn.neighbors import kneighbors_graph
from sklearn.multiclass import OneVsRestClassifier
from sklearn.naive_bayes import MultinomialNB
import lightgbm as lgb
classifiers={
    'MLP':MLPClassifier(),
    'RandomForest':RandomForestClassifier(),
    'DecisionTree':DecisionTreeClassifier(),
    'KNeighbors':KNeighborsClassifier(),
    'Logistic':LogisticRegression(),
    'AdaBoost':AdaBoostClassifier(),
    'GaussianNB':GaussianNB(),
    'LinearDis':LinearDiscriminantAnalysis(),
    'QuadraticDis':QuadraticDiscriminantAnalysis(),
    'SVM':SVC(),
    'ExtraTrees':ExtraTreesClassifier(),
    'Bagging':BaggingClassifier(DecisionTreeClassifier(6),
    'XGB':XGBClassifier(),
    'GradientBoosting':GradientBoostingClassifier(),
    'GradientBoosting':GradientBoostingClassifier(),
    'CatBoost':CatBoostClassifier(),
    'OneVsRest':OneVsRestClassifier(AdaBoostClassifier()),
    'Multion':MultinomialNB(),
    'BalBagg':BalancedBaggingClassifier(base_estimator=DecisionTreeClassifier())                           
}

1.2 分类器评价

分类任务常见的评价指标有准确率(Accuracy)、精确率(Precision)、召回率(Recall)、F1 score、ROC曲线(Receiver Operating Characteristic Curve)等。

把正例正确分类为正例,表示为TP(true positive),把正例错误分类为负例,表示为FN(false negative),

把负例正确分类为负例,表示为TN(true negative), 把负例错误分类为正例,表示为FP(false positive)

精确率和召回率可以从混淆矩阵中计算而来,precision = TP/(TP + FP), recall = TP/(TP +FN)。

准确率是分类正确的样本占总样本个数的比例

准确率是分类问题中最简单直观的评价指标,但存在明显的缺陷。比如如果样本中有99%的样本为正样本,那么分类器只需要一直预测为正,就可以得到99%的准确率,但其实际性能是非常低下的。也就是说,当不同类别样本的比例非常不均衡时,占比大的类别往往成为影响准确率的最主要因素。

在多分类问题中一般不直接使用整体的分类准确率,而是使用每个类别下的样本准确率的算术平均作为模型的评估指标。

1.2.1 准确率

准确率是分类问题中最简单直观的评价指标,但存在明显的缺陷。比如如果样本中有99%的样本为正样本,那么分类器只需要一直预测为正,就可以得到99%的准确率,但其实际性能是非常低下的。也就是说,当不同类别样本的比例非常不均衡时,占比大的类别往往成为影响准确率的最主要因素。

在多分类问题中一般不直接使用整体的分类准确率,而是使用每个类别下的样本准确率的算术平均作为模型的评估指标。

import numpy as np
from sklearn.metrics import accuracy_score
y_pred = [0, 2, 1, 3]
y_true = [0, 1, 2, 3]
print(accuracy_score(y_true, y_pred))  # 0.5
print(accuracy_score(y_true, y_pred, normalize=False))  # 2

1.2.2 精确率

精确率指模型预测为正的样本中实际也为正的样本占被预测为正的样本的比例。

Macro Average: 宏平均是指在计算均值时使每个类别具有相同的权重,最后结果是每个类别的指标的算术平均值。 Micro Average: 微平均是指计算多分类指标时赋予所有类别的每个样本相同的权重,将所有样本合在一起计算各个指标。

‘weighted’: 为每个标签计算指标,并通过各类占比找到它们的加权均值(每个标签的正例数).它解决了’macro’的标签不平衡问题;它可以产生不在精确率和召回率之间的F-score.

‘samples’: 为每个实例计算指标,找到它们的均值(只在多标签分类的时候有意义,并且和函数accuracy_score不同).

当average参数为None时,得到的结果是每个类别的precision。

from sklearn.metrics import precision_score
y_true = [0, 1, 2, 0, 1, 2]
y_pred = [0, 2, 1, 0, 0, 1]
print(precision_score(y_true, y_pred, average='macro')) 
print(precision_score(y_true, y_pred, average='micro')) 
print(precision_score(y_true, y_pred, average='weighted'))
print(precision_score(y_true, y_pred, average=None))

1.2.3 召回率

召回率指实际为正的样本中被预测为正的样本所占实际为正的样本的比例。 sklearn中recall_score方法和precision_score方法的参数说明都是一样的

Recall和Precision只有计算公式不同,它们average参数为’macro’,‘micro’,’weighted’和None时的计算方式都是相同的,

from sklearn.metrics import recall_score
y_true = [0, 1, 2, 0, 1, 2]
y_pred = [0, 2, 1, 0, 0, 1]
print(recall_score(y_true, y_pred, average='macro')) 
print(recall_score(y_true, y_pred, average='micro')) 
print(recall_score(y_true, y_pred, average='weighted')) 
print(recall_score(y_true, y_pred, average=None)) 

1.2.4 混淆矩阵

横为true label 竖为predict

y_pred = [0, 2, 1, 3,9,9,8,5,8]
y_true = [0, 1, 2, 3,2,6,3,5,9]
# 混淆矩阵
from sklearn.metrics import confusion_matrix
confusion_matrix(y_true, y_pred)

例:

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets
from sklearn.model_selection import train_test_split
from sklearn.metrics import plot_confusion_matrix
# import some data to play with
iris = datasets.load_iris()
X = iris.data
y = iris.target
class_names = iris.target_names
# Split the data into a training set and a test set
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
# Run classifier, using a model that is too regularized (C too low) to see
# the impact on the results
classifier = svm.SVC(kernel='linear', C=0.01).fit(X_train, y_train)
np.set_printoptions(precision=2)
# Plot non-normalized confusion matrix
titles_options = [("Confusion matrix, without normalization", None),
                  ("Normalized confusion matrix", 'true')]
for title, normalize in titles_options:
    disp = plot_confusion_matrix(classifier, X_test, y_test,
                                 display_labels=class_names,
                                 cmap=plt.cm.Blues,
                                 normalize=normalize)
    disp.ax_.set_title(title)
    print(title)
    print(disp.confusion_matrix)
plt.show()

1.2.5 分类报告

包含:precision/recall/fi-score/均值/分类个数

# 分类报告:precision/recall/fi-score/均值/分类个数
from sklearn.metrics import classification_report
y_true = [0, 2, 2, 2, 0,1,4,3,4]
y_pred = [0, 0, 2, 2, 0,2,3,4,2]
target_names = ['class 0', 'class 1', 'class 2', 'class 5', 'class 6']
print(classification_report(y_true, y_pred, target_names=target_names))

1.2.6 F1

F1 score是精确率和召回率的调和平均值,F1 score越高,说明模型越稳健。

F1 score的最好值为1,最差值为0. 精确率和召回率对F1 score的相对贡献是相等的. F1 score的计算公式为: F1 = 2 * (precision * recall) / (precision + recall)

from sklearn.metrics import f1_score
y_true = [0, 1, 2, 0, 1, 2]
y_pred = [0, 2, 1, 0, 0, 1]
print(f1_score(y_true, y_pred, average='macro'))  # 0.26666666666666666
print(f1_score(y_true, y_pred, average='micro'))  # 0.3333333333333333
print(f1_score(y_true, y_pred, average='weighted'))  # 0.26666666666666666
print(f1_score(y_true, y_pred, average=None))  # [0.8 0.  0. ]

1.2.7 kappa score

kappa score是一个介于(-1, 1)之间的数. score>0.8意味着好的分类;0或更低意味着不好(实际是随机标签)

from sklearn.metrics import cohen_kappa_score
y_true = [2, 0, 2, 2, 0, 1]
y_pred = [0, 0, 2, 2, 0, 2]
cohen_kappa_score(y_true, y_pred)

1.2.8 海明距离

from sklearn.metrics import hamming_loss
y_pred = [1, 2, 3, 4]
y_true = [2, 2, 3, 4]
hamming_loss(y_true, y_pred)

1.2.9 Jaccard距离

import numpy as np
from sklearn.metrics import jaccard_score
y_pred = [0, 2, 1, 3]
y_true = [0, 1, 2, 3]
jaccard_score(y_true, y_pred,average='macro')

1.2.10 铰链损失

铰链损失(Hinge loss)一般用来使“边缘最大化”(maximal margin)。损失取值在0~1之间,当取值为0,表示多分类模型分类完全准确,取值为1表明完全不起作用。

from sklearn.metrics import hinge_loss
y_pred = [0, 1,0,1,1]
y_true = [0, 1,1,0,0]
hinger = hinge_loss(y_true,y_pred)
print(hinger)

1.2.11 P-R曲线

评价一个模型的好坏,不能仅靠精确率或者召回率,最好构建多组精确率和召回率,绘制出模型的P-R曲线。 下面说一下P-R曲线的绘制方法。P-R曲线的横轴是召回率,纵轴是精确率。P-R曲线上的一个点代表着,在某一阈值下,模型将大于该阈值的结果判定为正样本,小于该阈值的结果判定为负样本,此时返回结果对应的召回率和精确率。整条P-R曲线是通过将阈值从高到低移动而生成的。原点附近代表当阈值最大时模型的精确率和召回率。

import numpy as np
from sklearn.metrics import precision_recall_curve
y_true = np.array([0, 0, 1, 1])
y_scores = np.array([0.1, 0.4, 0.35, 0.8])
precision, recall, thresholds = precision_recall_curve( y_true, y_scores)
print(precision,'\n', recall,'\n', thresholds)

例一:二分类

#二分类曲线
from sklearn.datasets import make_classification
from sklearn.metrics import (precision_recall_curve,PrecisionRecallDisplay)
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
from sklearn.metrics import average_precision_score
X, y = make_classification(random_state=0)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
clf = SVC(random_state=0)
clf.fit(X_train, y_train)
SVC(random_state=0)
predictions = clf.predict(X_test)
precision, recall, _ = precision_recall_curve(y_test, predictions)
disp = PrecisionRecallDisplay(precision=precision, recall=recall)
disp.plot() 
average_precision = average_precision_score(y_test, predictions)
print('Average precision-recall score: {0:0.2f}'.format(
      average_precision))

例二:多分类

def plot_precision_recall_curve(Y_test,y_score,n_classes):
    # For each class
    precision = dict()
    recall = dict()
    average_precision = dict()
    print("be sure of y_true.shape=(,n_classes),y_pred.shape=(,n_classes)")
    for i in range(n_classes):
        precision[i], recall[i], _ = precision_recall_curve(Y_test[:, i],
                                                            y_score[:, i])
        average_precision[i] = average_precision_score(Y_test[:, i], y_score[:, i])
       
    # A "micro-average": quantifying score on all classes jointly
    precision["micro"], recall["micro"], _ = precision_recall_curve(Y_test.ravel(),
        y_score.ravel())
    average_precision["micro"] = average_precision_score(Y_test, y_score,
                                                         average="micro")
    print('Average precision score, micro-averaged over all classes: {0:0.2f}'
          .format(average_precision["micro"]))
    
    #Plot the micro-averaged Precision-Recall curve
    plt.figure()
    plt.step(recall['micro'], precision['micro'], where='post')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.ylim([0.0, 1.05])
    plt.xlim([0.0, 1.0])
    plt.title(
        'Average precision score, micro-averaged over all classes: AP={0:0.2f}'
        .format(average_precision["micro"]))
    
    #Plot Precision-Recall curve for each class and iso-f1 curves
    
    plt.figure(figsize=(7, 8))
    f_scores = np.linspace(0.2, 0.8, num=4)
    lines = []
    labels = []
    for f_score in f_scores:
        x = np.linspace(0.01, 1)
        y = f_score * x / (2 * x - f_score)
        l, = plt.plot(x[y >= 0], y[y >= 0], color='gray', alpha=0.2)
        plt.annotate('f1={0:0.1f}'.format(f_score), xy=(0.9, y[45] + 0.02))

    lines.append(l)
    labels.append('iso-f1 curves')
    l, = plt.plot(recall["micro"], precision["micro"], color='gold', lw=2)
    lines.append(l)
    labels.append('micro-average Precision-recall (area = {0:0.2f})'
                  ''.format(average_precision["micro"]))
    
    '''
    # setup plot details
    colors = cycle(['navy', 'turquoise', 'darkorange', 'cornflowerblue', 'teal'])
    for i, color in zip(range(n_classes), colors):
        l, = plt.plot(recall[i], precision[i], color=color, lw=2)
        lines.append(l)
        labels.append('Precision-recall for class {0} (area = {1:0.2f})'
                      ''.format(i, average_precision[i]))
    '''
    for i in range(n_classes):
        l, = plt.plot(recall[i], precision[i],lw=2)
        lines.append(l)
        labels.append('Precision-recall for class {0} (area = {1:0.2f})'
                      ''.format(i, average_precision[i]))
        
    fig = plt.gcf()
    fig.subplots_adjust(bottom=0.25)
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.title('Extension of Precision-Recall curve to multi-class')
    plt.legend(lines, labels, loc=(0, -.38), prop=dict(size=14))
    plt.show()

在这里插入图片描述

1.2.12 ROC曲线和AUC

例一:二分类

#!/usr/bin/python
# -*- coding:utf-8 -*-
 
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegressionCV
from sklearn import metrics
from sklearn.preprocessing import label_binarize
from sklearn import datasets
if __name__ == '__main__':
    np.random.seed(0)
    
    iris = datasets.load_iris()
    x = iris.data
    y = iris.target
    print(x.shape,y.shape)
    print(y[:10])
    print(set(y))
    
    print(y.shape)
    print(y[:10])
    
    x_train, x_test, y_train, y_test = train_test_split(x, y, train_size = 0.6, random_state = 0)
    Y = label_binarize(y_test, classes=[0, 1, 2])
    n_classes = Y.shape[1]
    alpha = np.logspace(-2, 2, 20)  #设置超参数范围
    model = LogisticRegressionCV(Cs = alpha, cv = 3, penalty = 'l2')  #使用L2正则化
    model.fit(x_train, y_train)
    print('超参数:', model.C_)
    # 计算属于各个类别的概率,返回值的shape = [n_samples, n_classes]
    y_score = model.predict_proba(x_test)
    # 1、调用函数计算micro类型的AUC
    print('调用函数auc:', metrics.roc_auc_score(Y, y_score, average='micro'))
    # 2、手动计算micro类型的AUC
    #首先将矩阵y_one_hot和y_score展开,然后计算假正例率FPR和真正例率TPR
    fpr, tpr, thresholds = metrics.roc_curve(Y.ravel(),y_score.ravel())
    auc = metrics.auc(fpr, tpr)
    print('手动计算auc:', auc)
    #绘图
    mpl.rcParams['font.sans-serif'] = u'SimHei'
    mpl.rcParams['axes.unicode_minus'] = False
    #FPR就是横坐标,TPR就是纵坐标
    plt.plot(fpr, tpr, c = 'r', lw = 2, alpha = 0.7, label = u'AUC=%.3f' % auc)
    plt.plot((0, 1), (0, 1), c = '#808080', lw = 1, ls = '--', alpha = 0.7)
    plt.xlim((-0.01, 1.02))
    plt.ylim((-0.01, 1.02))
    plt.xticks(np.arange(0, 1.1, 0.1))
    plt.yticks(np.arange(0, 1.1, 0.1))
    plt.xlabel('False Positive Rate', fontsize=13)
    plt.ylabel('True Positive Rate', fontsize=13)
    plt.grid(b=True, ls=':')
    plt.legend(loc='lower right', fancybox=True, framealpha=0.8, fontsize=12)
    plt.title(u'鸢尾花数据Logistic分类后的ROC和AUC', fontsize=17)
    plt.show()

例二:多分类

def plot_roc_curve(y_test,y_score,n_classes,y_prob=None):
    # Compute ROC curve and ROC area for each class
    fpr = dict()
    tpr = dict()
    roc_auc = dict()
    for i in range(n_classes):
        fpr[i], tpr[i], _ = roc_curve(y_test[:, i], y_score[:, i])
        roc_auc[i] = auc(fpr[i], tpr[i])

    # Compute micro-average ROC curve and ROC area
    fpr["micro"], tpr["micro"], _ = roc_curve(y_test.ravel(), y_score.ravel())
    roc_auc["micro"] = auc(fpr["micro"], tpr["micro"])
    lw = 2
    '''
    #Plot of a ROC curve for a specific class
    plt.figure()
    lw = 2
    plt.plot(fpr[2], tpr[2], color='darkorange',
             lw=lw, label='ROC curve (area = %0.2f)' % roc_auc[2])
    plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Receiver operating characteristic example')
    plt.legend(loc="lower right")
    plt.show()
    '''
    
    #Plot ROC curves for the multilabel problem
    # First aggregate all false positive rates
    all_fpr = np.unique(np.concatenate([fpr[i] for i in range(n_classes)]))

    # Then interpolate all ROC curves at this points
    mean_tpr = np.zeros_like(all_fpr)
    for i in range(n_classes):
        mean_tpr += interp(all_fpr, fpr[i], tpr[i])

    # Finally average it and compute AUC
    mean_tpr /= n_classes

    fpr["macro"] = all_fpr
    tpr["macro"] = mean_tpr
    roc_auc["macro"] = auc(fpr["macro"], tpr["macro"])

    # Plot all ROC curves
    plt.figure()
    plt.plot(fpr["micro"], tpr["micro"],
             label='micro-average ROC curve (area = {0:0.2f})'
                   ''.format(roc_auc["micro"]),
             color='deeppink', linestyle=':', linewidth=4)

    plt.plot(fpr["macro"], tpr["macro"],
             label='macro-average ROC curve (area = {0:0.2f})'
                   ''.format(roc_auc["macro"]),
             color='navy', linestyle=':', linewidth=4)

    colors = cycle(['aqua', 'darkorange', 'cornflowerblue'])
    for i in range(n_classes):
        plt.plot(fpr[i], tpr[i], lw=2,
                 label='ROC curve of class {0} (area = {1:0.2f})'
                 ''.format(i, roc_auc[i]))

    plt.plot([0, 1], [0, 1], 'k--', lw=lw)
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Some extension of Receiver operating characteristic to multi-class')
    plt.legend(loc="lower right")
    plt.show()
    
    #Area under ROC for the multiclass problem
    #y_prob = classifier.predict_proba(X_test)
    if np.all(y_prob):
        macro_roc_auc_ovo = roc_auc_score(y_test, y_prob, multi_class="ovo",
                                          average="macro")
        weighted_roc_auc_ovo = roc_auc_score(y_test, y_prob, multi_class="ovo",
                                             average="weighted")
        macro_roc_auc_ovr = roc_auc_score(y_test, y_prob, multi_class="ovr",
                                          average="macro")
        weighted_roc_auc_ovr = roc_auc_score(y_test, y_prob, multi_class="ovr",
                                             average="weighted")
        print("One-vs-One ROC AUC scores:\n{:.6f} (macro),\n{:.6f} "
              "(weighted by prevalence)"
              .format(macro_roc_auc_ovo, weighted_roc_auc_ovo))
        print("One-vs-Rest ROC AUC scores:\n{:.6f} (macro),\n{:.6f} "
              "(weighted by prevalence)"
              .format(macro_roc_auc_ovr, weighted_roc_auc_ovr))

在这里插入图片描述

2 回归器

2.1 常用回归器

from sklearn.ensemble import VotingRegressor
#pip install xgboost -i https://pypi.tuna.tsinghua.edu.cn/simple
from xgboost.sklearn import XGBRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble import  AdaBoostRegressor
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.ensemble import BaggingRegressor
from catboost import CatBoostRegressor
from mlxtend.regressor import StackingRegressor
from sklearn.neural_network import MLPRegressor
from GCForest import GCForest
from sklearn.linear_model import ARDRegression#ARD贝叶斯ARD回归
from sklearn.linear_model import BayesianRidge#BayesianRidge贝叶斯岭回归
from sklearn.linear_model import TheilSenRegressor#TheilSen泰尔森估算
from sklearn.linear_model import RANSACRegressor#RANSAC随机抽样一致性算法
from sklearn.svm import SVR
from sklearn.linear_model import LinearRegression, LassoCV
from sklearn.model_selection import KFold
from sklearn.preprocessing import PowerTransformer
from sklearn.preprocessing import PolynomialFeatures
regressors={
    'RandomForest':RandomForestRegressor(n_estimators=130,random_state=3),
    'KNeighbors':KNeighborsRegressor(n_neighbors=4),
    'SVM':SVR(),
    'ARD':ARDRegression(),
    'BayesianRidge':BayesianRidge(),
    'TheilSen':TheilSenRegressor(),
    'RANSAC':RANSACRegressor(),
    'ExtraTrees':ExtraTreesRegressor(),
    'DecisionTree':DecisionTreeRegressor(random_state=3),
    'GradientBoosting':GradientBoostingRegressor(),
    'Bagging':BaggingRegressor(),
    'XGB':XGBRegressor(),
    'CatBoost':CatBoostRegressor(verbose=0),
    'AdaBoost':AdaBoostRegressor( base_estimator=BayesianRidge()),
    'MLP':MLPRegressor(hidden_layer_sizes=(110, 110, 110, 110, 110)),
    'Linear':LinearRegression()
}

2.2 回归器评价

1.可释方差分数(explain variance score)

2.平均绝对误差(mean absolute error)

3.均方误差(mean squared error)

4.均方对数误差(mean squared logarithmic error),适用于具有指数增长的趋势的目标。

5.中值绝对误差(median absolute error),该函数不支持多输出。

6.决定系数分数( score)

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.metrics import median_absolute_error
from sklearn.metrics import mean_squared_log_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import explained_variance_score
from sklearn.metrics import r2_score
import numpy as np

bos_house = datasets.load_boston()
bos_house_data = bos_house['data']
bos_house_target = bos_house['target']

x_train,x_test,y_train,y_test = train_test_split(bos_house_data,bos_house_target,random_state=41)
forest_reg = RandomForestRegressor(random_state=41)
forest_reg.fit(x_train,y_train)
y_pred = forest_reg.predict(x_test)

#mean_squared_error
print('MSE为:',mean_squared_error(y_test,y_pred))
print('MSE为(直接计算):',np.mean((y_test-y_pred)**2))

print('RMSE为:',np.sqrt(mean_squared_error(y_test,y_pred)))

#median_absolute_error
print(np.median(np.abs(y_test-y_pred)))
print(median_absolute_error(y_test,y_pred))

#mean_absolute_error
print(np.mean(np.abs(y_test-y_pred)))
print(mean_absolute_error(y_test,y_pred))

#mean_squared_log_error
print(mean_squared_log_error(y_test,y_pred))
print(np.mean((np.log(y_test+1)-np.log(y_pred+1))**2))

#explained_variance_score
print(explained_variance_score(y_test,y_pred))
print(1-np.var(y_test-y_pred)/np.var(y_test))

#r2_score
print(r2_score(y_test,y_pred))
print(1-(np.sum((y_test-y_pred)**2))/np.sum((y_test -np.mean(y_test))**2))

3 无监督学习

无监督学习主要算法

1 基于划分聚类算法(partition clustering)

优点:计算时间短,速度快。结果容易解释,一般聚类效果还算不错;

缺点:对异常值非常敏感,需要提前确定好k值

2 基于层次聚类算法

优点:采用随机抽样与分割相结合的办法来提高算法的空间和时间效率,并且在算法中用了堆和K-d树结构来提高了算法效率,使其可以高效的处理大量数据。

缺点:对异常数据比较脆弱。

CURE:采用抽样技术先对数据集D随机抽取样本,再采用分区技术对样本进行分区,然后对每个分区局部聚类,最后对局部聚类进行全局聚类。

3 基于密度聚类算法

优点:聚类簇的形状没有偏倚,不需要输入要划分的聚类个数。

缺点:DBSCAN算法对参数Eps及Minpts非常敏感,且这两个参数很难确定。

DBSCAN:DBSCAN算法是一种典型的基于密度的聚类算法,该算法采用空间索引技术来搜索对象的邻域,引入了“核心对象”和“密度可达”等概念,从核心对象出发,把所有密度可达的对象组成一个簇。

3.1 例:

这个例子演示不同的聚类算法用于6个标准数据集(toy datasets)聚类上的特性。这6个数据集分别是:

noisy_circles: 数据点围成大小两个同心圈状。

noisy_moons: 数据点构成两个交错的半圆。

blobs: 数据点形如团状高斯块。

varied: 可变方差的数据块。

aniso: 各向异性分布的数据块。

no_structure: 均匀分布的数据。

使用的聚类算法是:

MiniBatchKMeans,AP聚类,MeanShift,谱聚类,Ward聚类,层次聚类,DBSCAN聚类,Birch聚类,高斯混合模型聚类

除了最后一个数据集no_structure, 每一对数据集-算法的参数都被优化到产生最佳的聚类结果。我们发现,一些算法对参数是敏感的。no_structure数据集是一个null情形的例子。数据是同分布的,所以没有好的聚类结果。

import time
import warnings
import numpy as np
import matplotlib.pyplot as plt
from sklearn import cluster, datasets, mixture
from sklearn.neighbors import kneighbors_graph
from sklearn.preprocessing import StandardScaler
from itertools import cycle, islice
from sklearn.metrics import silhouette_score
np.random.seed(0)

# ============
# Generate datasets. We choose the size big enough to see the scalability
# of the algorithms, but not too big to avoid too long running times
# ============
n_samples = 1500
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5,
                                      noise=.05)
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)
no_structure = np.random.rand(n_samples, 2), None

# Anisotropicly distributed data
random_state = 170
X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
transformation = [[0.6, -0.6], [-0.4, 0.8]]
X_aniso = np.dot(X, transformation)
aniso = (X_aniso, y)

# blobs with varied variances
varied = datasets.make_blobs(n_samples=n_samples,
                             cluster_std=[1.0, 2.5, 0.5],
                             random_state=random_state)

# ============
# Set up cluster parameters
# ============
plt.figure(figsize=(9 * 2 + 3, 12.5))
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05,
                    hspace=.01)

plot_num = 1

default_base = {'quantile': .3,
                'eps': .3,
                'damping': .9,
                'preference': -200,
                'n_neighbors': 10,
                'n_clusters': 3,
                'min_samples': 20,
                'xi': 0.05,
                'min_cluster_size': 0.1}

datasets = [
    (noisy_circles, {'damping': .77, 'preference': -240,
                     'quantile': .2, 'n_clusters': 2,
                     'min_samples': 20, 'xi': 0.25}),
    (noisy_moons, {'damping': .75, 'preference': -220, 'n_clusters': 2}),
    (varied, {'eps': .18, 'n_neighbors': 2,
              'min_samples': 5, 'xi': 0.035, 'min_cluster_size': .2}),
    (aniso, {'eps': .15, 'n_neighbors': 2,
             'min_samples': 20, 'xi': 0.1, 'min_cluster_size': .2}),
    (blobs, {}),
    (no_structure, {})]

for i_dataset, (dataset, algo_params) in enumerate(datasets):
    # update parameters with dataset-specific values
    params = default_base.copy()
    params.update(algo_params)
    X, y = dataset

    # normalize dataset for easier parameter selection
    X = StandardScaler().fit_transform(X)

    # estimate bandwidth for mean shift
    bandwidth = cluster.estimate_bandwidth(X, quantile=params['quantile'])

    # connectivity matrix for structured Ward
    connectivity = kneighbors_graph(
        X, n_neighbors=params['n_neighbors'], include_self=False)
    # make connectivity symmetric
    connectivity = 0.5 * (connectivity + connectivity.T)
    # ============
    # Create cluster objects
    # ============
    ms = cluster.MeanShift(bandwidth=bandwidth, bin_seeding=True)
    two_means = cluster.MiniBatchKMeans(n_clusters=params['n_clusters'])
    ward = cluster.AgglomerativeClustering(
        n_clusters=params['n_clusters'], linkage='ward',
        connectivity=connectivity)
    spectral = cluster.SpectralClustering(
        n_clusters=params['n_clusters'], eigen_solver='arpack',
        affinity="nearest_neighbors")
    dbscan = cluster.DBSCAN(eps=params['eps'])
    optics = cluster.OPTICS(min_samples=params['min_samples'],
                            xi=params['xi'],
                            min_cluster_size=params['min_cluster_size'])
    affinity_propagation = cluster.AffinityPropagation(
        damping=params['damping'], preference=params['preference'],random_state=0)
    average_linkage = cluster.AgglomerativeClustering(
        linkage="average", affinity="cityblock",
        n_clusters=params['n_clusters'], connectivity=connectivity)
    birch = cluster.Birch(n_clusters=params['n_clusters'])
    gmm = mixture.GaussianMixture(
        n_components=params['n_clusters'], covariance_type='full')

    clustering_algorithms = (
        ('MiniBatchKMeans', two_means),
        ('AffinityPropagation', affinity_propagation),
        ('MeanShift', ms),
        ('SpectralClustering', spectral),
        ('Ward', ward),
        ('AgglomerativeClustering', average_linkage),
        ('DBSCAN', dbscan),
        ('OPTICS', optics),
        ('Birch', birch),
        ('GaussianMixture', gmm)
    )

    for name, algorithm in clustering_algorithms:
        t0 = time.time()

        # catch warnings related to kneighbors_graph
        with warnings.catch_warnings():
            warnings.filterwarnings(
                "ignore",
                message="the number of connected components of the " +
                "connectivity matrix is [0-9]{1,2}" +
                " > 1. Completing it to avoid stopping the tree early.",
                category=UserWarning)
            warnings.filterwarnings(
                "ignore",
                message="Graph is not fully connected, spectral embedding" +
                " may not work as expected.",
                category=UserWarning)
            algorithm.fit(X)

        t1 = time.time()
        if hasattr(algorithm, 'predict'):
            y_pre = algorithm.predict(X)
        else:
            y_pre = algorithm.fit_predict(X)

        plt.subplot(len(datasets), len(clustering_algorithms), plot_num)
        if i_dataset == 0:
            plt.title(name, size=18)
        colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a',
                                             '#f781bf', '#a65628', '#984ea3',
                                             '#999999', '#e41a1c', '#dede00']),
                                      int(max(y_pred) + 1))))
        # add black color for outliers (if any)
        colors = np.append(colors, ["#000000"])
        plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])

        plt.xlim(-2.5, 2.5)
        plt.ylim(-2.5, 2.5)
        plt.xticks(())
        plt.yticks(())
        plt.text(.99, .01, ('%.2fs' % (t1 - t0)).lstrip('0'),
                 transform=plt.gca().transAxes, size=15,
                 horizontalalignment='right')
        
        try:
            plt.text(.99, .70, ('sil:%.2f' % silhouette_score(X, y_pre)).lstrip('0'),
                 transform=plt.gca().transAxes, size=15,
                 horizontalalignment='right')
        except:
            pass
        plot_num += 1
plt.show()

在这里插入图片描述

3.2 例 :Fuzzy c-means clustering

from skfuzzy.cluster import cmeans
X,y=noisy_circles
alldata=np.c_[noisy_circles].T
print(alldata.shape)
center,u_orig, u0, d, jm, p, fpc = cmeans(alldata, m=2, c=4, error=0.005, maxiter=1000)

# Show n-cluster model
fig2, ax2 = plt.subplots()
ax2.set_title('Trained model')
for j in range(4):
    ax2.plot(alldata[0, u_orig.argmax(axis=0) == j],
             alldata[1, u_orig.argmax(axis=0) == j], 'o',
             label='series ' + str(j))
ax2.legend()

在这里插入图片描述

3.3 HMM隐马尔科夫模型

#coding=utf-8
'''
Created on 2017-12-4

本例为天气和行为的关系
'''
import numpy as np
import matplotlib.pyplot as plt
# hmmlearn可以在安装numpy以后,再使用pip install hmmlearn安装
from hmmlearn import hmm

states = ["Rainy", "Sunny"]##隐藏状态
n_states = len(states)##隐藏状态长度

observations = ["walk", "shop", "clean"]##可观察的状态
n_observations = len(observations)##可观察序列的长度

start_probability = np.array([0.6, 0.4])##开始转移概率,即开始是Rainy和Sunny的概率
##隐藏间天气转移混淆矩阵,即Rainy和Sunny之间的转换关系,例如[0,0]表示今天Rainy,明天Rainy的概率
transition_probability = np.array([
  [0.7, 0.3],
  [0.4, 0.6]
])
##隐藏状态天气和可视行为混淆矩阵,例如[0,0]表示今天Rainy,walk行为的概率为0.1
emission_probability = np.array([
  [0.1, 0.4, 0.5],
  [0.6, 0.3, 0.1]
])

#构建了一个MultinomialHMM模型,这模型包括开始的转移概率,隐藏间天气转换混淆矩阵(transmat),隐藏状态天气和可视行为混淆矩阵emissionprob,对模型参数初始化
model = hmm.MultinomialHMM(n_components=n_states)
model.startprob_= start_probability
model.transmat_ = transition_probability
model.emissionprob_ = emission_probability

#给出一个可见序列
bob_Actions = np.array([[2, 0, 1, 1, 2, 0]]).T

# 解决问题1,解码问题,已知模型参数和X,估计最可能的Z; 维特比算法 
logprob, weathers = model.decode(bob_Actions, algorithm="viterbi")
#print("Bob Actions:", ", ".join(map(lambda x: observations[x], bob_Actions)))
#print("weathers:", ", ".join(map(lambda x: states[x], weathers)))
print(logprob)#该参数反映模型拟合的好坏,数值越大越好
# 解决问题2,概率问题,已知模型参数和X,估计X出现的概率, 向前-向后算法 
score = model.score(bob_Actions, lengths=None)
#最后输出结果
print(score)
# 结果为-6.892170869,其实真正的概率是以自然数e为底数ln(P) = -6.892170869,所以概率P = 0.00101570648021
# import math
# print math.exp(-6.892170869)

# 解决问题3,学习问题,仅给出X,估计模型参数,鲍姆-韦尔奇算法,其实就是基于EM算法的求解
# 解决这个问题需要X的有一定的数据量,然后通过model.fit(X, lengths=None)来进行训练然后自己生成一个模型
# 并不需要设置model.startprob_,model.transmat_,model.emissionprob_
# 例如:
# import numpy as np
# from hmmlearn import hmm
#  
# states = ["Rainy", "Sunny"]##隐藏状态
# n_states = len(states)##隐藏状态长度
#  
# observations = ["walk", "shop", "clean"]##可观察的状态
# n_observations = len(observations)##可观察序列的长度
#  
# model = hmm.MultinomialHMM(n_components=n_states, n_iter=1000, tol=0.01)
#  
# X = np.array([[2, 0, 1, 1, 2, 0],[0, 0, 1, 1, 2, 0],[2, 1, 2, 1, 2, 0]])
# model.fit(X)
# print model.startprob_
# print model.transmat_
# print model.emissionprob_
## [[  1.11111111e-01   2.22222222e-01   6.66666667e-01]
##  [  5.55555556e-01   4.44444444e-01   6.27814351e-28]]
# print model.score(X)
# model.fit(X)
# print model.startprob_
# print model.transmat_
# print model.emissionprob_
# 和第一次fit(X)得到的行顺序不一样
## [[  5.55555556e-01   4.44444444e-01   9.29759770e-28]
##  [  1.11111111e-01   2.22222222e-01   6.66666667e-01]]
# print model.score(X)
# model.fit(X)
# print model.startprob_
# print model.transmat_
# print model.emissionprob_
# print model.score(X)
# # 可以进行多次fit,然后拿评分最高的模型,就可以预测了
# print model.predict(bob_Actions, lengths=None)
# # 预测最可能的隐藏状态
# # 例如:
# # [0 1 0 0 0 1]
# print model.predict_proba(bob_Actions, lengths=None)# 预测各个隐藏状态的概率
# # 例如:
# # [[ 0.82770645  0.17229355]
# #  [ 0.27361913  0.72638087]
# #  [ 0.58700959  0.41299041]
# #  [ 0.69861348  0.30138652]
# #  [ 0.81799813  0.18200187]
# #  [ 0.24723966  0.75276034]]
# # 在生成的模型中,可以随机生成随机生成一个模型的Z和X
# X,Z = model.sample(n_samples=5, random_state=None)
# print "Bob Actions:", ", ".join(map(lambda x: observations[x], X))
# print "weathers:", ", ".join(map(lambda x: states[x], Z))


# # 保存模型
# import pickle
# output = open('D:\\xxx\\data1111.pkl', 'wb')
# s = pickle.dump(model, output)
# output.close()
# # 调用模型
# input = open('D:\\xxx\\data.pkl', 'rb')
# model = pickle.load(model)
# input.close()
# model.predict(X)

3.4 原型聚类之学习向量量化(Learning Vector Quantization)—-有标签

import numpy as np
import copy
from sklearn.datasets import make_moons

from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
class LVQ():
    def __init__(self, max_iter=10000, eta=0.1, e=0.01):
        self.max_iter = max_iter
        self.eta = eta
        self.e = e

    def dist(self, x1, x2):
        return np.linalg.norm(x1 - x2)

    def get_mu(self, X, Y):
        k = len(set(Y))
        index = np.random.choice(X.shape[0], 1, replace=False)
        mus = []
        mus.append(X[index])
        mus_label = []
        mus_label.append(Y[index])
        for _ in range(k - 1):
            max_dist_index = 0
            max_distance = 0
            for j in range(X.shape[0]):
                min_dist_with_mu = 999999

                for mu in mus:
                    dist_with_mu = self.dist(mu, X[j])
                    if min_dist_with_mu > dist_with_mu:
                        min_dist_with_mu = dist_with_mu

                if max_distance < min_dist_with_mu:
                    max_distance = min_dist_with_mu
                    max_dist_index = j
            mus.append(X[max_dist_index])
            mus_label.append(Y[max_dist_index])

        mus_array = np.array([])
        for i in range(k):
            if i == 0:
                mus_array = mus[i]
            else:
                mus[i] = mus[i].reshape(mus[0].shape)
                mus_array = np.append(mus_array, mus[i], axis=0)
        mus_label_array = np.array(mus_label)
        return mus_array, mus_label_array

    def get_mu_index(self, x):
        min_dist_with_mu = 999999
        index = -1

        for i in range(self.mus_array.shape[0]):
            dist_with_mu = self.dist(self.mus_array[i], x)
            if min_dist_with_mu > dist_with_mu:
                min_dist_with_mu = dist_with_mu
                index = i

        return index

    def fit(self, X, Y):
        self.mus_array, self.mus_label_array = self.get_mu(X, Y)
        iter = 0

        while(iter < self.max_iter):
            old_mus_array = copy.deepcopy(self.mus_array)
            index = np.random.choice(Y.shape[0], 1, replace=False)

            mu_index = self.get_mu_index(X[index])
            if self.mus_label_array[mu_index] == Y[index]:
                self.mus_array[mu_index] = self.mus_array[mu_index] + \
                    self.eta * (X[index] - self.mus_array[mu_index])
            else:
                self.mus_array[mu_index] = self.mus_array[mu_index] - \
                    self.eta * (X[index] - self.mus_array[mu_index])

            diff = 0
            for i in range(self.mus_array.shape[0]):
                diff += np.linalg.norm(self.mus_array[i] - old_mus_array[i])
            if diff < self.e:
                print('迭代{}次退出'.format(iter))
                return
            iter += 1
        print("迭代超过{}次,退出迭代".format(self.max_iter))


if __name__ == '__main__':

    fig = plt.figure(1)

    plt.subplot(221)
    center = [[1, 1], [-1, -1], [1, -1]]
    cluster_std = 0.35
    X1, Y1 = make_blobs(n_samples=1000, centers=center,
                        n_features=2, cluster_std=cluster_std, random_state=1)
    plt.scatter(X1[:, 0], X1[:, 1], marker='o', c=Y1)

    plt.subplot(222)
    lvq1 = LVQ()
    lvq1.fit(X1, Y1)
    mus = lvq1.mus_array
    plt.scatter(X1[:, 0], X1[:, 1], marker='o', c=Y1)
    plt.scatter(mus[:, 0], mus[:, 1], marker='^', c='r')

    plt.subplot(223)
    X2, Y2 = make_moons(n_samples=1000, noise=0.1)
    plt.scatter(X2[:, 0], X2[:, 1], marker='o', c=Y2)

    plt.subplot(224)
    lvq2 = LVQ()
    lvq2.fit(X2, Y2)
    mus = lvq2.mus_array
    plt.scatter(X2[:, 0], X2[:, 1], marker='o', c=Y2)
    plt.scatter(mus[:, 0], mus[:, 1], marker='^', c='r')
    plt.show()

在这里插入图片描述

3.5 竞争网络

import numpy as np
import matplotlib.pyplot as plt
 
#定义激活函数
def sigmoid(x):
    return 1/(1+np.exp(-x))
 
#对一维向量归一化
def normalization(M):
    """
    M行向量
    """
    M = M/(np.sqrt(np.dot(M,M.T)))
    return M
 
#对矩阵进行归一化
def normalization_all(N):
    """
    对矩阵进行归一化
    N代表mxn的矩阵
    返回一个归一化以后的矩阵
    """
    N_all=[]
    for i in range(len(N)):   #len(N):N的行数
        K = normalization(N[i])
        N_all.append(K)
    return N_all
 
class competitive_network(object):
    """
    竞争神经网络
    """
    def __init__(self, num_in, num_out, lr):
        """
        初始化参数分别为输入层节点数、输出层或者是竞争层节点数、学习率
        """
        w = np.random.rand(num_out, num_in)*(-2) + 1   #随机生成一个权值矩阵,取值范围是[-1,1];
        self.w = normalization_all(w)
        self.lr = lr
        
    #信号向前传播
    def forward_propagation(self, x):
        x= x.reshape(1, x.shape[0])   #shape函数,输出一维的值也就是X的行数;
                                      #reshape()是数组对象中的方法,用于改变数组的形状。
                                      #reshape函数创建一个行数为1,列数为X矩阵的行数。
        #寻找获胜神经元
        y = np.dot(self.w, x.T)   #输入矩阵和权值向量的点积
        o = sigmoid(y)   #神经元输出值
        argmax = np.where(o == np.amax(o))[0][0]
        return argmax     #返回获胜神经元的位置
    
    #反向传播,求权值矩阵;
    def back_propagation(self, argmax, x):   #x是经过归一化处理的
        self.w[argmax] = self.lr * (x - self.w[argmax])
        self.w[argmax] = normalization(self.w[argmax])
        self.lr -= self.decay   #学习率逐渐减小
        
    #训练权重
    def train(self, x, num_item):   #num_item:迭代次数,x:输入矩阵
        x = np.array(x)
        self.decay = self.lr/num_item
        for i in range(num_item):   #迭代一定次数
            for j in range(x.shape[0]):   #对于输入模式的行数进行循环
                ar = self.forward_propagation(x[j])
                self.back_propagation(ar, x[j]) 
        for i in range(2):
            for j in range(2):
                print(self.w[i][j])   #输出权值                                       
    
    def prediction(self, x):
        argmax = self.forward_propagation(x)
        return argmax
    
datamat = np.random.rand(100,2) * (-2) + 1    #随机生成的输入矩阵
print(datamat)
 
#获胜神经元
c = []
 
CN = competitive_network(2, 2, 0.1)   #建立神经网络
CN.train(datamat, 1000)    #训练权重向量
 
for i in range(len(datamat)):
    prediction=CN.prediction(datamat[i])
    c.append(prediction*1)
print(c)   #输出100个输入数据的输出值0或者1;
 
datamat=normalization_all(datamat)
datamat=np.array(datamat)   #参考:https://blog.csdn.net/qq_38237214/article/details/76851107
plt.figure()
plt.scatter(datamat[:,0],datamat[:,1], c=c)
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.show()

3.6 自组织特征映射神经网络SOM

例一:
import numpy as np
import pylab as pl
 
#神经网络 
class SOM(object):
  def __init__(self, X, output, iteration, batch_size):
    """
    :param X: 形状是N*D,输入样本有N个,每个D维
    :param output: (n,m)一个元组,为输出层的形状是一个n*m的二维矩阵
    :param iteration:迭代次数
    :param batch_size:每次迭代时的样本数量
    初始化一个权值矩阵,形状为D*(n*m),即有n*m权值向量,每个D维
    """
    self.X = X
    self.output = output
    self.iteration = iteration
    self.batch_size = batch_size
    self.W = np.random.rand(X.shape[1], output[0] * output[1])
    print (self.W.shape)
    print("权值矩阵:",self.W)
 
  def GetN(self, t):
    """
    :param t:时间t, 这里用迭代次数来表示时间
    :return: 返回一个整数,表示拓扑距离,时间越大,拓扑邻域越小
    """
    a = min(self.output)
    return int(a-float(a)*t/self.iteration)
 
  #求学习率
  def Geteta(self, t, n):
    """
    :param t: 时间t, 这里用迭代次数来表示时间
    :param n: 拓扑距离
    :return: 返回学习率,
    """
    return np.power(np.e, -n)/(t+2)
 
#更新权值矩阵 
  def updata_W(self, X, t, winner):
    N = self.GetN(t)   #表示随时间变化的拓扑距离
    for x, i in enumerate(winner):
      to_update = self.getneighbor(i[0], N)
      for j in range(N+1):
        e = self.Geteta(t, j)   #表示学习率
        for w in to_update[j]:
          self.W[:, w] = np.add(self.W[:,w], e*(X[x,:] - self.W[:,w]))
 
  def getneighbor(self, index, N):
    """
    :param index:获胜神经元的下标
    :param N: 邻域半径
    :return ans: 返回一个集合列表,分别是不同邻域半径内需要更新的神经元坐标
    """
    a, b = self.output
    length = a*b
    def distence(index1, index2):
      i1_a, i1_b = index1 // a, index1 % b   #//:向下取整; %:返回除法的余数;
      i2_a, i2_b = index2 // a, index2 % b
      return np.abs(i1_a - i2_a), np.abs(i1_b - i2_b)   #abs() 函数返回数字的绝对值。
 
    ans = [set() for i in range(N+1)]
    for i in range(length):
      dist_a, dist_b = distence(i, index)
      if dist_a <= N and dist_b <= N: ans[max(dist_a, dist_b)].add(i)
    return ans

  def train(self):
    """
    train_Y:训练样本与形状为batch_size*(n*m)
    winner:一个一维向量,batch_size个获胜神经元的下标
    :return:返回值是调整后的W
    """
    count = 0
    while self.iteration > count:
      train_X = self.X[np.random.choice(self.X.shape[0], self.batch_size)]
      normal_W(self.W)
      normal_X(train_X)
      train_Y = train_X.dot(self.W)
      winner = np.argmax(train_Y, axis=1).tolist()
      self.updata_W(train_X, count, winner)
      count += 1
    return self.W
 
  def train_result(self):
    normal_X(self.X)
    train_Y = self.X.dot(self.W)
    winner = np.argmax(train_Y, axis=1).tolist()
    print (winner)
    return winner
 
def normal_X(X):
  """
  :param X:二维矩阵,N*D,N个D维的数据
  :return: 将X归一化的结果
  """
  N, D = X.shape
  for i in range(N):
    temp = np.sum(np.multiply(X[i], X[i]))
    X[i] /= np.sqrt(temp)
  return X
def normal_W(W):
  """
  :param W:二维矩阵,D*(n*m),D个n*m维的数据
  :return: 将W归一化的结果
  """
  for i in range(W.shape[1]):
    temp = np.sum(np.multiply(W[:,i], W[:,i]))
    W[:, i] /= np.sqrt(temp)
  return W
 
#画图
def draw(C):
  colValue = ['r', 'y', 'g', 'b', 'c', 'k', 'm']
  for i in range(len(C)):
    coo_X = []  #x坐标列表
    coo_Y = []  #y坐标列表
    for j in range(len(C[i])):
      coo_X.append(C[i][j][0])
      coo_Y.append(C[i][j][1])
    pl.scatter(coo_X, coo_Y, marker='x', color=colValue[i%len(colValue)], label=i)
 
  pl.legend(loc='upper right')   #图例位置
  pl.show()
 
#数据集:每三个是一组分别是西瓜的编号,密度,含糖量
data = """
1,0.697,0.46,2,0.774,0.376,3,0.634,0.264,4,0.608,0.318,5,0.556,0.215,
6,0.403,0.237,7,0.481,0.149,8,0.437,0.211,9,0.666,0.091,10,0.243,0.267,
11,0.245,0.057,12,0.343,0.099,13,0.639,0.161,14,0.657,0.198,15,0.36,0.37,
16,0.593,0.042,17,0.719,0.103,18,0.359,0.188,19,0.339,0.241,20,0.282,0.257,
21,0.748,0.232,22,0.714,0.346,23,0.483,0.312,24,0.478,0.437,25,0.525,0.369,
26,0.751,0.489,27,0.532,0.472,28,0.473,0.376,29,0.725,0.445,30,0.446,0.459"""
 
a = data.split(',')   #split() 通过指定分隔符对字符串进行切片,就是以逗号为分割依据,进行分割;
dataset = np.mat([[float(a[i]), float(a[i+1])] for i in range(1, len(a)-1, 3)])  #用mat函数转换为矩阵之后可以才进行一些线性代数的操作。
dataset_old = dataset.copy()
print("输入数据:", dataset)
 
som = SOM(dataset, (5, 5), 1, 30)
som.train()
res = som.train_result()
classify = {}
for i, win in enumerate(res):
  if not classify.get(win[0]):
    classify.setdefault(win[0], [i])
  else:
    classify[win[0]].append(i)
C = []#未归一化的数据分类结果
D = []#归一化的数据分类结果
for i in classify.values():
  C.append(dataset_old[i].tolist())
  D.append(dataset[i].tolist())
 
draw(C)
draw(D)
例二:
import numpy as np
import random
np.random.seed(22)
class CyrusSOM(object):
    def __init__(self,net=[[1,1],[1,1]],epochs = 50,r_t = [None,None],eps=1e-6):
        """
        :param net: 竞争层的拓扑结构,支持一维及二维,1表示该输出节点存在,0表示不存在该输出节点
        :param epochs: 最大迭代次数
        :param r_t:   [C,B]    领域半径参数,r = C*e**(-B*t/eoochs),其中t表示当前迭代次数
        :param eps: learning rate的阈值
        """
        self.epochs = epochs
        self.C = r_t[0]
        self.B = r_t[1]
        self.eps = eps
        self.output_net = np.array(net)
        if len(self.output_net.shape) == 1:
            self.output_net = self.output_net.reshape([-1,1])
        self.coord = np.zeros([self.output_net.shape[0],self.output_net.shape[1],2])
        for i in range(self.output_net.shape[0]):
            for j in range(self.output_net.shape[1]):
                self.coord[i,j] = [i,j]
        print(self.coord)


    def __r_t(self,t):
        if not self.C:
            return 0.5
        else:
            return self.C*np.exp(-self.B*t/self.epochs)

    def __lr(self,t,distance):
        return (self.epochs-t)/self.epochs*np.exp(-distance)
    def standard_x(self,x):
        x = np.array(x)
        for i in range(x.shape[0]):
            x[i,:] = [value/(((x[i,:])**2).sum()**0.5) for value in x[i,:]]
        return x
    def standard_w(self,w):
        for i in range(w.shape[0]):
            for j in range(w.shape[1]):
                w[i,j,:] = [value/(((w[i,j,:])**2).sum()**0.5) for value in w[i,j,:]]
        return w
    def cal_similar(self,x,w):
        similar = (x*w).sum(axis=2)
        coord = np.where(similar==similar.max())
        return [coord[0][0],coord[1][0]]

    def update_w(self,center_coord,x,step):
        for i in range(self.coord.shape[0]):
            for j in range(self.coord.shape[1]):
                distance = (((center_coord-self.coord[i,j])**2).sum())**0.5
                if distance <= self.__r_t(step):
                    self.W[i,j] = self.W[i,j] + self.__lr(step,distance)*(x-self.W[i,j])

    def transform_fit(self,x):
        self.train_x = self.standard_x(x)
        self.W = np.zeros([self.output_net.shape[0],self.output_net.shape[1],self.train_x.shape[1]])
        for i in range(self.W.shape[0]):
            for j in range(self.W.shape[1]):
                self.W[i,j,:] = self.train_x[random.choice(range(self.train_x.shape[0])),:]
        self.W = self.standard_w(self.W)
        for step in range(int(self.epochs)):
            j = 0
            if self.__lr(step,0) <= self.eps:
                break
            for index in range(self.train_x.shape[0]):
                print("*"*8,"({},{})/{} W:\n".format(step,j,self.epochs),self.W)
                center_coord = self.cal_similar(self.train_x[index,:],self.W)
                self.update_w(center_coord,self.train_x[index,:],step)
                self.W = self.standard_w(self.W)
                j += 1
        label = []
        for index in range(self.train_x.shape[0]):
            center_coord = self.cal_similar(self.train_x[index, :], self.W)
            label.append(center_coord[1]*self.coord.shape[1] + center_coord[0])
        class_dict = {}
        for index in range(self.train_x.shape[0]):
            if label[index] in class_dict.keys():
                class_dict[label[index]].append(index)
            else:
                class_dict[label[index]] = [index]
        cluster_center = {}
        for key,value in class_dict.items():
            cluster_center[key] = np.array([x[i, :] for i in value]).mean(axis=0)
        self.cluster_center = cluster_center
        return label

    def fit(self,x):
        self.train_x = self.standard_x(x)
        self.W = np.random.rand(self.output_net.shape[0], self.output_net.shape[1], self.train_x.shape[1])
        self.W = self.standard_w(self.W)
        for step in range(int(self.epochs)):
            j = 0
            if self.__lr(step,0) <= self.eps:
                break
            for index in range(self.train_x.shape[0]):
                #print("*"*8,"({},{})/{} W:\n".format(step, j, self.epochs), self.W)
                center_coord = self.cal_similar(self.train_x[index, :], self.W)
                self.update_w(center_coord, self.train_x[index, :], step)
                self.W = self.standard_w(self.W)
                j += 1
        label = []
        for index in range(self.train_x.shape[0]):
            center_coord = self.cal_similar(self.train_x[index, :], self.W)
            label.append(center_coord[1] * self.coord.shape[1] + center_coord[1])
        class_dict = {}
        for index in range(self.train_x.shape[0]):
            if label[index] in class_dict.keys():
                class_dict[label[index]].append(index)
            else:
                class_dict[label[index]] = [index]
        cluster_center = {}
        for key, value in class_dict.items():
            cluster_center[key] = np.array([x[i, :] for i in value]).mean(axis=0)
        self.cluster_center = cluster_center

    def predict(self,x):
        self.pre_x = self.standard_x(x)
        label = []
        for index in range(self.pre_x.shape[0]):
            center_coord = self.cal_similar(self.pre_x[index, :], self.W)
            label.append(center_coord[1] * self.coord.shape[1] + center_coord[1])
        return label
例三:
from numpy import *
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
class SOMnet(object):
    def __init__(self):      # 设计网络参数初始值
        self.lratemax = 0.8; # 最大学习率--欧式距离
        self.lratemin = 0.05 # 最小学习率--欧式距离
        self.rmax = 5.0      # 最大聚类半径--根据数据集
        self.rmin = 0.5      # 最小聚类半径--根据数据集
        self.Steps = 1000    # 迭代次数
        self.lratrlist = []  # 学习率收敛曲线
        self.rlist = []      # 学习半径收敛曲线
        self.w = []          # 权重向量组
        self.M = 2           # M*N 表示聚类总数
        self.N = 2           # M,N表示领域的参数
        self.dataMat = []    # 外部导入数据集
        self.classLabel = [] # 聚类后的类别标签
    """
    def loadDataSet(self,filename): # 加载数据文件
        numFeat = len(open(filename).readline().split('\t ')) -1
        fr = open(filename)
        for line in fr.readline():
            lineArr = []
            curLine = line.strip().split("\t")
            lineArr.append(float(curLine[0]))
            lineArr.append(float(curLine[1]))
            #print(shape(lineArr))
            self.dataMat.append(lineArr)
        self.dataMat = mat(self.dataMat)
    """
    def loadDataSet(self, iris):  # 加载数据集
        # 将sklearn中的数据集datasets.load_iris()的特征作分类,每班样品50,样品总数150,维数4
        self.dataMat = iris["data"]
        self.dataMat = mat(self.dataMat)

    """
    def file2matrix(self,path,delimiter):
        recordlsit = []
        fp = open(path)
        content = fp.read()
        fp.close()
        rowlist = content.splitlines()#按行转换为一维表
        # 逐行遍历 结果按分隔符分隔为行向量
        recordlsit = list(map(eval,row.split(delimiter)) for row in rowlist if row.strip())
        self.dataMat = mat(recordlsit)
    """
    def normalize(self, dataMat):  # 数据归一化
        [m, n] = shape(dataMat)
        for i in arange(n - 1):
            dataMat[:, i] = (dataMat[:, i] - mean(dataMat[:, i])) / (std(dataMat[:, i]) + 1.0e-10)
        return dataMat

    def distEclud(self, matA, matB):  # 计算欧式距离
        ma,na = shape(matA)
        mb,nb = shape(matB)
        rtnmat = zeros((ma,nb))
        for i in arange(ma):
            for j in arange(nb):
                rtnmat[i,j]=np.linalg.norm(matA[i,:]-matB[:,j].T)
        return rtnmat


    def init_grid(self):  # 初始化第二层网络
        k=0 # 构建第二层网络模型
        grid = mat(zeros((self.M * self.N,2)))
        for i in arange(self.M):
            for j in arange(self.N):
                grid[k,:] = [i,j]
                k += 1
        return grid

    def ratecalc(self,i): # 学习率 和 聚类半径
        lrate = self.lratemax - ((i+1.0) * (self.lratemax - self.lratemin)) / self.Steps
        r = self.rmax - ((i+1.0) * (self.rmax - self.rmin)) / self.Steps
        return lrate,r

    def trainSOM(self): # SOM网络的实现
        dm,dn = shape(self.dataMat) # 1. 构建输入层网络
        normDataset = self.normalize(self.dataMat) # 归一化数据
        grid = self.init_grid() # 2. 初始化第二层分类网络
        self.w = np.random.rand(dn,self.M*self.N) # 3. 随机初始化两层之间的权值向量
        distM = self.distEclud # 确定距离公式
        # 4. 迭代求解
        if self.Steps < 5*dm:self.Steps = 5*dm  # 设定最小迭代次数
        for i in arange(self.Steps):
            lrate,r = self.ratecalc(i) # 1) 计算当前迭代次数下的学习率和分类半径
            self.lratrlist.append(lrate);self.rlist.append(r)
            # 2) 随机生成样本索引,并抽取一个样本
            k = np.random.randint(0,dm)
            mySample = normDataset[k,:]
            # 3) 计算最优节点:返回最小距离的索引值
            minIndx = (distM(mySample,self.w)).argmin()
            # 4) 计算领域
            d1 = ceil(minIndx/self.M) # 计算此节点在第二层矩阵中的位置
            d2 = mod(minIndx,self.M)
            distMat = distM(mat([d1,d2]),grid.T)
            #nodelindx = (distMat < r).nonzeor()[1] # 获取领域内的所有节点
            nodelindx = np.nonzero((distMat < r))[1] # 获取领域内的所有节点
            for j in arange(shape(self.w)[1]): # 5) 案列更新权重
                if sum(nodelindx == j):
                    self.w[:,j] = self.w[:,j] + lrate * (mySample[0]-self.w[:,j])
        self.classLabel = list(range(dm)) # 分配和存储聚类后的类别标签
        for i in arange(dm):
            dist=distM(normDataset[i, :], self.w)
            print(type(dist))
            #self.classLabel[i] = argmin(distM(normDataset[i,:],self.w))# np.argmin()求最小值的坐标
            #self.classLabel[i] = distM(normDataset[i,:],self.w).argmin()# np.argmin()求最小值的坐标
            self.classLabel[i] = np.argmin(dist)
        self.classLabel = mat(self.classLabel)

    def showCluster(self,plt): # 绘图  显示聚类结果
        lst = unique(self.classLabel.tolist()[0]) # 去除
        i = 0
        for cindx in lst:
            myclass = nonzero(self.classLabel == cindx)[1]
            xx = self.dataMat[myclass].copy()
            if i == 0:plt.plot(xx[:,0],xx[:,1],'bo')
            if i == 1:plt.plot(xx[:,0],xx[:,1],'rd')
            if i == 2: plt.plot(xx[:, 0], xx[:, 1], 'gD')
            if i == 3: plt.plot(xx[:, 0], xx[:, 1], 'c^')
            i += 1
        plt.show()
if __name__ == "__main__":
    # 加载
    SOMNet =  SOMnet()
    iris = load_iris()
    SOMNet.loadDataSet(iris)
    #SOMNet.file2matrix("test.txt",'/t')
    SOMNet.trainSOM()
    print(SOMNet.w)
SOMNet.showCluster(plt)

4 集成学习

4.1 基础的集成技术—最大投票法

from sklearn.ensemble import VotingClassifier
model1 = LogisticRegression(random_state=1)
model2 = tree.DecisionTreeClassifier(random_state=1)
model = VotingClassifier(estimators=[('lr', model1), ('dt', model2)], voting='hard')
model.fit(x_train,y_train)
model.score(x_test,y_test)

4.2 基础的集成技术—平均法

4.2.1 自定义平均法

例一:

model1 = tree.DecisionTreeClassifier()
model2 = KNeighborsClassifier()
model3= LogisticRegression()
model1.fit(x_train,y_train)
model2.fit(x_train,y_train)
model3.fit(x_train,y_train)
pred1=model1.predict_proba(x_test)
pred2=model2.predict_proba(x_test)
pred3=model3.predict_proba(x_test)
finalpred=(pred1+pred2+pred3)/3

例二:

#取平均
# get predictions for test 
result = pipeline.mean().execute()
# or Validate 
_ = pipeline.mean().validate(mean_absolute_error,10)
#自定义
result = pipeline.apply(lambda x: np.max(x,axis=0)).execute()

4.2.1 自定义加权平均法

例一:

model1 = tree.DecisionTreeClassifier()
model2 = KNeighborsClassifier()
model3= LogisticRegression()
model1.fit(x_train,y_train)
model2.fit(x_train,y_train)
model3.fit(x_train,y_train)
pred1=model1.predict_proba(x_test)
pred2=model2.predict_proba(x_test)
pred3=model3.predict_proba(x_test)
finalpred=(pred1*0.3+pred2*0.3+pred3*0.4)

例二:

from heamy.dataset import Dataset
from heamy.estimator import Regressor, Classifier
from heamy.pipeline import ModelsPipeline
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error
from sklearn.neighbors import KNeighborsRegressor
data = load_boston()
X, y = data['data'], data['target']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, random_state=111)
#创建数据集
dataset = Dataset(X_train,y_train,X_test)

model_rf = Regressor(dataset=dataset, estimator=RandomForestRegressor, parameters={'n_estimators': 151},name='rf')
model_lr = Regressor(dataset=dataset, estimator=LinearRegression, parameters={'normalize': True},name='lr')
model_knn = Regressor(dataset=dataset, estimator=KNeighborsRegressor, parameters={'n_neighbors': 15},name='knn')

pipeline = ModelsPipeline(model_rf,model_lr,model_knn)

weights = pipeline.find_weights(mean_absolute_error)
result = pipeline.weight(weights)
results=result.execute() ##预测的结果
print(results.shape,y_test.shape)
print('mean_absolute_error:',mean_absolute_error(y_test,results))
print(sum(results==y_test)/y_test.shape[0])
print(y_test[:10],results[:10])

4.3 Stacking

高级的集成技术—堆叠: 第一步:把训练集分成10份 第二步:基础模型(假设是决策树)在其中9份上拟合,并对第10份进行预测。 第三步:对训练集上的每一份如此做一遍。 第四步:然后将基础模型(此处是决策树)拟合到整个训练集上。 第五步:使用此模型,在测试集上进行预测。 第六步:对另一个基本模型(比如knn)重复步骤2到4,产生对训练集和测试集的另一组预测。 第七步:训练集预测被用作构建新模型的特征。 第八步:该新模型用于对测试预测集(test prediction set)进行最终预测。

4.3.1 自定义stacking

from sklearn.model_selection import StratifiedKFold
class stack_model():
    def __init__(self,classifiers,stack_clf):
        self.classifiers=classifiers
        self.stack_clf=stack_clf
    def fit(self,X,y,df_init=False,n_folds=5):
        train_df=[]
        y_df=[]
        if np.all(df_init):
            train_df=[pd.DataFrame(X)]
            y_df=[pd.DataFrame(y)]
        skf = StratifiedKFold(n_folds)   # 初始化
        skf.get_n_splits(X, y)
        for name, clf in self.classifiers.items():
            '''依次训练各个单模型'''
            pred_list=[]
            score_sum=0
            for i, (train, test) in enumerate(skf.split(X, y)):
                '''使用第i个部分作为预测,剩余的部分来训练模型,获得其预测的输出作为第i部分的新特征。'''
                X_train, y_train, X_test, y_test = X[train], y[train], X[test], y[test]
                clf.fit(X_train, y_train)
                score_sum+=clf.score(X_test, y_test)
                pred_list.extend(clf.predict(X_test))
            train_df.append(pd.DataFrame(pred_list))
            pred_list=[]
            score_list=0
            print(name,'score:',score_sum/n_folds)
        df_train=pd.concat(train_df,axis=1)
        self.stack_clf.fit(df_train,y)
        self.pred=self.stack_clf.predict(df_train)
        self.score_value=sum(self.pred==y)/len(y)
        print('train score:',self.stack_clf.score(df_train,y))
    def predict(self,x,df_init=False):
        train_df=[]
        if np.all(df_init):
            train_df=[pd.DataFrame(x)]
        for name, clf in self.classifiers.items():
            train_df.append(pd.DataFrame(clf.predict(x)))
        df_train=pd.concat(train_df,axis=1)
        x_pred=self.stack_clf.predict(df_train)
        return x_pred
    def score(self,X_train,y_train):
        self.fit(X_train,y_train)
        return self.score_value

4.3.2 StackingClassifier(可多层)

stacking严格来说并不是一种算法,而是精美而又复杂的,对模型集成的一种策略。

1、首先我们会得到两组数据:训练集和测试集。将训练集分成5份:train1,train2,train3,train4,train5。

2、选定基模型。这里假定我们选择了xgboost, lightgbm 和 randomforest 这三种作为基模型。比如xgboost模型部分:依次用train1,train2,train3,train4,train5作为验证集,其余4份作为训练集,进行5折交叉验证进行模型训练;再在测试集上进行预测。这样会得到在训练集上由xgboost模型训练出来的5份predictions,和在测试集上的1份预测值B1。将这五份纵向重叠合并起来得到A1。lightgbm和randomforest模型部分同理。

3、三个基模型训练完毕后,将三个模型在训练集上的预测值作为分别作为3个”特征”A1,A2,A3,使用LR模型进行训练,建立LR模型。

4、使用训练好的LR模型,在三个基模型之前在测试集上的预测值所构建的三个”特征”的值(B1,B2,B3)上,进行预测,得出最终的预测类别或概率。

做stacking,首先需要安装mlxtend库。安装方法:进入Anaconda Prompt,输入命令 pip install mlxtend 即可。

stacking主要有几种使用方法: 1、最基本的使用方法,即使用基分类器所产生的预测类别作为meta-classifier“特征”的输入数据——>use_probas = False

2、这一种是使用第一层所有基分类器所产生的类别概率值作为meta-classfier的输入。需要在StackingClassifier 中增加一个参数设置:use_probas = True。

另外,还有一个参数设置average_probas = True,那么这些基分类器所产出的概率值将按照列被平均,否则会拼接。

from sklearn import datasets
iris = datasets.load_iris()
X, y = iris.data[:, 1:3], iris.target
from sklearn import model_selection
from sklearn.linear_model import LogisticRegression
#pip install xgboost -i https://pypi.tuna.tsinghua.edu.cn/simple
from xgboost.sklearn import XGBClassifier
import lightgbm as lgb
from sklearn.ensemble import RandomForestClassifier
from mlxtend.classifier import StackingClassifier
import numpy as np
basemodel1 = XGBClassifier()
basemodel2 = lgb.LGBMClassifier()
basemodel3 = RandomForestClassifier(random_state=1)
lr = LogisticRegression()
sclf = StackingClassifier(classifiers=[basemodel1, basemodel2, basemodel3], 
                           use_probas=True,
                           average_probas=False,
                           meta_classifier=lr)
 
print('5-fold cross validation:\n')
for basemodel, label in zip([basemodel1, basemodel2, basemodel3, sclf], 
                      ['xgboost', 
                       'lightgbm', 
                       'Random Forest',
                       'StackingClassifier']):
 
    scores = model_selection.cross_val_score(basemodel,X, y, 
                                              cv=5, scoring='accuracy')
    print("Accuracy: %0.2f (+/- %0.2f) [%s]" 
          % (scores.mean(), scores.std(), label))

注:

这一种方法是对基分类器训练的特征维度进行操作的,并不是给每一个基分类器全部的特征,而是赋予不同的基分类器不同的特征。比如:基分类器1训练前半部分的特征,基分类器2训练后半部分的特征。这部分的操作是通过sklearn中的pipelines实现。最终通过StackingClassifier组合起来。

from sklearn.datasets import load_iris
from mlxtend.classifier import StackingClassifier
from mlxtend.feature_selection import ColumnSelector
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
from xgboost.sklearn import XGBClassifier
from sklearn.ensemble import RandomForestClassifier
iris = load_iris()
X = iris.data
y = iris.target
#基分类器1:xgboost
pipe1 = make_pipeline(ColumnSelector(cols=(0, 2)),
                      XGBClassifier())
#基分类器2:RandomForest
pipe2 = make_pipeline(ColumnSelector(cols=(1, 2, 3)),
                     RandomForestClassifier())
sclf = StackingClassifier(classifiers=[pipe1, pipe2], 
                          meta_classifier=LogisticRegression())
sclf.fit(X, y)

4.3.3 heamy(可多层)

from sklearn.linear_model import LogisticRegression
from heamy.dataset import Dataset
from heamy.estimator import Regressor, Classifier
from heamy.pipeline import ModelsPipeline
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error
#加载数据集
from sklearn.datasets import load_boston
from sklearn import datasets
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import VotingClassifier
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
import numpy as np
iris = datasets.load_iris()
X, Y = iris.data[:, 1:3], iris.target
# Split the data into a training set and a test set
X_train, x_test, Y_train, y_test = train_test_split(X, Y,test_size=0.2, random_state=0,stratify=Y)

#创建数据集
dataset = Dataset(X_train,Y_train,x_test)
# 1st level
model_rf = Classifier(dataset=dataset, estimator=RandomForestClassifier, parameters={'n_estimators': 151},name='rf')
#model_lr = Classifier(dataset=dataset, estimator=XGBClassifier, parameters={'use_label_encoder': False},name='xgb')
model_knn = Classifier(dataset=dataset, estimator=KNeighborsClassifier, parameters={'n_neighbors': 15},name='knn')
 
pipeline = ModelsPipeline(model_rf,model_knn)
stack_ds = pipeline.stack(k=5,seed=111)
 
# 2nd level
stack_rf = Classifier(dataset=stack_ds, estimator=RandomForestClassifier, parameters={'n_estimators': 15},name='rf')
stack_lr = Classifier(dataset=stack_ds, estimator=KNeighborsClassifier, parameters={'n_neighbors': 10},name='kb')
stack_pipeline = ModelsPipeline(stack_rf,stack_lr)
stack_ms=pipeline.stack(k=5,seed=0)
# 3rd level
stacker = Classifier(dataset=stack_ds, estimator=KNeighborsClassifier)
results = stacker.predict()
results=np.argmax(results,axis=1)

# 使用10折交叉验证结果
#results10 = stacker.validate(k=10,scorer=mean_absolute_error)

print('mean_absolute_error:',mean_absolute_error(y_test,results))
print(sum(results==y_test)/y_test.shape[0])

4.4 Blending

高级的集成技术–混合 混合遵循与堆叠相同的方法,但仅使用来自训练集的一个留出(holdout)/验证集来进行预测。换句话说,与堆叠不同,预测仅在留出集上进行。留出集和预测用于构建在测试集上运行的模型。以下是混合过程的详细说明: 第一步:原始训练数据被分为训练集合验证集。 第二步:在训练集上拟合模型。 第三步:在验证集和测试集上进行预测。 第四步:验证集及其预测用作构建新模型的特征。 第五步:该新模型用于对测试集和元特征(meta-features)进行最终预测。

例一:

#我们将在训练集上建立两个模型,决策树和knn,以便对验证集进行预测。
model1 = tree.DecisionTreeClassifier()
model1.fit(x_train, y_train)
val_pred1=model1.predict(x_val)
test_pred1=model1.predict(x_test)
val_pred1=pd.DataFrame(val_pred1)
test_pred1=pd.DataFrame(test_pred1)
model2 = KNeighborsClassifier()
model2.fit(x_train,y_train)
val_pred2=model2.predict(x_val)
test_pred2=model2.predict(x_test)
val_pred2=pd.DataFrame(val_pred2)
test_pred2=pd.DataFrame(test_pred2)
#结合元特征和验证集,构建逻辑回归模型以对测试集进行预测。
df_val=pd.concat([x_val, val_pred1,val_pred2],axis=1)
df_test=pd.concat([x_test, test_pred1,test_pred2],axis=1)
model = LogisticRegression()
model.fit(df_val,y_val)
model.score(df_test,y_test)

例二:

from heamy.dataset import Dataset
from heamy.estimator import Regressor, Classifier
from heamy.pipeline import ModelsPipeline
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error
#加载数据集
from sklearn.datasets import load_boston
data = load_boston()
X, y = data['data'], data['target']
X_train, X_test, y_train, y_test =train_test_split(X, y, test_size=0.1, random_state=111)
#创建数据集
dataset = Dataset(X_train,y_train,X_test)
#创建RF模型和LR模型
model_rf = Regressor(dataset=dataset, estimator=RandomForestRegressor, parameters={'n_estimators': 50},name='rf')
model_lr = Regressor(dataset=dataset, estimator=LinearRegression, parameters={'normalize': True},name='lr')
# Blending两个模型
# Returns new dataset with out-of-fold predictions
pipeline = ModelsPipeline(model_rf,model_lr)
stack_ds = pipeline.blend(proportion=0.2,seed=111)
#第二层使用lr模型stack
stacker = Regressor(dataset=stack_ds, estimator=LinearRegression)
results = stacker.predict()
# 使用10折交叉验证结果
results10 = stacker.validate(k=10,scorer=mean_absolute_error)

4.5 同异质集成模型

4.5.1 Bagging

Bagging背后的想法是结合多个模型的结果(例如,所有决策树)来获得泛化的结果。这有一个问题:如果在同样一组数据上创建所有模型并将其组合起来,它会有用吗?这些模型极大可能会得到相同的结果,因为它们获得的输入相同。那我们该如何解决这个问题呢?其中一种技术是自举(bootstrapping)。 Bootstrapping是一种采样技术,我们有放回的从原始数据集上创建观察子集,子集的大小与原始集的大小相同。 Bagging(或Bootstrap Aggregating)技术使用这些子集(包)来获得分布的完整概念(完备集)。为bagging创建的子集的大小也可能小于原始集。 第一步:从原始数据集有放回的选择观测值来创建多个子集。 第二步:在每一个子集上创建一个基础模型(弱模型)。 第三步:这些模型同时运行,彼此独立。 第四步:通过组合所有模型的预测来确定最终预测。

Bagging 算法: Bagging meta-estimator, 随机森林

在我们进一步讨论之前,这里有另一个问题:如果第一个模型错误地预测了某一个数据点,然后接下来的模型(可能是所有模型),将预测组合起来会提供更好的结果吗?Boosting就是来处理这种情况的。 Boosting是一个顺序过程,每个后续模型都会尝试纠正先前模型的错误。后续的模型依赖于之前的模型。接下来一起看看boosting的工作方式: 第一步:从原始数据集创建一个子集。 第二步:最初,所有数据点都具有相同的权重。 第三步:在此子集上创建基础模型。 第四步:该模型用于对整个数据集进行预测。 第五步:使用实际值和预测值计算误差。 第六步:预测错误的点获得更高的权重。(这里,三个错误分类的蓝色加号点将被赋予更高的权重) 第七步:创建另一个模型并对数据集进行预测(此模型尝试更正先前模型中的错误)。 第八步:类似地,创建多个模型,每个模型校正先前模型的错误。 第九步:最终模型(强学习器)是所有模型(弱学习器)的加权平均值。

Boosting算法: AdaBoost, GBM, XGBM, Light GBM, CatBoost

Bagging meta-estimator是一种集成算法,可用于分类(BaggingClassifier)和回归(BaggingRegressor)问题。它采用典型的bagging技术进行预测。以下是Bagging meta-estimator算法的步骤:

第一步:从原始数据集(Bootstrapping)创建随机子集。

第二步:数据集的子集包括所有特征。

第三步:用户指定的基础估计器在这些较小的集合上拟合。

第四步:将每个模型的预测结合起来得到最终结果。

4.5.2 Boosting

自适应增强或AdaBoost

是最简单的boosting算法之一。通常用决策树来建模。创建多个顺序模型,每个模型都校正上一个模型的错误。AdaBoost为错误预测的观测值分配权重,后续模型来正确预测这些值。

以下是执行AdaBoost算法的步骤:

第一步:最初,数据集中的所有观察值都具有相同的权重。

第二步:在数据子集上建立一个模型。

第三步:使用此模型,可以对整个数据集进行预测。

第四步:通过比较预测值和实际值来计算误差。

第五步:在创建下一个模型时,会给预测错误的数据点赋予更高的权重。

第六步:可以使用误差值确定权重。例如,误差越大,分配给观察值的权重越大。

第七步:重复该过程直到误差函数没有改变,或达到估计器数量的最大限制

Gradient Boosting(梯度提升GBM)

Gradient Boosting或GBM是另一种集成机器学习算法,适用于回归和分类问题。GBM使用boosting技术,结合了许多弱学习器,以形成一个强大的学习器。回归树用作基础学习器,每个后续的树都是基于前一棵树计算的错误构建的。

我们将使用一个简单的例子来理解GBM算法。我们会使用以下数据预测一群人的年龄:

第一步:假设平均年龄是数据集中所有观测值的预测值。

第二步:使用该平均预测值和年龄的实际值来计算误差:

第三步:使用上面计算的误差作为目标变量创建树模型。我们的目标是找到最佳分割以最小化误差。

第四步:该模型的预测与预测1相结合:

第五步:上面计算的这个值是新的预测。

第六步:使用此预测值和实际值计算新误差:

第七步:重复步骤2到6,直到最大迭代次数(或误差函数不再改变)

XGBoost

XGBoost(extreme Gradient Boosting)是梯度提升算法的高级实现。实践证明,XGBoost是一种高效的ML算法,广泛应用于机器学习竞赛和黑客马拉松。 XGBoost具有很高的预测能力,几乎比其他梯度提升技术快10倍。它还包括各种正规化,可减少过拟合并提高整体性能。因此,它也被称为“正则化提升”技术。

让我们看看XGBoost为何比其他技术更好:

正则化:

标准GBM实现没有像XGBoost那样的正则化,因此,XGBoost还有助于减少过拟合

并行处理:

XGBoost实现并行处理,并且比GBM更快

XGBoost还支持Hadoop上的实现

高灵活性:

XGBoost允许用户自定义优化目标和评估标准,为模型添加全新维度

处理缺失值:

XGBoost有一个内置的例程来处理缺失值

树剪枝:

XGBoost先进行分割,直到指定的max_depth,然后开始向后修剪树并删除没有正向增益的分割

内置交叉验证:

XGBoost允许用户在提升过程的每次迭代中运行交叉验证,因此很容易在一次运行中获得精确的最佳提升迭代次数

Light GBM

在讨论Light GBM如何工作之前,先理解为什么在我们有如此多其他算法时(例如我们上面看到的算法)我们还需要这个算法。当数据集非常大时,Light GBM会击败所有其他算法。与其他算法相比,Light GBM在较大的数据集上运行所需的时间较短。

LightGBM是一个梯度提升框架,它使用基于树的算法并遵循逐叶子的方式(leaf-wise),而其他算法以逐层级(level-wise)模式工作。下图帮助你更好地理解二者差异:

逐叶子方式可能在较小的数据集上导致过拟合,但可以通过使用’max_depth’参数来避免这种情况。你可以在本文中阅读有关Light GBM及其与XGB比较的更多信息。

import lightgbm as lgb
train_data=lgb.Dataset(x_train,label=y_train)
#define parameters
params = {'learning_rate':0.001}
model= lgb.train(params, train_data, 100)
y_pred=model.predict(x_test)
for i in range(0,185):
    if y_pred[i]>=0.5:
        y_pred[i]=1
else:
    y_pred[i]=0
#回归问题示例代码:
import lightgbm as lgb
train_data=lgb.Dataset(x_train,label=y_train)
params = {'learning_rate':0.001}
model= lgb.train(params, train_data, 100)
from sklearn.metrics import mean_squared_error
rmse=mean_squared_error(y_pred,y_test)**0.5

CatBoost

处理类别型变量是一个繁琐的过程,尤其是你有大量此类变量时。当你的类别变量有很多标签(即它们是高度基数)时,对它们执行one-hot编码会指数级的增加维度,会让数据集的使用变得非常困难。

第三步:使用此模型,可以对整个数据集进行预测。

第四步:通过比较预测值和实际值来计算误差。

第五步:在创建下一个模型时,会给预测错误的数据点赋予更高的权重。

第六步:可以使用误差值确定权重。例如,误差越大,分配给观察值的权重越大。

第七步:重复该过程直到误差函数没有改变,或达到估计器数量的最大限制

Gradient Boosting(梯度提升GBM)

Gradient Boosting或GBM是另一种集成机器学习算法,适用于回归和分类问题。GBM使用boosting技术,结合了许多弱学习器,以形成一个强大的学习器。回归树用作基础学习器,每个后续的树都是基于前一棵树计算的错误构建的。

我们将使用一个简单的例子来理解GBM算法。我们会使用以下数据预测一群人的年龄:

第一步:假设平均年龄是数据集中所有观测值的预测值。

第二步:使用该平均预测值和年龄的实际值来计算误差:

第三步:使用上面计算的误差作为目标变量创建树模型。我们的目标是找到最佳分割以最小化误差。

第四步:该模型的预测与预测1相结合:

第五步:上面计算的这个值是新的预测。

第六步:使用此预测值和实际值计算新误差:

第七步:重复步骤2到6,直到最大迭代次数(或误差函数不再改变)

XGBoost

XGBoost(extreme Gradient Boosting)是梯度提升算法的高级实现。实践证明,XGBoost是一种高效的ML算法,广泛应用于机器学习竞赛和黑客马拉松。 XGBoost具有很高的预测能力,几乎比其他梯度提升技术快10倍。它还包括各种正规化,可减少过拟合并提高整体性能。因此,它也被称为“正则化提升”技术。

让我们看看XGBoost为何比其他技术更好:

正则化:

标准GBM实现没有像XGBoost那样的正则化,因此,XGBoost还有助于减少过拟合

并行处理:

XGBoost实现并行处理,并且比GBM更快

XGBoost还支持Hadoop上的实现

高灵活性:

XGBoost允许用户自定义优化目标和评估标准,为模型添加全新维度

处理缺失值:

XGBoost有一个内置的例程来处理缺失值

树剪枝:

XGBoost先进行分割,直到指定的max_depth,然后开始向后修剪树并删除没有正向增益的分割

内置交叉验证:

XGBoost允许用户在提升过程的每次迭代中运行交叉验证,因此很容易在一次运行中获得精确的最佳提升迭代次数

Light GBM

在讨论Light GBM如何工作之前,先理解为什么在我们有如此多其他算法时(例如我们上面看到的算法)我们还需要这个算法。当数据集非常大时,Light GBM会击败所有其他算法。与其他算法相比,Light GBM在较大的数据集上运行所需的时间较短。

LightGBM是一个梯度提升框架,它使用基于树的算法并遵循逐叶子的方式(leaf-wise),而其他算法以逐层级(level-wise)模式工作。下图帮助你更好地理解二者差异:

逐叶子方式可能在较小的数据集上导致过拟合,但可以通过使用’max_depth’参数来避免这种情况。你可以在本文中阅读有关Light GBM及其与XGB比较的更多信息。

import lightgbm as lgb
train_data=lgb.Dataset(x_train,label=y_train)
#define parameters
params = {'learning_rate':0.001}
model= lgb.train(params, train_data, 100)
y_pred=model.predict(x_test)
for i in range(0,185):
    if y_pred[i]>=0.5:
        y_pred[i]=1
else:
    y_pred[i]=0
#回归问题示例代码:
import lightgbm as lgb
train_data=lgb.Dataset(x_train,label=y_train)
params = {'learning_rate':0.001}
model= lgb.train(params, train_data, 100)
from sklearn.metrics import mean_squared_error
rmse=mean_squared_error(y_pred,y_test)**0.5

CatBoost

处理类别型变量是一个繁琐的过程,尤其是你有大量此类变量时。当你的类别变量有很多标签(即它们是高度基数)时,对它们执行one-hot编码会指数级的增加维度,会让数据集的使用变得非常困难。

CatBoost可以自动处理类别型变量,并且不需要像其他机器学习算法那样进行大量数据预处理。这篇文章详细解释了CatBoost。

     

如果觉得小墨的文章对您有用,请随意打赏。您的支持将鼓励我继续创作!

3条评论

发表评论