使用Python進(jìn)行線性回歸

線性回歸是最簡(jiǎn)單同時(shí)也是最常用的一個(gè)統(tǒng)計(jì)模型。線性回歸具有結(jié)果易于理解，計(jì)算量小等優(yōu)點(diǎn)。如果一個(gè)簡(jiǎn)單的線性回歸就能取得非常不錯(cuò)的預(yù)測(cè)效果，那么就沒(méi)有必要采用復(fù)雜精深的模型了。
今天，我們一起來(lái)學(xué)習(xí)使用Python實(shí)現(xiàn)線性回歸的幾種方法：
    通過(guò)公式編寫矩陣運(yùn)算程序；
    通過(guò)使用機(jī)器學(xué)習(xí)庫(kù)sklearn;
    通過(guò)使用statmodels庫(kù)。
這里，先由簡(jiǎn)至繁，先使用sklearn實(shí)現(xiàn)，再講解矩陣推導(dǎo)實(shí)現(xiàn)。
1.使用scikit-learn進(jìn)行線性回歸
設(shè)置工作路徑
#
import os
os.getcwd()
os.chdir('D:\my_python_workfile\Project\Writting')

加載擴(kuò)展包

import pandas as pd
import numpy as np
import pylab as pl
import matplotlib.pyplot as plt

載入數(shù)據(jù)并可視化分析

這里，為了簡(jiǎn)單起見，使用sklearn中自帶的數(shù)據(jù)集鳶尾花數(shù)據(jù)iris進(jìn)行分析，探索『花瓣寬』和『花瓣長(zhǎng)』之間的線性關(guān)系。

from sklearn.datasets import load_iris
# load data
iris = load_iris()
# Define a DataFrame
df = pd.DataFrame(iris.data, columns = iris.feature_names)
# take a look
df.head()
#len(df)

# correlation
df.corr()

  

# rename the column name

df.columns = ['sepal_length','sepal_width','petal_length','petal_width']
df.columns

Index([u'sepal_length', u'sepal_width', u'petal_length', u'petal_width'], dtype='object')

plt.matshow(df.corr())

由上面分析可知，花瓣長(zhǎng)sepal length和花瓣寬septal width有著非常顯著的相關(guān)性。
下面，通過(guò)線性回歸進(jìn)一步進(jìn)行驗(yàn)證。
# save image
fig,ax = plt.subplots(nrows = 1, ncols = 1)
ax.matshow(df.corr())
fig.savefig('./image/iris_corr.png')

建立線性回歸模型

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

lr = LinearRegression()
X = df[['petal_length']]
y = df['petal_width']
lr.fit(X,y)
# print the result
lr.intercept_,lr.coef_

(-0.3665140452167297, array([ 0.41641913]))

# get y-hat
yhat = lr.predict(X = df[['petal_length']])

# MSE
mean_squared_error(df['petal_width'],yhat)

# lm plot
plt.scatter(df['petal_length'],df['petal_width'])
plt.plot(df['petal_length'],yhat)

#save image
plt.savefig('./image/iris_lm_fit.png')

2.使用statmodels庫(kù)

#import statsmodels.api as sm
import statsmodels.formula.api as sm

linear_model = sm.OLS(y,X)

results = linear_model.fit()

results.summary()

OLS Regression Results 

3.使用公式推導(dǎo)
線性回歸，即是使得如下目標(biāo)函數(shù)最小化：

使用最小二乘法，不難得到β的估計(jì)：

從而，我們可以根據(jù)此公式，編寫求解β^的函數(shù)。
from numpy import *

#########################
# 定義相應(yīng)的函數(shù)進(jìn)行矩陣運(yùn)算求解。
def standRegres(xArr, yArr):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    xTx = xMat.T * xMat
    if linalg.det(xTx) == 0.0:
        print "this matrix is singular, cannot do inverse!"
        return NA
    else :
        ws = xTx.I * (xMat.T * yMat)
        return ws

# test
x0 = np.ones((150,1))
x0 = pd.DataFrame(x0)
X0 = pd.concat([x0,X],axis  = 1)
standRegres(X0,y)

matrix([[-0.36651405],
        [ 0.41641913]])

結(jié)果一致。