声明：版权所有，转载请联系作者并注明出处  http://blog.csdn.net/u013719780?viewmode=contents

博主简介：风雪夜归子（英文名：Allen），机器学习算法攻城狮，喜爱钻研Meachine Learning的黑科技，对Deep Learning和Artificial Intelligence充满兴趣，经常关注Kaggle数据挖掘竞赛平台，对数据、Machine Learning和Artificial Intelligence有兴趣的童鞋可以一起探讨哦，个人CSDN博客：http://blog.csdn.net/u013719780?viewmode=contents

分类模型处理表示类别的离散变量，而回归模型则处理可以取任意实数的目标变量。但是二者基本的原则类似，都是通过确定一个模型，将输入特征映射到预测的输出。回归模型和分类模型都是监督学习的一种形式。

回归模型可以用在如下一些场景：

 预测股票收益和其他经济相关的因素；

 预测贷款违约造成的损失（可以和分类模型相结合，分类模型预测违约概率，回归模型预测违约损失）；

 推荐系统（博文 spark机器学习笔记：（三）用Spark Python构建推荐系统 中的交替最小二乘分解模型在每次迭代时都使用了线性回归）；

 基于用户的行为和消费模式，预测顾客对于零售、移动或者其他商业形态的存在价值。

本文主要内容如下：

 介绍MLlib中的各种回归模型；

 讨论回归模型的特征提取和目标变量的变换；

 使用MLlib训练回归模型；

 介绍如何用训练好的模型做预测；

 使用交叉验证研究设置不同的参数对性能的影响。

1 特征提取

因为回归的基础模型和分类模型一样，所以我们可以使用同样的方法来处理输入的特征。实际中唯一的不同是，回归模型的预测目标是实数变量，而分类模型的预测目标是类别编号。为了满足两种情况，MLlib中的LabeledPoint类已经考虑了这一点，类中的label字段使用Double类型。本文选择bike sharing数据集做实验。这个数据集记录了bike

sharing系统每小时自行车的出租次数。另外还包括日期、时间、天气、季节和节假日等相关信息（bike sharing数据集下载地址：http://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset）。各个字段的含义如下：

- instant: record index

- dteday : date

- season : season (1:springer, 2:summer, 3:fall, 4:winter)

- yr : year (0: 2011, 1:2012)

- mnth : month ( 1 to 12)

- hr : hour (0 to 23)

- holiday : weather day is holiday or not (extracted from http://dchr.dc.gov/page/holiday-schedule)

- weekday : day of the week

- workingday : if day is neither weekend nor holiday is 1, otherwise is 0.

+ weathersit : 

- 1: Clear, Few clouds, Partly cloudy, Partly cloudy

- 2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist

- 3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds

- 4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog

- temp : Normalized temperature in Celsius. The values are divided to 41 (max)

- atemp: Normalized feeling temperature in Celsius. The values are divided to 50 (max)

- hum: Normalized humidity. The values are divided to 100 (max)

- windspeed: Normalized wind speed. The values are divided to 67 (max)

- casual: count of casual users

- registered: count of registered users

- cnt: count of total rental bikes including both casual and registered

输入命令 sed 1d hour.csv > hour_noheader.csv，将文件hour.csv的第一行头文件去掉并保存到新的 hour_noheader.csv。

import numpy as np
path = '/Users/youwei.tan/Desktop/Bike-Sharing-Dataset/hour_noheader.csv'
raw_data = sc.textFile(path) 
num_data = raw_data.count() 
records =raw_data.map(lambda x: x.split(',')) 
first = records.first() 
print '数据的第一行:',first 
print '数据样本数:',num_data

输出结果：

数据的第一行: [u'1', u'2011-01-01', u'1', u'0', u'1', u'0', u'0', u'6', u'0', u'1', u'0.24', u'0.2879', u'0.81', u'0', u'3', u'13', u'16']

数据样本数: 17379

结果显示，数据集中共有17 379个小时的记录。接下来的实验，我们会忽略记录中的instant和dteday 。忽略两个记录次数的变量casual 和registered ， 只保留cnt （ casual 和

registered的和）。最后就剩下12个变量，其中前8个是类型变量，后4个是归一化后的实数变量。对其中8个类型变量，我们使用之前提到的二元编码，剩下4个实数变量不做处理。

因为变量records下文经常要用到，此处对其进行缓存：

records.cache()

为了将类型特征表示成二维形式，我们将特征值映射到二元向量中非0的位置。下面定义这样一个映射函数：

def get_mapping(rdd, idx):
    return rdd.map(lambda fields: fields[idx]).distinct().zipWithIndex().collectAsMap()

上面的函数首先将第idx列的特征值去重，然后对每个值使用zipWithIndex函数映射到一个唯一的索引，这样就组成了一个RDD的键值映射，键是变量，值是索引。上述索引便是特征在二元向量中对应的非0位置，最后我们将这个RDD表示成Python的字典类型。

下面，我们用特征矩阵的第三列（索引2）来测试上面的映射函数：

print '第三个特征的类别编码： %s '%get_mapping(records,2) 

输出结果：

第三个特征的类别编码： {u'1': 0, u'3': 1, u'2': 2, u'4': 3} 

接着，对是类型变量的列（第2~9列）应用该函数：

mappings = [get_mapping(records, i) for i in range(2,10)]   #对类型变量的列（第2~9列）应用映射函数
print '类别特征打编码字典:',mappings
cat_len = sum(map(len,[i for i in mappings]))        #类别特征的个数                 
#cat_len = sum(map(len, mappings))
#print map(len,mappings)

num_len = len(records.first()[11:15])                      #数值特征的个数
total_len = num_len+cat_len                                  #所有特征的个数
print '类别特征的个数： %d'% cat_len
print '数值特征的个数： %d'% num_len
print '所有特征的个数:：%d' % total_len
          
          

输出结果：

类别特征的编码字典: [{u'1': 0, u'3': 1, u'2': 2, u'4': 3}, {u'1': 0, u'0': 1}, {u'11': 0, u'10': 1, u'12': 2, u'1': 3, u'3': 4, u'2': 5, u'5': 6, u'4': 7, u'7': 8, u'6': 9, u'9': 10, u'8': 11}, {u'20': 0, u'21': 1, u'22': 2, u'23': 3, u'1': 4, u'0': 5, u'3': 6, u'2': 7, u'5': 8, u'4': 9, u'7': 10, u'6': 11, u'9': 12, u'8': 13, u'11': 14, u'10': 15, u'13': 16, u'12': 17, u'15': 18, u'14': 19, u'17': 20, u'16': 21, u'19': 22, u'18': 23}, {u'1': 0, u'0': 1}, {u'1': 0, u'0': 1, u'3': 2, u'2': 3, u'5': 4, u'4': 5, u'6': 6}, {u'1': 0, u'0': 1}, {u'1': 0, u'3': 1, u'2': 2, u'4': 3}]

类别特征的个数： 57

数值特征的个数： 4

所有特征的个数:：61

1.1 为线性模型创建特征向量

接下来用上面的映射函数将所有类型特征转换为二元编码的特征。为了方便对每条记录提取特征和标签，我们分别定义两个辅助函数extract_features和extract_label。如下为代码实现，注意需要引入numpy和MLlib的LabeledPoint对特征向量和目标变量进行封装：

from pyspark.mllib.regression import LabeledPoint
import numpy as np

def extract_features(record): 
    cat_vec = np.zeros(cat_len) 
    step = 0
    for i,raw_feature in enumerate(record[2:9]):
        dict_code = mappings[i]
        index = dict_code[raw_feature]
        cat_vec[index+step] = 1
        step = step+len(dict_code)
    num_vec = np.array([float(raw_feature) for raw_feature in record[10:14]])
    return np.concatenate((cat_vec, num_vec))

def extract_label(record):
    return float(record[-1])

在extract_features函数中，我们遍历了数据的每一行每一列，根据已经创建的映射对每个特征进行二元编码。其中step变量用来确保非0特征在整个特征向量中位于正确的位置（另外一种实现方法是将若干较短的二元向量拼接在一起）。数值向量直接对之前已经被转换成浮点数的数据用numpy的array进行封装。最后将二元向量和数值向量拼接起来。定义extract_label函数将数据中的最后一列cnt的数据转换成浮点数。

下面对数据进行特征提取：

data = records.map(lambda point: LabeledPoint(extract_label(point),extract_features(point)))
first_point = data.first()

print '原始特征向量:' +str(first[2:])
print '标签:' + str(first_point.label)
print '对类别特征进行独热编码之后的特征向量: \n' + str(first_point.features)
print '对类别特征进行独热编码之后的特征向量长度:' + str(len(first_point.features))

输出结果：

原始特征向量:[u'1', u'0', u'1', u'0', u'0', u'6', u'0', u'1', u'0.24', u'0.2879', u'0.81', u'0', u'3', u'13', u'16']

标签:16.0

对类别特征进行独热编码之后的特征向量: 

[1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.24,0.2879,0.81,0.0]

对类别特征进行独热编码之后的特征向量长度:61

1.2 为决策树创建特征向量

我们已经知道，决策树模型可以直接使用原始数据（不需要将类型数据用二元向量表示）。因此，只需要创建一个分割函数简单地将所有数值转换为浮点数，最后用numpy的array封装：

def extract_features_dt(record):
    return np.array(map(float, record[2:14]))
data_dt = records.map(lambda point: LabeledPoint(extract_label(point), extract_features_dt(point)))
first_point_dt = data_dt.first()
print '决策树特征向量: '+str(first_point_dt.features)
print '决策树特征向量长度: '+str(len(first_point_dt.features))

输出结果：

决策树特征向量: [1.0,0.0,1.0,0.0,0.0,6.0,0.0,1.0,0.24,0.2879,0.81,0.0]

决策树特征向量长度: 12

2 回归模型的训练和应用

使用决策树和线性模型训练回归模型的步骤和使用分类模型相同，都是简单将训练数据封装在LabeledPoint的RDD中，并送到相关的train方法上进行训练。注意在Scala中，如果要自定义不同的模型参数（比如SGD优化的正则化和步长），就需要初始化一个新的模型实例，使用实例的optimizer变量访问和设置参数。

Python提供了方便我们访问所有模型参数的方法，因此只要使用相关方法即可。可以通过引入相关模块，并调用train方法中的help函数查看这些方法的具体细节：

from pyspark.mllib.regression import LinearRegressionWithSGD
from pyspark.mllib.tree import DecisionTree
help(LinearRegressionWithSGD.train)

输出结果：

Help on method train in module pyspark.mllib.regression:

train(cls, data, iterations=100, step=1.0, miniBatchFraction=1.0, initialWeights=None, regParam=0.0, regType=None, intercept=False, validateData=True, convergenceTol=0.001) method of __builtin__.type instance

    Train a linear regression model using Stochastic Gradient

    Descent (SGD).

    This solves the least squares regression formulation

        f(weights) = 1/(2n) ||A weights - y||^2,

    which is the mean squared error.

    Here the data matrix has n rows, and the input RDD holds the

    set of rows of A, each with its corresponding right hand side

    label y. See also the documentation for the precise formulation.

    :param data:              The training data, an RDD of

                              LabeledPoint.

    :param iterations:        The number of iterations

                              (default: 100).

    :param step:              The step parameter used in SGD

                              (default: 1.0).

    :param miniBatchFraction: Fraction of data to be used for each

                              SGD iteration (default: 1.0).

    :param initialWeights:    The initial weights (default: None).

    :param regParam:          The regularizer parameter

                              (default: 0.0).

    :param regType:           The type of regularizer used for

                              training our model.

                              :Allowed values:

                                 - "l1" for using L1 regularization (lasso),

                                 - "l2" for using L2 regularization (ridge),

                                 - None for no regularization

                                 (default: None)

    :param intercept:         Boolean parameter which indicates the

                              use or not of the augmented representation

                              for training data (i.e. whether bias

                              features are activated or not,

                              default: False).

    :param validateData:      Boolean parameter which indicates if

                              the algorithm should validate data

                              before training. (default: True)

    :param convergenceTol:    A condition which decides iteration termination.

                              (default: 0.001)

    .. versionadded:: 0.9.0

调用决策树模型的trainRegressor方法查看帮助信息

help(DecisionTree.trainRegressor)

输出结果：

Help on method trainRegressor in module pyspark.mllib.tree:

trainRegressor(cls, data, categoricalFeaturesInfo, impurity='variance', maxDepth=5, maxBins=32, minInstancesPerNode=1, minInfoGain=0.0) method of __builtin__.type instance

    Train a DecisionTreeModel for regression.

    :param data: Training data: RDD of LabeledPoint.

                 Labels are real numbers.

    :param categoricalFeaturesInfo: Map from categorical feature

             index to number of categories.

             Any feature not in this map is treated as continuous.

    :param impurity: Supported values: "variance"

    :param maxDepth: Max depth of tree.

             E.g., depth 0 means 1 leaf node.

             Depth 1 means 1 internal node + 2 leaf nodes.

    :param maxBins: Number of bins used for finding splits at each

             node.

    :param minInstancesPerNode: Min number of instances required at

             child nodes to create the parent split

    :param minInfoGain: Min info gain required to create a split

    :return: DecisionTreeModel

    Example usage:

    >>> from pyspark.mllib.regression import LabeledPoint

    >>> from pyspark.mllib.tree import DecisionTree

    >>> from pyspark.mllib.linalg import SparseVector

    >>>

    >>> sparse_data = [

    ...     LabeledPoint(0.0, SparseVector(2, {0: 0.0})),

    ...     LabeledPoint(1.0, SparseVector(2, {1: 1.0})),

    ...     LabeledPoint(0.0, SparseVector(2, {0: 0.0})),

    ...     LabeledPoint(1.0, SparseVector(2, {1: 2.0}))

    ... ]

    >>>

    >>> model = DecisionTree.trainRegressor(sc.parallelize(sparse_data), {})

    >>> model.predict(SparseVector(2, {1: 1.0}))

    1.0

    >>> model.predict(SparseVector(2, {1: 0.0}))

    0.0

    >>> rdd = sc.parallelize([[0.0, 1.0], [0.0, 0.0]])

    >>> model.predict(rdd).collect()

    [1.0, 0.0]

    .. versionadded:: 1.1.0

在bike sharing 数据上训练回归模型

我们已经从bike sharing数据中提取了用于训练模型的特征，下面进行具体的训练。首先训练线性模型并测试该模型在训练数据上的预测效果：

linear_model = LinearRegressionWithSGD.train(data, iterations=10, step=0.1, intercept =False)
true_vs_predicted = data.map(lambda point:(point.label,linear_model.predict(point.features)))
print '线性回归模型对前5个样本的预测值: '+ str(true_vs_predicted.take(5))

输出结果：

线性回归模型对前5个样本的预测值: [(16.0, 117.89250386724844), (40.0, 116.22496123192109), (32.0, 116.02369145779232), (13.0, 115.67088016754431), (1.0, 115.56315650834314)]

上述代码中我们没有使用默认的迭代次数和步长，而是使用较小的迭代次数以缩短训练时间，关于步长的设置我们稍后会详细介绍。

接下来，我们在trainRegressor中使用默认参数来训练决策树模型（相当于深度为5的树）。注意，这里训练数据集是从原始特征中提取的，名为data_dt（不同于之前线性模型中使用的二元编码的特征）。另外，我们还需要为categoricalFeaturesInfo传入一个字典参数，这个字典参数将类型特征的索引映射到特征中类型的数目。如果某个特征值不在这个字典中，则将其映射设置为空：

dt_model = DecisionTree.trainRegressor(data_dt,{})
preds = dt_model.predict(data_dt.map(lambda p: p.features))
actual = data.map(lambda p:p.label)
true_vs_predicted_dt = actual.zip(preds)
print '决策树回归模型对前5个样本的预测值: '+str(true_vs_predicted_dt.take(5)) 
print '决策树模型的深度: ' + str(dt_model.depth()) 
print '决策树模型的叶子节点个数: '+str(dt_model.numNodes())

输出结果：

决策树回归模型对前5个样本的预测值: [(16.0, 54.913223140495866), (40.0, 54.913223140495866), (32.0, 53.171052631578945), (13.0, 14.284023668639053), (1.0, 14.284023668639053)]

决策树模型的深度: 5

决策树模型的叶子节点个数: 63

3 评估回归模型的性能

上一篇博文评估分类模型仅仅关注预测输出的类别和实际类别。特别是对于所有预测的二元结果，某个样本预测的正确与否并不重要，我们更关心预测结果中正确或者错误的总数。对回归模型而言，因为目标变量是任一实数，所以我们的模型不大可能精确预测到目标变量。然而，我们可以计算预测值和实际值的误差，并用某种度量方式进行评估。

一些用于评估回归模型的方法包括：均方误差（MSE，Mean Squared Error）、均方根误差（RMSE，Root Mean Squared Error）、平均绝对误差（MAE，Mean Absolute Error）、R-平方系数（R-squared coefficient）等。

3.1 均方误差和均方根误差

MSE是平方误差的均值，用作最小二乘回归的损失函数，公式如下：

这个公式计算的是所有样本预测值和实际值平方差之和，最后除以样本总数。而RMSE是MSE的平方根。MSE的公式类似平方损失函数，会进一步放大误差。

为了计算模型预测的平均误差，我们首先预测RDD实例LabeledPoint中每个特征向量，然后计算预测值与实际值的误差并组成一个Double数组的RDD，最后使用mean方法计算所有Double值的平均值。计算平方误差函数实现如下：

def squared_error(actual, pred): 
    return (pred-actual)**2 

3.2 平均绝对误差

MAE是预测值和实际值的差的绝对值的平均值。公式如下：

MAE和MSE大体类似，区别在于MAE对大的误差没有惩罚。计算MAE的代码如下：

def abs_error(actual, pred): 
    return np.abs(pred-actual) 

3.3 均方根对数误差

这个度量方法虽然没有MSE和MAE使用得广，但被用于Kaggle中以bike sharing作为数据集的比赛。RMSLE可以认为是对预测值和目标值进行对数变换后的RMSE。这个度量方法适用于目标变量值域很大，并且没有必要对预测值和目标值的误差进行惩罚的情况。另外，它也适用于计算误差的百分率而不是误差的绝对值。计算RMSLE的代码如下：

def squared_log_error(pred, actual): 
    return (np.log(pred+1)-np.log(actual+1))**2

3.4 R-平方系数

R-平方系数，也称判定系数，用来评估模型拟合数据的好坏，常用于统计学中。R-平方系数具体测量目标变量的变异度（degree of variation），最终结果为0到1的一个值，1表示模型能够完美拟合数据。

3.5 计算不同度量下的性能

根据上面定义的函数，我们在bike sharing数据集上计算不同度量下的性能。

3.5.1 线性模型

我们的方法对RDD的每一条记录应用相关的误差函数，其中线性模型的误差函数为true_vs_predicted，相关代码实现如下:

mse = true_vs_predicted.map(lambda (t, p): squared_error(t, p)).mean() 
mae = true_vs_predicted.map(lambda (t, p): abs_error(t, p)).mean() 
rmsle = np.sqrt(true_vs_predicted.map(lambda (t, p): squared_log_error(t, p)).mean()) 
print 'Linear Model - Mean Squared Error: %2.4f' % mse 
print 'Linear Model - Mean Absolute Error: %2.4f' % mae 
print 'Linear Model - Root Mean Squared Log Error: %2.4f' % rmsle

输出结果：

Linear Model - Mean Squared Error: 30679.4539

Linear Model - Mean Absolute Error: 130.6429

Linear Model - Root Mean Squared Log Error: 1.4653

3.5.2 决策树

决策树的误差函数为true_vs_predicted_dt，相关代码如下:

mse_dt = true_vs_predicted_dt.map(lambda (t, p): squared_error(t, p)).mean()
mae_dt = true_vs_predicted_dt.map(lambda (t, p): abs_error(t, p)).mean()
rmsle_dt = np.sqrt(true_vs_predicted_dt.map(lambda (t, p): squared_log_error(t, p)).mean())
print 'Decision Tree - Mean Squared Error: %2.4f' % mse_dt
print 'Decision Tree - Mean Absolute Error: %2.4f' % mae_dt
print 'Decision Tree - Root Mean Squared Log Error: %2.4f' %rmsle_dt

输出结果：

Decision Tree - Mean Squared Error: 11560.7978

Decision Tree - Mean Absolute Error: 71.0969

Decision Tree - Root Mean Squared Log Error: 0.6259

%matplotlib inline
import matplotlib
from matplotlib.pyplot import hist
import numpy as np 
import matplotlib.pyplot as plt 

targets = records.map(lambda r: float(r[-1])).collect() 
hist(targets, bins=40, color='lightblue', normed=True) 
fig = matplotlib.pyplot.gcf() 
fig.set_size_inches(12, 6)
plt.show()

log_targets = records.map(lambda r : np.log(float(r[-1]))).collect() 
plt.hist(log_targets, bins = 40, color ='lightblue', normed =True) 
fig = plt.gcf()
fig.set_size_inches(12,6)
plt.show()

sqrt_targets = records.map(lambda r: np.sqrt(float(r[-1]))).collect() 
plt.hist(sqrt_targets, bins=40, color='lightblue', normed=True) 
fig = matplotlib.pyplot.gcf() 
fig.set_size_inches(12, 6) 

data_log = data.map(lambda lp:LabeledPoint(np.log(lp.label),lp.features)) 
model_log =LinearRegressionWithSGD.train(data_log, iterations=10, step=0.1) 
true_vs_predicted_log = data_log.map(lambda p:(np.exp(p.label),np.exp(model_log.predict(p.features))))

#计算模型的MSE,MAE,RMSLE
mse_log = true_vs_predicted_log.map(lambda(t, p): squared_error(t,p)).mean() 
mae_log = true_vs_predicted_log.map(lambda(t, p): abs_error(t, p)).mean() 
rmsle_log = np.sqrt(true_vs_predicted_log.map(lambda(t, p): squared_log_error(t, p)).mean())

print'Linear Model —— Mean Squared Error:%2.4f'% mse_log 
print'Linear Model —— Mean Absolue Error:%2.4f'% mae_log 
print'Linear Model —— Root Mean Squared Log Error:%2.4f'% rmsle_log 
print'Linear Model —— Non log-transformed predictions:\n'+ str(true_vs_predicted.take(3))
print'Linear Model —— Log-transformed predictions:\n'+ str(true_vs_predicted_log.take(3))

data_dt_log = data_dt.map(lambda lp:LabeledPoint(np.log(lp.label), lp.features)) 
dt_model_log = DecisionTree.trainRegressor(data_dt_log,{}) 
preds_log = dt_model_log.predict(data_dt_log.map(lambda p:p.features)) 
actual_log = data_dt_log.map(lambda p: p.label) 
true_vs_predicted_dt_log = actual_log.zip(preds_log).map(lambda(t,p):(np.exp(t), np.exp(p)))

#计算模型的MSE,MAE,RMSLE
mse_log_dt = true_vs_predicted_dt_log.map(lambda(t, p): squared_error(t, p)).mean() 
mae_log_dt = true_vs_predicted_dt_log.map(lambda(t, p): abs_error(t,p)).mean() 
rmsle_log_dt = np.sqrt(true_vs_predicted_dt_log.map(lambda(t, p): 
squared_log_error(t, p)).mean())

print'Decision Tree —— Mean Squared Error:%2.4f'% mse_log_dt 
print'Decision Tree —— Mean Absolue Error:%2.4f'% mae_log_dt 
print'Decision Tree —— Root Mean Squared Log Error:%2.4f'% rmsle_log_dt 
print'Decision Tree —— Non log-transformed predictions:\n'+ str(true_vs_predicted_dt.take(3))
print'Decision Tree —— Log-transformed predictions:\n'+str(true_vs_predicted_dt_log.take(3))

data_with_idx = data.zipWithIndex().map(lambda(k,v):(v,k)) 
test = data_with_idx.sample(False,0.2,42) 
train = data_with_idx.subtractByKey(test)

train_test_data_split = data.randomSplit([0.8,0.2],123) 
train = train_test_data_split[0] 
test = train_test_data_split[1]
print '测试集的样本数:',test.count()
print '训练集的样本数:',train.count()

train_test_data_dt_split = data_dt.randomSplit([0.8,0.2],123) 
train_dt = train_test_data_dt_split[0] 
test_dt = train_test_data_dt_split[1]

def evaluate(train, test, iterations, step, regParam, regType, intercept):
    model =LinearRegressionWithSGD.train(train, iterations, step, regParam=regParam, 
                                         regType=regType,intercept=intercept)
    testLabel_vs_testPrediction = test.map(lambda point:(point.label, model.predict(point.features)))
    rmsle = np.sqrt(testLabel_vs_testPrediction.map(lambda(t,p):squared_log_error(t,p)).mean())
    return rmsle

params = [1, 5, 10, 20, 50, 100]
metrics = [evaluate(train, test, param, 0.01, 0.0, 'l2', False) for param in params]
print params
print metrics

plt.plot(params, metrics) 
fig = matplotlib.pyplot.gcf() 
fig.set_size_inches(12, 6) 
plt.xscale('log')
plt.show()

params=[0.01,0.025,0.05,0.1,0.5,1.0] 
metrics =[evaluate(train, test,10,param,0.0,'l2',False)for param in params]
for i in range(len(params)):
    print'the rmsle:%f when step :%f'%(metrics[i],params[i])

#绘制步长与RMSLE的关系图：
plt.plot(params, metrics) 
fig = matplotlib.pyplot.gcf() 
fig.set_size_inches(12, 6) 
plt.xscale('log')
plt.xlabel('step')
plt.ylabel('RMSLE')
plt.show()

params=[0.0,0.01,0.1,1.0,5.0,10.0,20.0]
metrics =[evaluate(train, test,10,0.1, param,'l2',False) for param in params]
for i in range(len(params)):
    print'the rmsle:%f when regParam :%f'%(metrics[i],params[i])

#绘制L2正则化系数与RMSLE的关系图：
plt.plot(params, metrics)
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(12, 8)
plt.xscale('log')
plt.xlabel('regParam')
plt.ylabel('RMSLE')
plt.show()

params=[0.0,0.01,0.1,1.0,10.0,100.0,1000.0]
metrics =[evaluate(train, test,10,0.1, param,'l1',False) for param in params]
for i in range(len(params)):
    print'the rmsle:%f when regParam :%f'%(metrics[i],params[i])

#绘制L2正则化系数与RMSLE的关系图：
plt.plot(params, metrics)
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(12, 8)
plt.xscale('log')
plt.xlabel('regParam')
plt.ylabel('RMSLE')
plt.show()

model_l1 = LinearRegressionWithSGD.train(train,10,0.1,regParam=1.0, regType='l1', intercept=False) 
model_l2 = LinearRegressionWithSGD.train(train,10,0.1,regParam=1.0, regType='l2', intercept=False) 
model_l1_10 = LinearRegressionWithSGD.train(train,10,0.1,regParam=10.0, regType='l1', intercept=False) 
model_l2_10 = LinearRegressionWithSGD.train(train,10,0.1,regParam=10.0, regType='l2', intercept=False) 
model_l1_100 = LinearRegressionWithSGD.train(train,10,0.1,regParam=100.0, regType='l1', intercept=False) 
model_l2_100 = LinearRegressionWithSGD.train(train,10,0.1,regParam=100.0, regType='l2', intercept=False)

print'L1 (1.0) number of zero weights:'+ str(sum(model_l1.weights.array == 0))
print'L2 (1.0) number of zero weights:'+ str(sum(model_l2.weights.array == 0))
print'L1 (10.0) number of zeros weights:'+ str(sum(model_l1_10.weights.array == 0))
print'L2 (10.0) number of zeros weights:'+ str(sum(model_l2_10.weights.array == 0))
print'L1 (100.0) number of zeros weights:'+ str(sum(model_l1_100.weights.array == 0))
print'L2 (100.0) number of zeros weights:'+ str(sum(model_l2_100.weights.array == 0))

params=[False, True]
metrics =[evaluate(train, test, 10, 0.1, 1.0,'l2', param) for param in params]
for i in range(len(params)):
    print'the rmsle:%f when intercept:%f'%(metrics[i],params[i])

#绘制L2正则化系数与RMSLE的关系图：
plt.bar(params, metrics, color='r')
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(12, 8)
plt.xlabel('intercept')
plt.ylabel('RMSLE')
plt.show()

def evaluate_dt(train, test, maxDepth, maxBins): 
    model =DecisionTree.trainRegressor(train,{},impurity='variance', maxDepth=maxDepth, maxBins=maxBins)
    predictions = model.predict(test.map(lambda point: point.features)) 
    actual = test.map(lambda point: point.label) 
    actual_vs_predictions = actual.zip(predictions) 
    rmsle = np.sqrt(actual_vs_predictions.map(lambda(t, p): squared_log_error(t,p)).mean())
    return rmsle

params=[1,2,3,4,5,10,20] 
metrics =[evaluate_dt(train_dt, test_dt, param,32) for param in params]
for i in range(len(params)):    
    print'the rmsle:%f when maxDepth :%d'%(metrics[i],params[i])

#绘制树的最大深度与RMSLE的关系图：
plt.plot(params, metrics)
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(12, 8)
plt.xlabel('maxDepth')
plt.ylabel('RMSLE')
plt.show()

params=[2,4,8,16,32,64,100] 
metrics =[evaluate_dt(train_dt, test_dt,5, param) for param in params]
for i in range(len(params)):    
    print'the rmsle:%f when maxBins :%d'%(metrics[i],params[i]) 

#绘制树的最大划分数与RMSLE的关系图：
plt.plot(params, metrics)
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(12, 8)
plt.xlabel('maxDepth')
plt.ylabel('RMSLE')
plt.show()