回归算法

回归算法是一种有监督算法。

回归算法是一种比较常用的机器学习算法，用来建立“解释”变量(自变量X)和观 测值(因变量Y)之间的关系；从机器学习的角度来讲，用于构建一个算法模型(函 数)来做属性(X)与标签(Y)之间的映射关系，在算法的学习过程中，试图寻找一个 函数，使得参数之间的关系拟合性最好。

回归算法中算法(函数)的最终结果是一个连续的数据值，输入值(属性值)是一个d 维度的属性/数值向量。

线性回归算法

损失函数及求解

 假设有以下线性函数

最终要求是计算出 θ 的值，并选择最优的 θ 值构成算法公式。

θ求解方式：最小二乘法(直接计算，目标函数是平方和损失函数)、梯度下降 (BGD\SGD\MBGD)

最大似然估计

 

1.似然函数

 

2.对数似然、目标函数

 

3.θ的求解过程

 

最小二乘法的参数最优解

目标函数(loss/cost function)

普通最小二乘法线性回归案例

Date;Time;Global_active_power;Global_reactive_power;Voltage;Global_intensity;Sub_metering_1;Sub_metering_2;Sub_metering_3
16/12/2006;17:24:00;4.216;0.418;234.840;18.400;0.000;1.000;17.0
16/12/2006;17:25:00;5.360;0.436;233.630;23.000;0.000;1.000;16.0
16/12/2006;17:26:00;5.374;0.498;233.290;23.000;0.000;2.000;17.0
16/12/2006;17:27:00;5.388;0.502;233.740;23.000;0.000;1.000;17.0
16/12/2006;17:28:00;3.666;0.528;235.680;15.800;0.000;1.000;17.0
16/12/2006;17:29:00;3.520;0.522;235.020;15.000;0.000;2.000;17.0
16/12/2006;17:30:00;3.702;0.520;235.090;15.800;0.000;1.000;17.0
16/12/2006;17:31:00;3.700;0.520;235.220;15.800;0.000;1.000;17.0
16/12/2006;17:32:00;3.668;0.510;233.990;15.800;0.000;1.000;17.0
16/12/2006;17:33:00;3.662;0.510;233.860;15.800;0.000;2.000;16.0
16/12/2006;17:34:00;4.448;0.498;232.860;19.600;0.000;1.000;17.0
16/12/2006;17:35:00;5.412;0.470;232.780;23.200;0.000;1.000;17.0
16/12/2006;17:36:00;5.224;0.478;232.990;22.400;0.000;1.000;16.0
16/12/2006;17:37:00;5.268;0.398;232.910;22.600;0.000;2.000;17.0
16/12/2006;17:38:00;4.054;0.422;235.240;17.600;0.000;1.000;17.0
16/12/2006;17:39:00;3.384;0.282;237.140;14.200;0.000;0.000;17.0
16/12/2006;17:40:00;3.270;0.152;236.730;13.800;0.000;0.000;17.0
16/12/2006;17:41:00;3.430;0.156;237.060;14.400;0.000;0.000;17.0
16/12/2006;17:42:00;3.266;0.000;237.130;13.800;0.000;0.000;18.0
16/12/2006;17:43:00;3.728;0.000;235.840;16.400;0.000;0.000;17.0
16/12/2006;17:44:00;5.894;0.000;232.690;25.400;0.000;0.000;16.0
16/12/2006;17:45:00;7.706;0.000;230.980;33.200;0.000;0.000;17.0
16/12/2006;17:46:00;7.026;0.000;232.210;30.600;0.000;0.000;16.0
16/12/2006;17:47:00;5.174;0.000;234.190;22.000;0.000;0.000;17.0
16/12/2006;17:48:00;4.474;0.000;234.960;19.400;0.000;0.000;17.0
16/12/2006;17:49:00;3.248;0.000;236.660;13.600;0.000;0.000;17.0
16/12/2006;17:50:00;3.236;0.000;235.840;13.600;0.000;0.000;17.0
16/12/2006;17:51:00;3.228;0.000;235.600;13.600;0.000;0.000;17.0
16/12/2006;17:52:00;3.258;0.000;235.490;13.800;0.000;0.000;17.0
16/12/2006;17:53:00;3.178;0.000;235.280;13.400;0.000;0.000;17.0
16/12/2006;17:54:00;2.720;0.000;235.060;11.600;0.000;0.000;17.0
16/12/2006;17:55:00;3.758;0.076;234.170;16.400;0.000;0.000;17.0
16/12/2006;17:56:00;4.342;0.090;233.770;18.400;0.000;0.000;16.0
16/12/2006;17:57:00;4.512;0.000;233.620;19.200;0.000;0.000;17.0
16/12/2006;17:58:00;4.058;0.200;234.680;17.600;0.000;0.000;17.0
16/12/2006;17:59:00;2.472;0.058;236.940;10.400;0.000;0.000;17.0
16/12/2006;18:00:00;2.790;0.180;237.520;11.800;0.000;0.000;18.0
16/12/2006;18:01:00;2.624;0.144;238.200;11.000;0.000;0.000;17.0
16/12/2006;18:02:00;2.772;0.118;238.280;11.600;0.000;0.000;17.0
16/12/2006;18:03:00;3.740;0.108;236.930;16.400;0.000;16.000;18.0
16/12/2006;18:04:00;4.928;0.202;235.010;21.000;0.000;37.000;16.0
16/12/2006;18:05:00;6.052;0.192;232.930;26.200;0.000;37.000;17.0
16/12/2006;18:06:00;6.752;0.186;232.120;29.000;0.000;36.000;17.0
16/12/2006;18:07:00;6.474;0.144;231.850;27.800;0.000;37.000;16.0
16/12/2006;18:08:00;6.308;0.116;232.250;27.000;0.000;36.000;17.0
16/12/2006;18:09:00;4.464;0.136;234.660;19.000;0.000;37.000;16.0
16/12/2006;18:10:00;3.396;0.148;236.200;15.000;0.000;22.000;18.0
16/12/2006;18:11:00;3.090;0.152;237.070;13.800;0.000;12.000;17.0
16/12/2006;18:12:00;3.730;0.144;235.780;16.400;0.000;27.000;17.0
16/12/2006;18:13:00;2.308;0.160;237.430;9.600;0.000;1.000;17.0
16/12/2006;18:14:00;2.388;0.158;237.260;10.000;0.000;1.000;17.0
16/12/2006;18:15:00;4.598;0.100;234.250;21.400;0.000;20.000;17.0
16/12/2006;18:16:00;4.524;0.076;234.200;19.600;0.000;9.000;17.0
16/12/2006;18:17:00;4.202;0.082;234.310;17.800;0.000;1.000;17.0
16/12/2006;18:18:00;4.472;0.000;233.290;19.200;0.000;1.000;16.0
16/12/2006;18:19:00;2.852;0.000;235.610;12.000;0.000;1.000;17.0
16/12/2006;18:20:00;2.928;0.000;235.250;12.400;0.000;1.000;17.0
16/12/2006;18:21:00;2.940;0.000;236.040;12.400;0.000;2.000;17.0
16/12/2006;18:22:00;2.934;0.000;235.510;12.400;0.000;1.000;17.0
16/12/2006;18:23:00;2.926;0.000;235.680;12.400;0.000;1.000;17.0
16/12/2006;18:24:00;3.452;0.000;235.200;15.200;0.000;1.000;17.0
16/12/2006;18:25:00;4.870;0.000;233.740;20.800;0.000;1.000;17.0
16/12/2006;18:26:00;4.868;0.000;233.840;20.800;0.000;1.000;17.0
16/12/2006;18:27:00;4.866;0.000;233.790;20.800;0.000;1.000;17.0
16/12/2006;18:28:00;3.176;0.000;235.500;13.800;0.000;1.000;17.0
16/12/2006;18:29:00;2.920;0.000;235.840;12.400;0.000;1.000;17.0
16/12/2006;18:30:00;2.930;0.000;236.150;12.400;0.000;1.000;17.0
16/12/2006;18:31:00;2.912;0.050;235.810;12.400;0.000;1.000;17.0
16/12/2006;18:32:00;2.608;0.052;235.410;11.000;0.000;1.000;17.0
16/12/2006;18:33:00;2.714;0.162;234.820;11.600;0.000;0.000;17.0
16/12/2006;18:34:00;3.538;0.086;233.760;15.600;0.000;1.000;16.0
16/12/2006;18:35:00;6.072;0.000;232.480;26.400;0.000;27.000;17.0
16/12/2006;18:36:00;4.536;0.000;233.540;19.400;0.000;1.000;17.0
16/12/2006;18:37:00;4.408;0.000;232.320;18.800;0.000;1.000;16.0
16/12/2006;18:38:00;2.912;0.048;234.020;13.000;0.000;1.000;17.0
16/12/2006;18:39:00;2.326;0.054;234.760;9.800;0.000;1.000;17.0
16/12/2006;18:40:00;2.264;0.054;234.670;9.600;0.000;1.000;17.0
16/12/2006;18:41:00;2.270;0.054;235.270;9.600;0.000;1.000;17.0
16/12/2006;18:42:00;2.258;0.054;235.120;9.600;0.000;1.000;17.0
16/12/2006;18:43:00;2.188;0.068;235.800;9.200;0.000;1.000;17.0
16/12/2006;18:44:00;2.978;0.166;234.810;13.200;0.000;1.000;17.0
16/12/2006;18:45:00;4.200;0.174;234.380;17.800;0.000;1.000;17.0
16/12/2006;18:46:00;4.204;0.186;234.200;17.800;0.000;1.000;16.0
16/12/2006;18:47:00;4.218;0.178;233.980;18.000;0.000;1.000;17.0
16/12/2006;18:48:00;2.786;0.188;234.990;12.000;0.000;2.000;17.0
16/12/2006;18:49:00;2.540;0.088;234.670;10.800;0.000;4.000;17.0
16/12/2006;18:50:00;2.496;0.080;233.920;10.600;0.000;3.000;17.0
16/12/2006;18:51:00;2.336;0.070;233.510;10.000;0.000;1.000;16.0
16/12/2006;18:52:00;2.322;0.000;233.440;9.800;0.000;0.000;17.0
16/12/2006;18:53:00;2.448;0.000;233.640;10.600;0.000;1.000;17.0
16/12/2006;18:54:00;4.298;0.000;232.390;18.400;0.000;1.000;16.0
16/12/2006;18:55:00;4.230;0.090;232.250;18.200;0.000;1.000;17.0
16/12/2006;18:56:00;4.230;0.090;232.320;18.200;0.000;2.000;16.0
16/12/2006;18:57:00;3.924;0.084;232.790;17.000;0.000;1.000;17.0
16/12/2006;18:58:00;4.218;0.090;232.090;18.000;0.000;1.000;17.0
16/12/2006;18:59:00;4.224;0.090;231.960;18.200;0.000;1.000;16.0
16/12/2006;19:00:00;4.070;0.088;231.990;17.400;0.000;1.000;17.0
16/12/2006;19:01:00;3.612;0.090;232.360;15.600;0.000;2.000;16.0
16/12/2006;19:02:00;3.458;0.090;232.710;14.800;0.000;1.000;17.0
16/12/2006;19:03:00;3.434;0.090;232.010;14.800;0.000;1.000;16.0
16/12/2006;19:04:00;3.366;0.000;232.780;14.400;0.000;1.000;17.0
16/12/2006;19:05:00;3.404;0.082;233.080;14.600;0.000;2.000;16.0
16/12/2006;19:06:00;3.382;0.074;233.690;14.400;0.000;0.000;17.0
16/12/2006;19:07:00;3.616;0.106;234.440;15.400;0.000;4.000;17.0
16/12/2006;19:08:00;3.476;0.092;234.600;14.800;0.000;3.000;17.0
16/12/2006;19:09:00;3.394;0.064;235.120;14.400;0.000;0.000;17.0
16/12/2006;19:10:00;3.400;0.000;234.280;14.400;0.000;1.000;17.0
16/12/2006;19:11:00;3.414;0.000;234.190;14.400;0.000;1.000;17.0
16/12/2006;19:12:00;3.432;0.000;234.690;14.600;0.000;1.000;16.0
16/12/2006;19:13:00;3.418;0.000;234.260;14.600;0.000;2.000;17.0
16/12/2006;19:14:00;3.392;0.000;233.300;14.400;0.000;1.000;17.0
16/12/2006;19:15:00;3.394;0.000;233.430;14.400;0.000;1.000;17.0
16/12/2006;19:16:00;3.394;0.000;233.450;14.400;0.000;1.000;16.0
16/12/2006;19:17:00;3.406;0.000;233.910;14.400;0.000;2.000;17.0
16/12/2006;19:18:00;3.420;0.000;234.430;14.600;0.000;1.000;17.0
16/12/2006;19:19:00;3.418;0.000;234.340;14.400;0.000;1.000;17.0
16/12/2006;19:20:00;3.422;0.052;234.580;14.600;0.000;1.000;17.0
16/12/2006;19:21:00;3.332;0.000;234.020;14.200;0.000;1.000;16.0
16/12/2006;19:22:00;3.316;0.046;234.290;14.000;0.000;1.000;17.0
16/12/2006;19:23:00;3.334;0.000;234.360;14.200;0.000;2.000;17.0
16/12/2006;19:24:00;3.262;0.052;232.640;14.000;0.000;0.000;17.0
16/12/2006;19:25:00;3.476;0.064;232.960;14.800;0.000;3.000;16.0
16/12/2006;19:26:00;3.620;0.090;233.470;15.400;0.000;5.000;17.0
16/12/2006;19:27:00;3.610;0.084;234.980;15.200;0.000;5.000;17.0
16/12/2006;19:28:00;3.646;0.138;234.340;15.400;0.000;5.000;17.0
16/12/2006;19:29:00;3.614;0.192;233.920;15.400;0.000;5.000;17.0
16/12/2006;19:30:00;3.408;0.162;233.940;14.400;0.000;2.000;16.0
16/12/2006;19:31:00;3.388;0.134;233.720;14.400;0.000;1.000;17.0
16/12/2006;19:32:00;3.372;0.132;233.200;14.400;0.000;2.000;17.0
16/12/2006;19:33:00;3.370;0.132;233.100;14.400;0.000;1.000;16.0
16/12/2006;19:34:00;3.342;0.120;232.290;14.400;0.000;1.000;17.0
16/12/2006;19:35:00;3.284;0.000;232.450;14.000;0.000;2.000;16.0
16/12/2006;19:36:00;3.276;0.000;231.970;14.000;0.000;1.000;17.0
16/12/2006;19:37:00;3.382;0.130;232.420;14.400;0.000;2.000;17.0
16/12/2006;19:38:00;3.418;0.160;233.380;14.600;0.000;1.000;16.0
16/12/2006;19:39:00;3.396;0.158;232.690;14.600;0.000;1.000;17.0
16/12/2006;19:40:00;3.406;0.158;233.200;14.600;0.000;2.000;16.0
16/12/2006;19:41:00;3.406;0.158;233.260;14.600;0.000;1.000;17.0
16/12/2006;19:42:00;3.392;0.156;232.840;14.400;0.000;1.000;17.0
16/12/2006;19:43:00;3.378;0.154;232.390;14.400;0.000;2.000;16.0
16/12/2006;19:44:00;3.368;0.154;232.140;14.400;0.000;1.000;17.0
16/12/2006;19:45:00;3.350;0.152;231.570;14.400;0.000;1.000;16.0
16/12/2006;19:46:00;3.352;0.152;231.660;14.400;0.000;2.000;17.0
16/12/2006;19:47:00;3.364;0.154;232.220;14.400;0.000;1.000;16.0
16/12/2006;19:48:00;3.380;0.156;232.810;14.400;0.000;1.000;17.0
16/12/2006;19:49:00;3.404;0.160;233.770;14.400;0.000;1.000;17.0
16/12/2006;19:50:00;3.378;0.156;232.820;14.400;0.000;2.000;16.0
16/12/2006;19:51:00;3.388;0.158;233.220;14.400;0.000;1.000;17.0
16/12/2006;19:52:00;3.376;0.156;232.840;14.400;0.000;1.000;17.0
16/12/2006;19:53:00;3.254;0.000;232.750;13.800;0.000;0.000;16.0
16/12/2006;19:54:00;3.214;0.078;232.640;13.800;0.000;0.000;17.0
16/12/2006;19:55:00;3.202;0.078;232.200;13.800;0.000;0.000;16.0
16/12/2006;19:56:00;3.202;0.078;232.200;13.800;0.000;0.000;17.0
16/12/2006;19:57:00;3.212;0.078;232.580;13.800;0.000;0.000;16.0
16/12/2006;19:58:00;3.206;0.078;232.400;13.800;0.000;0.000;17.0
16/12/2006;19:59:00;3.214;0.078;232.660;13.800;0.000;0.000;17.0
16/12/2006;20:00:00;3.206;0.078;232.430;13.800;0.000;0.000;16.0
16/12/2006;20:01:00;3.222;0.078;233.020;13.800;0.000;0.000;17.0
16/12/2006;20:02:00;3.216;0.076;232.800;13.800;0.000;0.000;16.0
16/12/2006;20:03:00;3.216;0.078;232.790;13.800;0.000;0.000;17.0
16/12/2006;20:04:00;3.194;0.078;231.980;13.600;0.000;0.000;16.0
16/12/2006;20:05:00;3.278;0.000;232.880;14.000;0.000;0.000;17.0
16/12/2006;20:06:00;3.282;0.000;232.520;14.000;0.000;0.000;17.0
16/12/2006;20:07:00;3.292;0.000;232.870;14.000;0.000;0.000;16.0
16/12/2006;20:08:00;3.316;0.000;233.780;14.000;0.000;0.000;17.0
16/12/2006;20:09:00;3.370;0.080;233.690;14.400;0.000;0.000;17.0
16/12/2006;20:10:00;3.420;0.144;234.040;14.600;0.000;0.000;17.0
16/12/2006;20:11:00;3.410;0.142;233.930;14.600;0.000;0.000;16.0
16/12/2006;20:12:00;3.400;0.140;233.620;14.400;0.000;0.000;17.0
16/12/2006;20:13:00;3.400;0.140;233.640;14.400;0.000;0.000;17.0
16/12/2006;20:14:00;3.388;0.140;233.260;14.400;0.000;0.000;16.0
16/12/2006;20:15:00;3.370;0.086;233.900;14.400;0.000;0.000;17.0
16/12/2006;20:16:00;3.332;0.000;233.680;14.200;0.000;0.000;17.0
16/12/2006;20:17:00;3.310;0.000;232.830;14.200;0.000;0.000;17.0
16/12/2006;20:18:00;3.312;0.000;233.000;14.200;0.000;0.000;16.0
16/12/2006;20:19:00;3.314;0.000;233.080;14.200;0.000;0.000;17.0
16/12/2006;20:20:00;3.302;0.000;232.690;14.200;0.000;0.000;16.0
16/12/2006;20:21:00;3.302;0.000;232.910;14.000;0.000;0.000;17.0
16/12/2006;20:22:00;3.300;0.000;232.870;14.000;0.000;0.000;17.0
16/12/2006;20:23:00;3.290;0.000;232.420;14.000;0.000;0.000;16.0
16/12/2006;20:24:00;3.286;0.000;232.310;14.000;0.000;0.000;17.0
16/12/2006;20:25:00;3.230;0.052;232.790;13.800;0.000;0.000;16.0
16/12/2006;20:26:00;3.218;0.078;233.050;13.800;0.000;0.000;17.0
16/12/2006;20:27:00;3.258;0.076;234.570;13.800;0.000;0.000;17.0
16/12/2006;20:28:00;3.264;0.076;234.790;13.800;0.000;0.000;17.0
16/12/2006;20:29:00;3.260;0.076;234.600;13.800;0.000;0.000;17.0
16/12/2006;20:30:00;3.262;0.076;234.540;13.800;0.000;0.000;17.0
16/12/2006;20:31:00;3.264;0.076;234.600;13.800;0.000;0.000;16.0
16/12/2006;20:32:00;3.238;0.076;233.770;13.800;0.000;0.000;17.0
16/12/2006;20:33:00;3.228;0.078;233.500;13.800;0.000;0.000;17.0
16/12/2006;20:34:00;3.222;0.078;233.270;13.800;0.000;0.000;17.0
16/12/2006;20:35:00;3.226;0.078;233.370;13.800;0.000;0.000;16.0
16/12/2006;20:36:00;3.206;0.078;232.640;13.600;0.000;0.000;17.0
16/12/2006;20:37:00;3.188;0.080;231.950;13.600;0.000;0.000;16.0
16/12/2006;20:38:00;3.196;0.080;232.210;13.600;0.000;0.000;17.0
16/12/2006;20:39:00;3.204;0.078;232.500;13.600;0.000;0.000;16.0
16/12/2006;20:40:00;3.232;0.078;233.630;13.800;0.000;0.000;17.0
16/12/2006;20:41:00;3.254;0.074;234.430;13.800;0.000;0.000;17.0
16/12/2006;20:42:00;3.376;0.050;234.630;14.400;0.000;0.000;17.0
16/12/2006;20:43:00;3.372;0.048;235.050;14.200;0.000;0.000;17.0

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False

# 1. 加载数据
# path = '../datas/household_power_consumption_200.txt'
path = '../datas/household_power_consumption_1000.txt'
df = pd.read_csv(path, sep=';')
# print(df.head(2))
# df.info()

# 2. 获取功率值作为作为特征属性X，电流作为目标属性Y
X = df.iloc[:, 2:4]
# print(X.head(2))
Y = df.iloc[:, 5]

# 3. 获取训练数据和测试数据
n = int(X.shape[0] * 0.8)
train_x = np.array(X[:n])
test_x = np.array(X[n:])
train_y = np.array(Y[:n])
test_y = np.array(Y[n:])
print("总样本数目:{}, 训练数据样本数目:{}, 测试数据样本数目:{}".format(X.shape, train_x.shape, test_x.shape))

# 4. 训练模型
# a. 训练数据转换为矩阵形式
x = np.mat(train_x)
y = np.mat(train_y).reshape(-1, 1)
# b. 训练模型参数θ值
theta = (x.T * x).I * x.T * y
print(theta.shape)
print("求解出来的theta值:{}".format(theta))

# 5. 模型效果评估
y_hat = np.mat(test_x) * theta
# 画图看一下效果
t = np.arange(len(test_x))
plt.figure(facecolor='w')
plt.plot(t, y_hat, 'g-', linewidth=2, label=u'预测值')
plt.plot(t, test_y, 'r-', linewidth=2, label=u'真实值')
plt.legend(loc='lower right')
plt.title('线性回归')
plt.show()

# 6. 模型的存储, 这里可以直接将θ值保存到数据库中，然后再需要的程序中，再讲θ加载到程序中
# TODO: 模拟一下加载的操作
theta1 = theta[0]
theta2 = theta[1]

# 产生预测值
global_active_power = 4.216
global_reactive_power = 0.418
print("当前的输入的特征属性值为:{}---------{}".format(global_active_power, global_reactive_power))
print("预测值为:{}".format(global_active_power * theta1 + global_reactive_power * theta2))

使用sklearn实现模型的训练

https://scikit-learn.org/0.19/modules/generated/sklearn.linear_model.LinearRegression.html#sklearn.linear_model.LinearRegression

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.externals import joblib

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False

# 一、加载数据
# path = '../datas/household_power_consumption_201.txt'
path = '../datas/household_power_consumption_1000.txt'
df = pd.read_csv(path, sep=';')
# df.info()


# 二、数据的清洗
"""
可选操作
比如异常数据的处理，缺省数据的填充....
"""
# inplace: 该参数的作用是给定，是否是在原始的DataFrame上做替换，设置为true表示在原始DataFrame上做替换，直接修改原始的DataFrame的数据格式；默认是False
# df = df.replace('?', np.nan)
df.replace('?', np.nan, inplace=True)
# 删除所有具有非数字的样本
# how: 指定如何删除数据，如果设置为any，默认值, 那么只要行或者列中存在nan值，那么就进行删除操作；如果设置为all，那边要求行或者列中的所有值均为nan的时候，才可以进行删除操作
df.dropna(axis=0, how='any', inplace=True)


# 三、基于业务提取最原始的特征属性X和目标属性Y
X = df[['Global_active_power', 'Global_reactive_power']]
# X.info()
Y = df['Global_intensity']
# print("原始X形状:{}， 原始Y形状:{}".format(X.shape, Y.shape))


# 四、数据的划分(将数据划分为训练集和测试集)
# train_size：给定划分中的训练集占比，一般情况下，我们训练集和测试集的占比一般为8:2，如果样本比较多的时候，可以考虑7:3，如果样本比较少的时候，可以考虑9:1
# test_size: 给定划分中的测试集的占比，其中train_size和test_size只能有且给定其中的一个。
x_train, x_test, y_train, y_test = train_test_split(X, Y, train_size=0.8, random_state=28)
# print("训练集样本形状:{}, 测试集样本形状:{}".format(x_train.shape, x_test.shape))


# 五、特征工程
"""
可选操作
主要就是一些：哑编码、TF-IDF、连续数据的离散化、标准化、归一化、特征选择、降维....
NOTE: 所有的特征工程其实和算法模型训练一样，特征工程的意思就是对特征属性做一个数据的转换操作，所以需要获取一个对应数据转换函数---> 需要从训练数据中训练出来
"""
"""
NOTE:
sklearn中，所有特征工程、算法模型的API基本类似，主要是以下几个API：
fit: 模型/函数训练，基于给定的训练集数据，训练模型对象
transform: 使用训练好的模型对给定的数据集做一个对应的数据转换(使用训练好的函数转换数据)，一般出现在特征工程的相关的类中间
predict：使用构建好的算法模型对样本数据做一个预测的操作，其实内部也就是使用训练好的函数转换特征属性到目标属性，一般出现在算法模型的相关的类中间
fit_transform: 也就是fit和transform的结合，底层先使用给定的训练数据训练模型，找出对应的转换函数或者模型对象，然后使用找到的函数或者模型对训练数据做一个转换操作
"""
# 做一下数据的标准化操作
"""
实现方式一：
# a. 创建标准化的对象
ss = StandardScaler()
# b. 训练标准化的转换函数
ss.fit(x_train, y_train)
# c. 使用训练好的标注化模型对数据做一个转换操作
x_train = ss.transform(x_train)
x_test = ss.transform(x_test)

"""
# 实现方式二：
# a. 创建标准化的对象
ss = StandardScaler()
# b. 训练标准化的转换函数并使用训练好的函数对训练数据做一个转换操作
x_train = ss.fit_transform(x_train, y_train)
# c. 使用训练好的标注化模型对数据做一个转换操作
x_test = ss.transform(x_test)


# 六、算法模型的选择/算法模型对象的构建
# fit_intercept：是否训练模型的截距项，默认为True，表示训练；如果设置为False，表示不训练。
algo = LinearRegression(fit_intercept=True)


# 七、算法模型的训练
algo.fit(x_train, y_train)

# 7.1 查看训练好的模型参数
print("线性回归的各个特征属性对应的权重参数θ:{}".format(algo.coef_))
print("线性回归的截距项的值:{}".format(algo.intercept_))


# 八、模型效果评估
y_hat = algo.predict(x_test)
print("在训练集上的模型效果(回归算法中为R2):{}".format(algo.score(x_train, y_train)))
print("在测试集上的模型效果(回归算法中为R2):{}".format(algo.score(x_test, y_test)))
print("在测试集上的MSE的值:{}".format(mean_squared_error(y_true=y_test, y_pred=y_hat)))
# 画图看一下效果
t = np.arange(len(x_test))
plt.figure(facecolor='w')
plt.plot(t, y_hat, 'g-', linewidth=2, label=u'预测值')
plt.plot(t, y_test, 'r-', linewidth=2, label=u'真实值')
plt.legend(loc='lower right')
plt.title('线性回归')
plt.show()


# 九、模型保存
"""
两种保存方式：
1. 直接将模型输出为二进制的磁盘文件
2. 直接将预测结果y_hat输出到数据库
"""
# 查看python源码的API说明的方式：键盘按着ctrl键，鼠标放到需要查看的API上面，然后左键点击，进入源码即可查看
joblib.dump(ss, './model/ss.pkl')
joblib.dump(algo, './model/lr.pkl')

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import time

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False


def date_format(dt):
    t = time.strptime(' '.join(dt), '%d/%m/%Y %H:%M:%S')
    return (t.tm_year, t.tm_mon, t.tm_mday, t.tm_hour, t.tm_min, t.tm_sec)


# 加载数据
# path = '../datas/household_power_consumption_200.txt'  ## 200行数据
path = '../datas/household_power_consumption_1000.txt'  ## 1000行数据
df = pd.read_csv(path, sep=';', low_memory=False)

# 日期、时间、有功功率、无功功率、电压、电流、厨房用电功率、洗衣服用电功率、热水器用电功率
names2 = df.columns
names = ['Date', 'Time', 'Global_active_power', 'Global_reactive_power', 'Voltage', 'Global_intensity',
         'Sub_metering_1', 'Sub_metering_2', 'Sub_metering_3']

# 异常数据处理(异常数据过滤)
new_df = df.replace('?', np.nan)
datas = new_df.dropna(axis=0, how='any')  # 只要有数据为空，就进行删除操作

# 获取x和y变量, 并将时间转换为数值型连续变量
X = datas[names[0:2]]
X = X.apply(lambda x: pd.Series(date_format(x)), axis=1)
Y = datas[names[4]]
X = X.astype(np.float)
Y = Y.astype(np.float)

# 对数据集进行测试集合训练集划分
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0)

# 1. 构建一个管道流对象，定义数据处理的顺序
model = Pipeline(steps=[
    ('Poly', PolynomialFeatures()),  # 给定第一步操作，名称为Poly
    ('Linear', LinearRegression(fit_intercept=False))  # 给定第二步操作，名称为Linear
])

# 1.2 Pipeline对象设置参数
# Poly__degree: Poly是定义Pipeline对象的时候给定的步骤名称，degree是对应步骤对象中的属性名称, 中间是两个连续的下划线
model.set_params(Poly__degree=4)
model.set_params(Linear__normalize=True)

# 2. 模型训练（先调用第一步进行数据处理，然后再调用第二步做模型训练）
# 假设我们的步骤是n步，那么前n-1步做的操作是: fit + transform， 最后一步做的操作是fit
model.fit(X_train, Y_train)
"""
model.fit等价于linear.fit(poly.fit_transform(x_train,y_train),y_train)
"""

print("多项式模型:{}".format(model.get_params()['Poly']))
print("线性回归模型:{}".format(model.get_params()['Linear']))
# 3. 预测值产生(先调用第一步的transform对数据转换，再调用第二步的predict对数据预测)
# 假设我们的步骤是n步，那么前n-1步做的操作是: transform， 最后一步做的操作是fit
y_hat = model.predict(X_test)

# 7.1 查看训练好的模型参数
linear_algo = model.get_params()['Linear']
print("线性回归的各个特征属性对应的权重参数θ:{}".format(linear_algo.coef_))
print("线性回归的截距项的值:{}".format(linear_algo.intercept_))

# 八、模型效果评估
print("在训练集上的模型效果(回归算法中为R2):{}".format(model.score(X_train, Y_train)))
print("在测试集上的模型效果(回归算法中为R2):{}".format(model.score(X_test, Y_test)))
print("在测试集上的MSE的值:{}".format(mean_squared_error(y_true=Y_test, y_pred=y_hat)))
print("在测试集上的RMSE的值:{}".format(np.sqrt(mean_squared_error(y_true=Y_test, y_pred=y_hat))))

# 画图看一下效果
t = np.arange(len(X_test))
plt.figure(facecolor='w')
plt.plot(t, y_hat, 'g-', linewidth=2, label=u'预测值')
plt.plot(t, Y_test, 'r-', linewidth=2, label=u'真实值')
plt.legend(loc='lower right')
plt.title('线性回归')
plt.show()

加载使用sklearn构建好的模型对数据做一个预测

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.externals import joblib

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False

# 1. 加载模型
ss = joblib.load('./model/ss.pkl')
algo = joblib.load('./model/lr.pkl')

# 2. 使用加载模型对数据做一个预测操作
# a. 加载数据
path = '../datas/household_power_consumption_201.txt'
df = pd.read_csv(path, sep=';')
# b. 数据的清洗
df.replace('?', np.nan, inplace=True)
df.dropna(axis=0, how='any', inplace=True)
# c. 得到需要预测的x
x = df[['Global_active_power', 'Global_reactive_power']]
y = df['Global_intensity']

# d. 使用模型对x做一个预测的操作
y_predict = algo.predict(ss.transform(x))
print("预测结果:")
print(y_predict)

# e. 看一下预测效果
print("预测的效果:{}".format(algo.score(ss.transform(x), y)))
# 画图看一下效果
t = np.arange(len(x))
plt.figure(facecolor='w')
plt.plot(t, y_predict, 'g-', linewidth=2, label=u'预测值')
plt.plot(t, y, 'r-', linewidth=2, label=u'真实值')
plt.legend(loc='lower right')
plt.title('线性回归')
plt.show()

模型优化

算法模型构建的过程中，需要注意一下两个方面：

• 算法模型的误差不能过大(预测值和实际值之间的误差不能过大)
• 算法模型不能太复杂

如果模型的误差较大，也就是预测值和实际值之间的差值比较大的情况下：

1. 如果是训练数据集存在该问题： --> 欠拟合
表示此时模型训练的不够好，模型没有学习到数据的特征信息，也就是没有找到特征属性和目标属性之间的映射关系。
问题描述：当前称模型存在欠拟合的问题
导致原因：

• 1. 算法的学习能力比较弱
• 2. 数据量样本少
• 3. 有用的特征属性数目太少

解决方案：

• 1. 换一种学习能力比较强的算法，eg：svm、集成算法....
• 2. 增加样本数据(基于现有数据产生模拟的虚假数据（多项式扩展）、收集更多的数据)
• 3. 增加特征属性的数目(结合业务产生更多的有效特征、将现有的特征属性做一个融合从而产生更多的特征属性)

2. 如果仅是在测试数据集上存在该问题： --> 过拟合

如果模型在训练集上的效果不错，但是在测试集上的效果非常差，或者说训练集上的效果和测试集上的效果不匹配：
这种情况下，我们认为模型存在过拟合；也就是说模型的学习能力太强了，把训练数据的特征学习的太细了，导致学习了某些只存在于训练集中，而在测试集或者真实数据中不存在的特征信息
直白来讲：模型的学习能力太强，就可能存在过拟合的情况
产生原因：

• 1. 数据样本少
• 2. 模型的学习能力太强(模型比较复杂)
• 3. 做了太多的特征的增维操作

解决方案：

• 1. 增加数据样本(基于现有数据产生模拟的虚假数据、收集更多的数据)
• 2. 换一个算法模型或者在模型训练过程中，加入正则化项系数，限制模型过拟合，正则化主要分为：L1和L2
• 3. 不要做太多的、太深的维度的增维操作

多项式扩展
属于增维的一种方式，通过这种方式可以将数据映射到高维度空间变成线性可分的数据
功能：将低维空间上的数据通过多项式的组合，映射到高维空间中
效果：可以将低维空间的非线性的数据转换到高维空间中变成线性数据
方式(sklearn中的多项式扩展的API的效果)：
原始数据: (2,3)

• 如果做一个最高次项为2的多项式变换，最终结果为: (1,2,3,4,6,9)
• 如果做一个最高次项为3的多项式转换，最终结果为：(1,2,3,4,6,9,8,12,18,27)
• 如果做一个最高次项为4的多项式转换，最终结果为：(1,2,3,4,6,9,8,12,18,27,16,24,36,54,81)

def poly_expand(data,depth):
    result=np.array([1])
    element_index = [0]
    data_len = len(data)
    new_result = np.array([])
    for n in range(depth):
        if n == 0:
            new_result = data
            result = np.hstack((result,new_result))
        elif n == 1:
            new_result = []
            for i in range(data_len):
                for j in range(i,data_len):
                    new_result.append(data[i]*data[j])
                element_index.append(element_index[-1]+data_len-i)
            new_result = np.array(new_result)
            result = np.hstack((result,new_result))
        else:
            new_result_stack = np.array([])
            new_index = [0]
            old_result = new_result
            for i in range(data_len):
                index = element_index[i]
                new_result = data[i]*old_result[index:]
                new_result_stack = np.hstack((new_result_stack,new_result))
                new_index.append(new_index[-1] + element_index[-1]-index)
            new_result = new_result_stack
            result = np.hstack((result, new_result))
            element_index = new_index
    return result

def array_poly_expand(datas,depth):
    shape = datas.shape
    datas = datas.reshape(-1,shape[-1])
    new_datas = []
    for data in datas:
        res = poly_expand(data,depth)
        new_datas.append(res)
    return np.array(new_datas)

多项式回归算法

https://scikit-learn.org/0.19/modules/generated/sklearn.preprocessing.PolynomialFeatures.html#sklearn.preprocessing.PolynomialFeatures

多项式扩展+线性回归 ---> 多项式线性回归：其实相当于首先对数据做一个维度的扩展，然后对扩展之后的特征属性矩阵X做一个线性回归的算法模型的训练。

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import time

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False


def date_format(dt):
    t = time.strptime(' '.join(dt), '%d/%m/%Y %H:%M:%S')
    return (t.tm_year, t.tm_mon, t.tm_mday, t.tm_hour, t.tm_min, t.tm_sec)


# 加载数据
# path = '../datas/household_power_consumption_200.txt'  ## 200行数据
path = '../datas/household_power_consumption_1000.txt'  ## 1000行数据
df = pd.read_csv(path, sep=';', low_memory=False)

# 日期、时间、有功功率、无功功率、电压、电流、厨房用电功率、洗衣服用电功率、热水器用电功率
names2 = df.columns

names = ['Date', 'Time', 'Global_active_power', 'Global_reactive_power', 'Voltage', 'Global_intensity',
         'Sub_metering_1', 'Sub_metering_2', 'Sub_metering_3']

# 异常数据处理(异常数据过滤)
new_df = df.replace('?', np.nan)
datas = new_df.dropna(axis=0, how='any')  # 只要有数据为空，就进行删除操作

# 获取x和y变量, 并将时间转换为数值型连续变量
X = datas[names[0:2]]
X = X.apply(lambda x: pd.Series(date_format(x)), axis=1)
Y = datas[names[4]]
X = X.astype(np.float)
Y = Y.astype(np.float)

# 对数据集进行测试集合训练集划分
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0)

# 1. 做一个维度多项式扩展的操作
# degree：给定模型做几阶的多项式扩展，也就是转换后的最高次项是多少
poly = PolynomialFeatures(degree=3)
# fit_transform：首先使用给定的数据集进行模型训练，找出模型的转换函数，然后使用找出的转换函数对给定的X数据做一个转换操作
X_train = poly.fit_transform(X_train, Y_train)
X_test = poly.transform(X_test)

# 2. 做一个线性回归
#  fit_intercept：是否训练模型的截距项，默认为True，表示训练；如果设置为False，表示不训练。
algo = LinearRegression(fit_intercept=True)

# 七、算法模型的训练
algo.fit(X_train, Y_train)

# 7.1 查看训练好的模型参数
print("线性回归的各个特征属性对应的权重参数θ:{}".format(algo.coef_))
print("线性回归的截距项的值:{}".format(algo.intercept_))

# 八、模型效果评估
y_hat = algo.predict(X_test)
print("在训练集上的模型效果(回归算法中为R2):{}".format(algo.score(X_train, Y_train)))
print("在测试集上的模型效果(回归算法中为R2):{}".format(algo.score(X_test, Y_test)))
print("在测试集上的MSE的值:{}".format(mean_squared_error(y_true=Y_test, y_pred=y_hat)))
# 画图看一下效果
t = np.arange(len(X_test))
plt.figure(facecolor='w')
plt.plot(t, y_hat, 'g-', linewidth=2, label=u'预测值')
plt.plot(t, Y_test, 'r-', linewidth=2, label=u'真实值')
plt.legend(loc='lower right')
plt.title('线性回归')
plt.show()

正则化

使用L2正则的线性回归模型就称为Ridge回归(岭回归)。使用L1正则的线性回归模型就称为LASSO回归(Least Absolute Shrinkage and Selection Operator)

• L2-norm中，由于对于各个维度的参数缩放是在一个圆内缩放的，不可能导致有维度参数变为0的情况，那么也就不会产生稀疏解；实际应用中，数据的维度中是存在噪音和冗余的，稀疏的解可以找到有用的维度并且减少冗余，提高回归预测的准确性和鲁棒性（减少了overfitting）(L1-norm可以达到最终解的稀疏 性的要求)。
• Ridge模型具有较高的准确性、鲁棒性以及稳定性(冗余特征已经被删除了)。
• LASSO模型具有较高的求解速度。
• 如果既要考虑稳定性也考虑求解的速度，就使用Elasitc Net。

稀疏就是训练出来得到的模型参数中有很多参数值都是0，这个我们就叫做稀疏解；稀疏解的作用主要是用于特征选择的，因为参数为0所对应的特征相当于没有决策能力，不会影响y的取值，既然这样，我们可以将这些为0所对应的特征删除。

Elasitc Net

同时使用L1正则和L2正则的线性回归模型就称为Elasitc Net算法(弹性网络算法)

from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge,ElasticNet,Lasso
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import time

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False


def date_format(dt):
    t = time.strptime(' '.join(dt), '%d/%m/%Y %H:%M:%S')
    return (t.tm_year, t.tm_mon, t.tm_mday, t.tm_hour, t.tm_min, t.tm_sec)


# 加载数据
path = '../datas/household_power_consumption_200.txt'  ## 200行数据
# path = '../datas/household_power_consumption_1000.txt'  ## 1000行数据
df = pd.read_csv(path, sep=';', low_memory=False)

# 日期、时间、有功功率、无功功率、电压、电流、厨房用电功率、洗衣服用电功率、热水器用电功率
names2 = df.columns
names = ['Date', 'Time', 'Global_active_power', 'Global_reactive_power', 'Voltage', 'Global_intensity',
         'Sub_metering_1', 'Sub_metering_2', 'Sub_metering_3']

# 异常数据处理(异常数据过滤)
new_df = df.replace('?', np.nan)
datas = new_df.dropna(axis=0, how='any')  # 只要有数据为空，就进行删除操作

# 获取x和y变量, 并将时间转换为数值型连续变量
X = datas[names[0:2]]
X = X.apply(lambda x: pd.Series(date_format(x)), axis=1)
Y = datas[names[4]]
X = X.astype(np.float)
Y = Y.astype(np.float)

# 对数据集进行测试集合训练集划分
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0)

# 1. 构建一个管道流对象，定义数据处理的顺序
"""
Ridge参数说明: 
  alpha=1.0: ppt上的λ，正则化系数, 
  fit_intercept=True: 模型是否训练截距项，默认为训练(True), 
  normalize=False：在做模型训练之前，是否做归一化操作，一般不改动,
  max_iter=None：模型求解的迭代最大次数，默认表示不限制, 
  tol=1e-3: 模型收敛条件，当损失函数的变化值小于该值的时候，介绍迭代更新, 
  solver="auto"：给定求解方式,
"""
model = Pipeline(steps=[
    ('Poly', PolynomialFeatures()),  # 给定第一步操作，名称为Poly
    ('Linear', Ridge(alpha=0.1, fit_intercept=False))  # 给定第二步操作，名称为Linear
])

# 1.2 Pipeline对象设置参数
# Poly__degree: Poly是定义Pipeline对象的时候给定的步骤名称，degree是对应步骤对象中的属性名称, 中间是两个连续的下划线
model.set_params(Poly__degree=4)
model.set_params(Linear__normalize=True)

# 2. 模型训练（先调用第一步进行数据处理，然后再调用第二步做模型训练）
# 假设我们的步骤是n步，那么前n-1步做的操作是: fit + transform， 最后一步做的操作是fit
model.fit(X_train, Y_train)
"""
model.fit等价于linear.fit(poly.fit_transform(x_train,y_train),y_train)
"""

print("多项式模型:{}".format(model.get_params()['Poly']))
print("线性回归模型:{}".format(model.get_params()['Linear']))
# 3. 预测值产生(先调用第一步的transform对数据转换，再调用第二步的predict对数据预测)
# 假设我们的步骤是n步，那么前n-1步做的操作是: transform， 最后一步做的操作是fit
y_hat = model.predict(X_test)

# 7.1 查看训练好的模型参数
linear_algo = model.get_params()['Linear']
print("线性回归的各个特征属性对应的权重参数θ:{}".format(linear_algo.coef_))
print("线性回归的截距项的值:{}".format(linear_algo.intercept_))

# 八、模型效果评估
print("在训练集上的模型效果(回归算法中为R2):{}".format(model.score(X_train, Y_train)))
print("在测试集上的模型效果(回归算法中为R2):{}".format(model.score(X_test, Y_test)))
print("在测试集上的MSE的值:{}".format(mean_squared_error(y_true=Y_test, y_pred=y_hat)))
print("在测试集上的RMSE的值:{}".format(np.sqrt(mean_squared_error(y_true=Y_test, y_pred=y_hat))))

# 画图看一下效果
t = np.arange(len(X_test))
plt.figure(facecolor='w')
plt.plot(t, y_hat, 'g-', linewidth=2, label=u'预测值')
plt.plot(t, Y_test, 'r-', linewidth=2, label=u'真实值')
plt.legend(loc='lower right')
plt.title('线性回归')
plt.show()

模型效果判断（评估）

MSE：误差平方和，越趋近于0表示模型越拟合训练数据。

RMSE：MSE的平方根，作用同MSE。

R2：取值范围(负无穷,1]，值越大表示模型越拟合训练数据；最优解是1；当模型 预测为随机值的时候，有可能为负；若预测值恒为样本期望，R2为0。

TSS：总平方和TSS(Total Sum of Squares)，表示样本之间的差异情况，是伪方 差的m倍。

RSS：残差平方和RSS（Residual Sum of Squares），表示预测值和样本值之 间的差异情况，是MSE的m倍。

机器学习调参

模型参数：需要在训练数据集上通过某种给定的方式找出的模型参数，也就是说这个模型参数的求解就是我们经常所说的模型学习，eg：线性回归中的θ....
超参：在模型训练中需要使用到的参数值，但是该参数值需要开发人员给定的。eg：Ridge API中的alpha....
给定超参数的方式：

• 可以根据算法的特性、业务背景以及经验，来给定一个比较适合的值
• 通过sklearn提供的交叉验证的方式来选择最优的参数
• 通过网格交叉验证的方式选择最优参数

在实际工作中，对于各种算法模型(线性回归)来讲，我们需要获取θ、λ、p的值； θ的求解其实就是算法模型的求解，一般不需要开发人员参与(算法已经实现)，主要需要求解的是λ和p的值，这个过程就叫做调参(超参)

交叉验证：将训练数据分为多份，其中一份进行数据验证并获取最优的超参：λ 和p；比如：十折交叉验证、五折交叉验证(scikit-learn中默认)等

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import warnings
import sklearn
from sklearn.linear_model import LinearRegression, LassoCV, RidgeCV, ElasticNetCV
from sklearn.preprocessing import PolynomialFeatures#数据预处理，标准化
from sklearn.pipeline import Pipeline
from sklearn.linear_model.coordinate_descent import ConvergenceWarning

## 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif']=['SimHei']
mpl.rcParams['axes.unicode_minus']=False
## 拦截异常
warnings.filterwarnings(action = 'ignore', category=ConvergenceWarning)

## 创建模拟数据
np.random.seed(100)
np.set_printoptions(linewidth=1000, suppress=True)#显示方式设置，每行的字符数用于插入换行符，是否使用科学计数法
N = 10
x = np.linspace(0, 6, N) + np.random.randn(N)
y = 1.8*x**3 + x**2 - 14*x - 7 + np.random.randn(N)
## 将其设置为矩阵
x.shape = -1,1
y.shape = -1,1

## RidgeCV和Ridge的区别是：前者可以进行交叉验证
models = [
    Pipeline([
            ('Poly', PolynomialFeatures(include_bias=True)),
            ('Linear', LinearRegression(fit_intercept=False))
        ]),
    Pipeline([
            ('Poly', PolynomialFeatures(include_bias=True)),
            # alpha给定的是Ridge算法中，L2正则项的权重值，也就是ppt中的兰姆达
            # alphas是给定CV交叉验证过程中，Ridge算法的alpha参数值的取值的范围
            ('Linear', RidgeCV(alphas=np.logspace(-3,2,50), fit_intercept=False))
        ]),
    Pipeline([
            ('Poly', PolynomialFeatures(include_bias=True)),
            ('Linear', LassoCV(alphas=np.logspace(0,1,10), fit_intercept=False))
        ]),
    Pipeline([
            ('Poly', PolynomialFeatures(include_bias=True)),
            # la_ratio：给定EN算法中L1正则项在整个惩罚项中的比例，这里给定的是一个列表；
            # 表示的是在CV交叉验证的过程中，EN算法L1正则项的权重比例的可选值的范围
            ('Linear', ElasticNetCV(alphas=np.logspace(0,1,10), l1_ratio=[.1, .5, .7, .9, .95, 1], fit_intercept=False))
        ])
]

## 线性模型过拟合图形识别
fig,axes = plt.subplots(nrows=4,ncols=3,figsize=(30,40))
degree = np.arange(1,N,4) # 阶
dm = degree.size
colors = [] # 颜色
for c in np.linspace(16711680, 255, dm):
    colors.append('#%06x' % int(c))
names = ['线性回归','L2正则','L1正则','弹性网络']
for j in range(len(models)):
    model = models[j]
    for i,d in enumerate(degree):
        axes[j,i].plot(x, y, 'ro',ms=10, zorder=N)
        # 设置阶数
        model.set_params(Poly__degree=d)
        # 模型训练
        model.fit(x, y.ravel())

        lin = model.get_params()['Linear']
        output = u'%d阶，系数为：' % (d)
        print (output, lin.coef_.ravel())

        x_hat = np.linspace(x.min(), x.max(), num=100) ## 产生模拟数据
        x_hat.shape = -1,1
        y_hat = model.predict(x_hat)
        s = model.score(x, y)

        z = N - 1 if (d == 4) else 0
        axes[j,i].plot(x_hat, y_hat, color=colors[i], lw=2, alpha=0.75,zorder=z)
        axes[j,i].set(title='%s:%d阶, 正确率=%.3f' % (names[j],d,s),xlabel='X',ylabel='Y')
        axes[j,i].grid(True)
#         axes[j,i].legend(loc = 'upper left')
#         plt.tight_layout(1, rect=(0,0,1,0.95))
plt.show()

## 线性回归、Lasso回归、Ridge回归、ElasticNet比较
plt.figure(facecolor='w')
degree = np.arange(1,N, 2) # 阶， 多项式扩展允许给定的阶数
dm = degree.size
colors = [] # 颜色
for c in np.linspace(16711680, 255, dm):
    colors.append('#%06x' % int(c))
titles = [u'线性回归', u'Ridge回归', u'Lasso回归', u'ElasticNet']

for t in range(4):
    model = models[t]#选择了模型--具体的pipeline(线性、Lasso、Ridge、EN)
    plt.subplot(2,2,t+1) # 选择具体的子图
    plt.plot(x, y, 'ro', ms=10, zorder=N - 1) # 在子图中画原始数据点； zorder：图像显示在第几层

    # 遍历不同的多项式的阶，看不同阶的情况下，模型的效果
    for i,d in enumerate(degree):
        # 设置阶数(多项式)
        model.set_params(Poly__degree=d)
        # 模型训练
        model.fit(x, y.ravel())

        # 获取得到具体的算法模型
        # model.get_params()方法返回的其实是一个dict对象，后面的Linear其实是dict对应的key
        # 也是我们在定义Pipeline的时候给定的一个名称值
        lin = model.get_params()['Linear']
        # 打印数据
        output = u'%s:%d阶，系数为：' % (titles[t],d)
        # 判断lin对象中是否有对应的属性
        if hasattr(lin, 'alpha_'): # 判断lin这个模型中是否有alpha_这个属性
            idx = output.find(u'系数')
            output = output[:idx] + (u'alpha=%.6f, ' % lin.alpha_) + output[idx:]
        if hasattr(lin, 'l1_ratio_'): # 判断lin这个模型中是否有l1_ratio_这个属性
            idx = output.find(u'系数')
            output = output[:idx] + (u'l1_ratio=%.6f, ' % lin.l1_ratio_) + output[idx:]
        # line.coef_：获取线性模型的参数列表，也就是我们ppt中的theta值，ravel()将结果转换为1维数据
        print (output, lin.coef_.ravel())

        # 产生模拟数据
        x_hat = np.linspace(x.min(), x.max(), num=100) ## 产生模拟数据
        x_hat.shape = -1,1
        # 数据预测
        y_hat = model.predict(x_hat)
        # 计算准确率
        s = model.score(x, y)

        # 当d等于5的时候，设置为N-1层，其它设置0层；将d=5的这条线凸显出来
        z = N + 1 if (d == 5) else 0
        label = u'%d阶, 正确率=%.3f' % (d,s)
        plt.plot(x_hat, y_hat, color=colors[i], lw=2, alpha=0.75, label=label, zorder=z)
    
    plt.legend(loc = 'upper left')
    plt.grid(True)
    plt.title(titles[t])
    plt.xlabel('X', fontsize=16)
    plt.ylabel('Y', fontsize=16)
plt.tight_layout(1, rect=(0,0,1,0.95))
plt.suptitle(u'各种不同线性回归过拟合显示', fontsize=22)
plt.show()

梯度下降求解θ参数

梯度方向

 

批量梯度下降算法(BGD)

import time
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False

# 1. 随机数据的产生
np.random.seed(1)
n = 1000
# 这里产生均值为-1，标准差为10的服从正太分布的随机数据
b_values = np.random.normal(loc=-1.0, scale=0.01, size=n)
c_values = np.random.normal(loc=0.0, scale=1.0, size=n)

def bgd(n, b_values, c_values, max_iter=1000, tol=0.00001, alpha=0.01, show_img=True):
    """
    计算当样本数目，对应的最小的y以及x的值
    :param b_values:
    :param c_values:
    :param max_iter:
    :param tol:
    :param alpha:
    :return:
    """
    # 检查一下数据是否是n组数据, assert的功能就是检查代码的运行，如果assert后的表达式执行为False，那么代码报错；如果执行为True，那么直接运行接下来的代码
    assert len(b_values) == n and len(c_values) == n

    def f(x, n, b_values, c_values):
        # 原始函数
        result = 0
        for i in range(n):
            # 当n足够大的时候，直接累加的这种方式可能会出现result溢出的情况，这里就为了解决这个文件，这里将累加更换为均值的操作
            y = x ** 2 + b_values[i] * x + c_values[i]
            result += y / n
        return result

    def h(x, n, b_values, c_values):
        # 原始函数对应的导函数
        result = 0
        for i in range(n):
            # 当n足够大的时候，直接累加的这种方式可能会出现result溢出的情况，这里就为了解决这个文件，这里将累加更换为均值的操作
            y = 2 * x + b_values[i]
            result += y / n
        return result

    step_channge = 1.0 + tol
    step = 0
    # 随机的给定一个初始的x值
    current_x = np.random.randint(low=-10, high=10)
    current_y = f(current_x, n, b_values, c_values)
    print("当前的样本数据为:{}".format(n))
    print("b的均值为:{}".format(np.mean(b_values)))
    print("c的均值为:{}".format(np.mean(c_values)))

    # 画图用相关变量
    y_values = []
    if show_img:
        y_values.append(current_y)
    y_change_values = []

    # 开始迭代
    t1 = time.time()
    while step_channge > tol and step < max_iter:
        # 1. 计算函数的梯度值
        current_df = h(current_x, n, b_values, c_values)
        # 2. 基于梯度下降更新模型参数(更新x)
        current_x = current_x - alpha * current_df
        # 3. 更新x所对应的y值
        tmp_y = f(current_x, n, b_values, c_values)
        # 4. 记录一下变化大小
        step_channge = np.abs(current_y - tmp_y)
        step += 1
        current_y = tmp_y

        # 5. 添加可视化内容
        if show_img:
            y_values.append(current_y)
            y_change_values.append(step_channge)

    t2 = time.time()
    print("BGD的执行时间:{}".format(t2 - t1))
    print("最终更新的次数为:{}, 最终的变化率为:{}".format(step, step_channge))
    print("最终的结果为:{}---->{}".format(current_x, current_y))

    # 可视化
    if show_img:
        plt.figure(facecolor='w')
        plt.subplot(1, 2, 1)
        plt.plot(range(step), y_change_values, 'r-')
        plt.xlabel('迭代次数')
        plt.ylabel('每次迭代的y值变化大小')
        plt.subplot(1, 2, 2)
        plt.plot(range(step + 1), y_values, 'g-')
        plt.xlabel('迭代次数')
        plt.ylabel('函数值')
        plt.suptitle('BGD的变化情况')
        plt.show()

随机梯度下降算法(SGD)

 

import time
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False

# 1. 随机数据的产生
np.random.seed(1)
n = 1000
# 这里产生均值为-1，标准差为10的服从正太分布的随机数据
b_values = np.random.normal(loc=-1.0, scale=0.01, size=n)
c_values = np.random.normal(loc=0.0, scale=1.0, size=n)

def sgd(n, b_values, c_values, max_iter=10, tol=0.0000001, alpha=0.01, show_img=True):
    """
    计算当样本数目，对应的最小的y以及x的值
    :param b_values:
    :param c_values:
    :param max_iter:
    :param tol:
    :param alpha:
    :return:
    """
    # 检查一下数据是否是n组数据, assert的功能就是检查代码的运行，如果assert后的表达式执行为False，那么代码报错；如果执行为True，那么直接运行接下来的代码
    assert len(b_values) == n and len(c_values) == n

    def f1(x, b, c):
        return x ** 2 + b * x + c

    def f(x, n, b_values, c_values):
        # 原始函数
        result = 0
        for i in range(n):
            # 当n足够大的时候，直接累加的这种方式可能会出现result溢出的情况，这里就为了解决这个文件，这里将累加更换为均值的操作
            y = f1(x, b_values[i], c_values[i])
            result += y / n
        return result

    def h1(x, b, c):
        return x * 2 + b

    step_channge = 1.0 + tol
    step = 0
    # 随机的给定一个初始的x值
    current_x = np.random.randint(low=-10, high=10)
    current_y = f(current_x, n, b_values, c_values)
    print("当前的样本数据为:{}".format(n))
    print("b的均值为:{}".format(np.mean(b_values)))
    print("c的均值为:{}".format(np.mean(c_values)))

    # 画图用相关变量
    change_numbers = 0
    y_values = []
    if show_img:
        y_values.append(current_y)
    y_change_values = []

    # 开始迭代
    t1 = time.time()
    while step_channge > tol and step < max_iter:
        """
        在随机梯度中，是使用每条样本更新一次模型参数(这里的模型参数就是x)；这里总共有：n条样本，那也就是说在一次epoch中，更新参数n次。
        """
        print(step)
        # 将更新的样本顺序随机化
        random_index = np.random.permutation(n)
        for index in random_index:
            # 使用index样本更新模型参数
            # 1. 计算函数的梯度值
            current_df = h1(current_x, b_values[index], c_values[index])
            # 2. 基于梯度下降更新模型参数(更新x)
            current_x = current_x - alpha * current_df
            # 3. 更新x所对应的y值
            tmp_y = f(current_x, n, b_values, c_values)
            # 4. 记录一下变化大小
            step_channge = np.abs(current_y - tmp_y)
            current_y = tmp_y
            change_numbers += 1

            # 5. 添加可视化内容
            if show_img:
                y_values.append(current_y)
                y_change_values.append(step_channge)

            # 如果模型效果不错的情况下，退出更新
            if step_channge < tol:
                break
        # 更新迭代次数
        step += 1

    t2 = time.time()
    print("SGD的执行时间:{}".format(t2 - t1))
    print("最终更新的次数为:{}, 参数的更新次数:{}, 最终的变化率为:{}".format(step, change_numbers, step_channge))
    print("最终的结果为:{}---->{}".format(current_x, current_y))

    # 可视化
    if show_img:
        plt.figure(facecolor='w')
        plt.subplot(1, 2, 1)
        plt.plot(range(change_numbers), y_change_values, 'r-')
        plt.xlabel('迭代次数')
        plt.ylabel('每次迭代的y值变化大小')
        plt.subplot(1, 2, 2)
        plt.plot(range(change_numbers + 1), y_values, 'g-')
        plt.xlabel('迭代次数')
        plt.ylabel('函数值')
        plt.suptitle('SGD的变化情况')
        plt.show()

SGD速度比BGD快(迭代次数少).

SGD在某些情况下(全局存在多个相对最优解/J(θ)不是一个二次)，SGD有可能跳 出某些小的局部最优解，所以不会比BGD坏；SGD在收敛的位置会存在J(θ)函数 波动的情况。

BGD一定能够得到一个局部最优解(在线性回归模型中一定是得到一个全局最优 解)，SGD由于随机性的存在可能导致最终结果比BGD的差。

注意：优先选择SGD

小批量梯度下降法(MBGD)

如果即需要保证算法的训练过程比较快，又需要保证最终参数训练的准确率，而 这正是小批量梯度下降法（Mini-batch Gradient Descent，简称MBGD）的初 衷。MBGD中不是每拿一个样本就更新一次梯度，而且拿b个样本(b一般为10)的 平均梯度作为更新方向。

梯度下降法调优策略

由于梯度下降法中负梯度方向作为变量的变化方向，所以有可能导 致最终求解的值是局部最优解，所以在使用梯度下降的时候，一般需 要进行一些调优策略：

• 学习率的选择：学习率过大，表示每次迭代更新的时候变化比较大，有可能 会跳过最优解；学习率过小，表示每次迭代更新的时候变化比较小，就会导 致迭代速度过慢，很长时间都不能结束；
• 算法初始参数值的选择：初始值不同，最终获得的最小值也有可能不同，因 为梯度下降法求解的是局部最优解，所以一般情况下，选择多次不同初始值 运行算法，并最终返回损失函数最小情况下的结果值；
• 标准化：由于样本不同特征的取值范围不同，可能会导致在各个不同参数上 迭代速度不同，为了减少特征取值的影响，可以将特征进行标准化操作。

 BGD、SGD、MBGD的区别：

• 当样本量为m的时候，每次迭代BGD算法中对于参数值更新一次，SGD算法 中对于参数值更新m次，MBGD算法中对于参数值更新m/n次，相对来讲 SGD算法的更新速度最快；
• SGD算法中对于每个样本都需要更新参数值，当样本值不太正常的时候，就 有可能会导致本次的参数更新会产生相反的影响，也就是说SGD算法的结果 并不是完全收敛的，而是在收敛结果处波动的；
• SGD算法是每个样本都更新一次参数值，所以SGD算法特别适合样本数据量 大的情况以及在线机器学习(Online ML)。

线性回归的扩展

线性回归针对的是θ而言是一种，对于样本本身而言，样本可以是非线性的 也就是说最终得到的函数f:x->y；函数f(x)可以是非线性的，比如：曲线等。

局部加权回归-损失函数

w(i)是权重，它根据要预测的点与数据集中的点的距离来为数据集中的点赋权值。 当某点离要预测的点越远，其权重越小，否则越大。常用值选择公式为：

该函数称为指数衰减函数，其中k为波长参数，它控制了权值随距离下降的速率。

注意：使用该方式主要应用到样本之间的相似性考虑，主要内容在SVM中再考虑 (核函数)。

回归分类

Logistic回归算法

 

 

极大似然函数的随机梯度

θ参数求解

 

极大似然估计与Logistic回归损失函数

 

import pandas as pd
import numpy as np
import matplotlib as mpl

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import recall_score, precision_score, accuracy_score

# 设置字符集，防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False

# 一、加载数据
names = ['id', 'Clump Thickness', 'Uniformity of Cell Size', 'Uniformity of Cell Shape',
         'Marginal Adhesion', 'Single Epithelial Cell Size', 'Bare Nuclei',
         'Bland Chromatin', 'Normal Nucleoli', 'Mitoses', 'Class']
path = '../datas/breast-cancer-wisconsin.data'
df = pd.read_csv(path, sep=',', header=None, names=names)
# df.info()

# 二、数据的清洗
"""
可选操作
比如异常数据的处理，缺省数据的填充....
"""
# inplace: 该参数的作用是给定，是否是在原始的DataFrame上做替换，设置为true表示在原始DataFrame上做替换，直接修改原始的DataFrame的数据格式；默认是False
# df = df.replace('?', np.nan)
df.replace('?', np.nan, inplace=True)
# 删除所有具有非数字的样本
# how: 指定如何删除数据，如果设置为any，默认值, 那么只要行或者列中存在nan值，那么就进行删除操作；如果设置为all，那边要求行或者列中的所有值均为nan的时候，才可以进行删除操作
df.dropna(axis=0, how='any', inplace=True)

# 三、基于业务提取最原始的特征属性X和目标属性Y
X = df[names[1:10]]
# X.info()
Y = df[names[10]]
print("Y的取值为:{}".format(np.unique(Y)))
print("原始X形状:{}， 原始Y形状:{}".format(X.shape, Y.shape))

# 四、数据的划分(将数据划分为训练集和测试集)
# train_size：给定划分中的训练集占比，一般情况下，我们训练集和测试集的占比一般为8:2，如果样本比较多的时候，可以考虑7:3，如果样本比较少的时候，可以考虑9:1
# test_size: 给定划分中的测试集的占比，其中train_size和test_size只能有且给定其中的一个。
x_train, x_test, y_train, y_test = train_test_split(X, Y, train_size=0.8, random_state=28)
# print("训练集样本形状:{}, 测试集样本形状:{}".format(x_train.shape, x_test.shape))


# 五、特征工程
"""
可选操作
主要就是一些：哑编码、TF-IDF、连续数据的离散化、标准化、归一化、特征选择、降维....
NOTE: 所有的特征工程其实和算法模型训练一样，特征工程的意思就是对特征属性做一个数据的转换操作，所以需要获取一个对应数据转换函数---> 需要从训练数据中训练出来
"""
"""
NOTE:
sklearn中，所有特征工程、算法模型的API基本类似，主要是以下几个API：
fit: 模型/函数训练，基于给定的训练集数据，训练模型对象
transform: 使用训练好的模型对给定的数据集做一个对应的数据转换(使用训练好的函数转换数据)，一般出现在特征工程的相关的类中间
predict：使用构建好的算法模型对样本数据做一个预测的操作，其实内部也就是使用训练好的函数转换特征属性到目标属性，一般出现在算法模型的相关的类中间
fit_transform: 也就是fit和transform的结合，底层先使用给定的训练数据训练模型，找出对应的转换函数或者模型对象，然后使用找到的函数或者模型对训练数据做一个转换操作
"""
# 做一下数据的标准化操作
"""
实现方式一：
# a. 创建标准化的对象
ss = StandardScaler()
# b. 训练标准化的转换函数
ss.fit(x_train, y_train)
# c. 使用训练好的标注化模型对数据做一个转换操作
x_train = ss.transform(x_train)
x_test = ss.transform(x_test)

"""
# 实现方式二：
# a. 创建标准化的对象
ss = StandardScaler()
# b. 训练标准化的转换函数并使用训练好的函数对训练数据做一个转换操作
x_train = ss.fit_transform(x_train, y_train)
# c. 使用训练好的标注化模型对数据做一个转换操作
x_test = ss.transform(x_test)

# 六、算法模型的选择/算法模型对象的构建
algo = LogisticRegression()

# 七、算法模型的训练
algo.fit(x_train, y_train)

# 7.1 查看训练好的模型参数
print("Logistic回归的各个特征属性对应的权重参数θ:{}".format(algo.coef_))
print("Logistic的截距项的值:{}".format(algo.intercept_))

# 八、模型效果评估
y_hat = algo.predict(x_test)
print("在训练集上的模型效果(分类算法中为准确率):{}".format(algo.score(x_train, y_train)))
print("在测试集上的模型效果(分类算法中为准确率):{}".format(algo.score(x_test, y_test)))
print("在测试集上的召回率的值:{}".format(recall_score(y_true=y_test, y_pred=y_hat, pos_label=4)))

# 九、输出分类中特有的一些API
print("=" * 100)
print("预测值:{}".format(algo.predict(x_test)))
print("=" * 100)
print("决策函数输出值(在Logistic回归中就是θx的值):{}".format(algo.decision_function(x_test)))
print("=" * 100)
print("属于各个类别的概率值:{}".format(algo.predict_proba(x_test)))
print("=" * 100)
print("属于各个类别的log转换之后的概率值:{}".format(algo.predict_log_proba(x_test)))
# 更新阈值，产生预测结果, 概率大于0.3，属于4，小于等于0.3属于2
a = 0.2
y_predict_proba = algo.predict_proba(x_test)[:, 1]
y_predict_proba[y_predict_proba > a] = 4
y_predict_proba[y_predict_proba <= a] = 2
# pos_label： 指定计算召回率和精确率的时候，正例的y值是等于多少的，默认为0
print("在测试集上的召回率和精确率的值（修改阈值）:{},{},{}".format(
    recall_score(y_true=y_test, y_pred=y_predict_proba, pos_label=4),
    precision_score(y_true=y_test, y_pred=y_predict_proba, pos_label=4),
    accuracy_score(y_true=y_test, y_pred=y_predict_proba)
))

Softmax回归算法

 softmax回归是logistic回归的一般化，适用于K分类的问题，第k类的参数为 向量θk，组成的二维矩阵为θk*n；

softmax函数的本质就是将一个K维的任意实数向量压缩（映射）成另一个K维 的实数向量，其中向量中的每个元素取值都介于（0，1）之间。

softmax回归概率函数为：

Softmax算法原理

Softmax算法损失函数

Softmax算法梯度下降法求解

 

sklearn中没有softmax算法api，依然是使用Logistic，只不过需要把参数multi_class设为：multi_class = 'multinomial'。

总结

 线性模型一般用于回归问题，Logistic和Softmax模型一般用于分类问题。

求θ的主要方式是梯度下降算法，梯度下降算法是参数优化的重要手段，主要 是SGD，适用于在线学习以及跳出局部极小值。

Logistic/Softmax回归是实践中解决分类问题的最重要的方法。

广义线性模型对样本要求不必要服从正态分布、只需要服从指数分布簇(二项 分布、泊松分布、伯努利分布、指数分布等)即可；广义线性模型的自变量可 以是连续的也可以是离散的。