advcase1.py-20180704

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jul  3 17:43:53 2018

@author: vicky
"""

# 导入第三方包
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from sklearn.cross_validation import train_test_split
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import seaborn as sns

#数据集中各变量的描述性统计分析
#data = pd.read_csv('/Users/vicky/Desktop/Advertising.csv')
data.describe()

#对比三个因子与y的散点图
sns.pairplot(data, x_vars=["TV", "Radio", "Newspaper"], y_vars="Sales",size=5, aspect=0.7, kind="reg")
sns.plt.show()

#相关系数
data.corr()
#相关图
plt.imshow(data.corr(), cmap=plt.cm.Blues, interpolation='nearest')
plt.colorbar()
tick_marks = [i for i in range(len(data.columns))]
plt.xticks(tick_marks, data.columns, rotation='vertical')
plt.yticks(tick_marks, data.columns)

#训练集和测试集二八分
Train,Test = train_test_split(data, train_size = 0.8, random_state=1234)

#建线性回归模型
fit = smf.ols('Sales~TV+Radio+Newspaper', data = Train).fit()
fit.summary()

#去掉newspaper
fit2 = smf.ols('Sales~TV+Radio', data = Train.drop('Newspaper', axis = 1)).fit()
fit2.summary()

#加交互作用
fit3 = smf.ols('Sales~TV+Radio+TV:Radio', data = Train.drop('Newspaper', axis = 1)).fit()
fit3.summary()

pred = fit.predict(exog = Test)
pred2 = fit2.predict(exog = Test.drop('Newspaper', axis = 1))
pred3 = fit3.predict(exog = Test.drop('Newspaper', axis = 1))

#均方根误差
RMSE = np.sqrt(mean_squared_error(Test.Sales, pred))
RMSE2 = np.sqrt(mean_squared_error(Test.Sales, pred2))
RMSE3 = np.sqrt(mean_squared_error(Test.Sales, pred3))
print('RMES=%.4f\n' %RMSE)
print('RMES=%.4f\n' %RMSE2)
print('RMES=%.4f\n' %RMSE3)

#画真实值与预测值的对比图
plt.style.use('ggplot')
plt.scatter(Test.Sales, pred,c='b',label = 'Observations')
plt.plot([Test.Sales.min(), Test.Sales.max()], [pred.min(), pred.max()], 'r--', lw=2, label = 'Fitted line')
plt.title('Real Values VS. Predict Values')
plt.xlabel('Real Values')
plt.ylabel('Prediction Values')
plt.legend(loc = 'upper left')
plt.show()

#画残差图
plt.style.use('ggplot')
plt.scatter(pred,Test.Sales-pred3,c='b',label = 'Residuals')
plt.title('Residual Plot')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.legend(loc = 'upper left')
plt.show()

#两因子3d图
normal = np.array([0.1904,0.0456,-1])
point  = np.array([2.9004,0,0])
# a plane is a*x + b*y +c*z + d = 0
# [a,b,c] is the normal. Thus, we have to calculate
# d and we're set
d = -np.sum(point*normal) # dot product
# create x,y
x, y = np.meshgrid(range(50), range(300))
# calculate corresponding z
z = (-normal[0]*x - normal[1]*y - d)*1./normal[2]
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
fig.suptitle('Regression: Sales ~ Radio + TV')
ax = Axes3D(fig)
ax.set_xlabel('Radio')
ax.set_ylabel('TV')
ax.set_zlabel('Sales')
ax.scatter(data.Radio, data.TV, data.Sales, c='red')
ax.plot_surface(x,y,z, color='cyan', alpha=0.3)