線性迴歸-最小二乘方法程式碼實現
線性迴歸-最小二乘方法
使用最小二乘的方法進行原始的計算方式編寫
先把該匯入的包全部匯入了
# 首先需要匯入對應的包
import pandas as pd # 資料處理
import numpy as np # 資料計算
import matplotlib.pyplot as plt # 畫圖
from sklearn.model_selection import train_test_split # 專門用於資料分割的函式,分割訓練集和測試集
# 因為畫圖,所以先設定字符集防止中文無法顯示的問題
plt.rcParams['font.sans-serif' 1 讀取資料
資料儲存在檔案中需要讀取,現在使用的資料是使用者的不同時間段家用電量的表格
dir = './data/household_power_consumption_1000.txt'
# 使用pandas進行資料匯入
# 使用的是;進行分隔,資料中可能存在大量對的混合型別,所以需要進行low_menorey設定
data_source = pd.read_csv(dir, sep=';', low_memory=False)
資料匯入完成後檢視下資料,大致分清需要哪些資料引數計算,由於是估計功率和電流之間的關係,所以取出資料的時候就需要取出對應的資料即可
2 取出資料
# 根據資料的內容而言,資料結構如下
# 日期、時間、有功功率、無功功率、電壓、電流、廚房用電功率、洗衣服用電功率、熱水器用電功率
# 尋找功率和電流之間的關係
# 這裡就需要對資料(可以看做是一個二維矩陣)進行切片選出需要的資料
X = data_source.iloc[:, 2:4] # 切片不包括最後一項
Y = data_source.iloc[:, 5]
pandas的資料結構中可以直接使用iloc進行切片,其切片的規則和numpy中的規則非常相似,因為資料本身是一個表格,這個表格有行有列,只有兩個資料維度,那麼每個樣本存在多個數據,每個屬性就是一列
3 分割資料
這裡採用的是隨機劃分,將資料劃分成為訓練集:測試集=4:1的形式,使用隨機抽樣的方法,當然分割資料肯定不止是這一種方式,這裡先使用這種隨機劃分的方式
# 劃分資料,test_size是按照比例,random_state是表示隨機種子
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0)
4 資料變化
直接的pandas資料不好計算,轉成numpy資料進行處理,之後使用最小二乘的時候對矩陣的求逆矩陣,求轉置等等,方便計算
# 建立兩個訓練資料的關係
x_train = np.mat(X_train)
y_train = np.mat(Y_train).reshape(-1,1) # 因為y實際上搞出來不是一個列向量,這裡進行轉一把
這裡有個坑,如果本身的pandas就是一列的話這裡直接轉回造成資料的矩陣型狀無法參與後面的計算,所以這裡需要轉變一下
5 計算出係數矩陣
theta的求解利用最小二乘的方式,根據理論,中間的計算推導步驟很麻煩,要用到矩陣求導,大量的矩陣變化,好在我們直接有最後的公式
這裡已經有X和Y了,直接帶進去計算
theta = (x_train.T * x_train).I * x_train.T * y_train
引數矩陣就已經得到
6 利用引數矩陣計算預測結果
引數矩陣已經得到後,帶入原來的y=kx+b這種線性方程直接可以計算出預測值,而測試集資料需要用x_test去和引數矩陣計算,計算的結果與y_test進行比較
# 再帶入theta進行驗證 使用測試集合獲得預測的數值
the_predict = np.mat(X_test) * theta
7 畫出對比圖
為了直觀的看到整個圖,這裡直接畫出圖,但是畫圖需要注意的是x軸的取值,因為資料本身是對比測試值和預測值的關係,那麼這兩個的樣本都是一一對應的,將樣本取值當做x軸顯然不現實,每個樣本本身就具有多個屬性,因此,這裡使用個數來當做x軸取值
# 接下來畫圖
# 建立圖布
fig = plt.figure()
# 設定標題
plt.title('功率和電流的預測與實際對比圖')
# 設定x軸和y軸
# 這裡因為資料樣品是離散的隨機取樣,所以x軸定為點的序號,比如第一個點第二個點
x = np.arange(len(X_test)) # 總共的樣本個數
plt.plot(x, Y_test, 'r-',label='真實值',) # 紅色的為真實值
plt.plot(x, the_predict, 'g-',label='預測值',) # 綠色的為預測值
plt.legend(loc = 'lower right') # 設定圖例顯示的位置
plt.show()
以上是最原始的最小二乘法,中間的模型是通過自己計算的到的結論,但是我們的前輩們早就已經計算好了對應的theta公式,直接用就可以了,但是這裡還存在特例,就是X轉置乘以X如果不是正定矩陣,則還不能這麼計算
使用模型導包的方式進行預測
下面是另外的方式進行計算用電的功率和時間之間的關係
在實際專案中,線性迴歸模型實際上已經是封裝好的模型
依然開始各種匯入
# coding=gbk
# 因為裡面存在中文,所以這裡直接新增後就不報錯了
# 匯入對應的庫
import numpy as np # 資料計算
import matplotlib.pyplot as plt # 畫圖
import pandas as pd # 資料處理
from sklearn.linear_model import LinearRegression # 線性迴歸模型
from sklearn.model_selection import train_test_split # 構建資料分割器
from sklearn.preprocessing import StandardScaler # 歸一化資料的模型
# 配置對應的中文顯示問題
plt.rcParams['font.sans-serif'] = [u'simHei']
plt.rcParams['axes.unicode_minus'] = False
這次因為裡面註釋寫的多,不知道哪裡有問題,反正不加第一個就報錯,既然網上查了說加上就不撥錯,親測有效,所以沒有為什麼,加上就好了
1 匯入資料
# 匯入檔案
data = pd.read_csv('./data/household_power_consumption_1000.txt', sep=';', low_memory=False)
2 擷取需要的資料
這裡考慮的是時間和功率之間的關係,但是時間本身是一個字串,在資料中存在兩列用來記錄時間
這裡需要將這兩列單獨取出來
data.iloc[:, :2]
依然使用的是切片的方式
取出來之後並不能直接參與計算,因為字串是需要進行轉換才能進行計算,時間有年月日時分秒六個單位組成,資料中確實也是記錄下了這六個單位,因此需要進行對應的處理把時間分解出來,使用6個屬性進行描述時間
將字串轉變為時間的函式
time.striptime(str, format)
這裡需要得到str,也就是時間字串,原來的時間字串是儲存在表中的兩列,雖然能夠切片進行切出來,但是實際上還是一個具有很多行的列表,這裡使用
pandas_data.apply(function, axis=1, args, kwargs)
這個函式的意義在於進行轉換,把pandas資料傳遞到function函式中進行處理,axis=1表示已一行一行的傳進去,很顯然這裡的function就是一個函式,但是這個函式不能直接在後面加上()傳參呼叫,這裡只能寫函式名
所以不管這個函式本身需要幾個引數,這裡預設的第一個引數就是pandas資料的一行,以列表的方式傳進去,那麼如果這個函式本身還需要其他引數的話就通過args和kwargs進行傳入,這兩個引數是留給function的,因為只寫函式名沒法指明需要的引數,所以這裡留下了兩個介面
既然apply把時間對應的兩列拆成一行一行的傳進去,那麼就應該由function這個函式來接收,傳進去的是時間字串列表,所以需要自定義函式處理
def date_reduce(dt):
import time
dt = time.striptime(' '.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)
這裡得到傳進來的列表後進行拆分,把列表轉成可以用的時間字串,之後使用格式化字串對原時間字串轉義,返回的是一個時間物件
前面說過需要將時間資料拆分成年月日時分秒六個屬性,所以這裡直接返回這六個屬性對應的值
apply函式非常漂亮的將返回的接收了這六個屬性值,但是它依然不是正確的能夠使用的資料各式,因此這裡使用
pd.Serires(data)
使用pd的資料型別Serires直接強行將返回的資料轉成Serires類,這就很像str()這種方式強行轉成字串一樣
因此整個時間的處理如果下
X = data.iloc[:, :2].apply(lambda x: pd.Series(date_reduce(x)), axis=1)
得到了可以用於計算的X樣本資料
下面再把需要的Y值取出來
Y = data.iloc[:, 2] # 就當做是numpy一樣直接取
3 使用拆分函式直接將資料集拆分稱為訓練集和測試集
# 分成訓練集和測試集
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0)
但是這裡只要檢視一下就會發現這個資料根本不能使用,X中存在六個屬性,而且使用方法進行檢視
print(X_train.describe())
你會發現這六個屬性中有三個屬性的std(標準差)為0,還有兩個標準差很大,標準差為0說明資料本身沒有變動,那就意味著整個一列都是完全相同的,這樣的資料在前面增加任何係數都沒有意義,標準差很大的資料表明波動性太大,因此前面的係數確定就很不均衡,容易給一些屬性過高的權重而使得某些屬性根本就沒有權重
為了解決這些問題,必須引入歸一化操作
4 歸一化資料
歸一化是利用模型將資料進行壓縮到標準差為1的區間,這樣使得每個有用的屬性都把標準差壓縮到1,大家權重相同,這樣訓練出來的模型可靠度要高一些
ss = StandardScaler()
直接匯入標準化模型物件,此時整個模型還是空的,還沒有給他資料進行訓練,將訓練資料扔進這個模型中,其訓練出來的資料就是經過標準化以後的資料了,並且將針對這些資料的標準化計算方法封裝在ss物件上
這樣做的好處在於如果後面訓練出來的預測模型效果還不錯,那麼就把這個模型和預測模型都儲存起來,以後來的新資料都會先進過ss進行標準化後再扔進預測模型進行計算
歸一化訓練資料和測試資料
X_train = ss.fit_transform(X_train)
X_test = ss.transform(X_test)
這裡上面和下面不一樣,這裡包含fit的方法存在著訓練的含義,單純的transform表示你轉變,上面已經fit一次表示已經訓練好了x_train資料然後再_transform,表示再把訓練資料轉變成訓練好的資料返回.
下面的X_test的意思就是直接對x_test進行轉變返回,這裡不能再訓練,因為一旦再訓練就會給覆蓋原來的訓練模型ss,那麼X_test與X_train又不是完全相同,這就會回使得獲得X_train資料時候的ss模型與獲得X_test資料的ss模型不載相同,後面一定會存在誤差.所以這裡訓練過的模型就不再訓練了
5 使用最小二乘進行計算theta值
這裡當然可以對資料進行numpy轉變,使用我們的推導結論公式帶入直接計算.這樣做沒有任何問題
但是這裡我們使用別人已經寫好的模型
model_train = LinearRegression(fit_intercept=True) # 建立模型物件
同ss一樣,這個只是建立了模型物件,還沒有進行訓練
model_train.fit(X_train, Y_train)
帶入資料進行訓練,這裡只有fit沒有後面的transform,表示這裡只是把資料扔進去進行模型訓練,不需要返回值,我們需要的是模型,執行完這一步之後模型就已經訓練完成了,theta什麼的已經整合到model_train這個物件身上了,後面只需要呼叫這個物件就行了
Y_predict = model_train.predict(X_test)
使用這個模型直接把測試資料扔進去進行計算,並把計算後的值反回出來,因為這裡扔進去的是測試資料,所以計算結果就是預測值,看吧這裡別人直接就算好了
現在我們手裡面有了預測資料和測試資料的真實值,我們也不知道準不準確,不要驚慌,別人封裝好的功能在model_train上面也有
# 檢視模型的r方 也叫做準確率
print(f'訓練樣本的R**2(準確率):{model_train.score(X_train,Y_train)}')
print(f'測試樣本的R**2(準確率):{model_train.score(X_test,Y_test)}')
print(f'模型引數θ:{model_train.coef_}') # coef_可以直接獲取到訓練後的引數
print(f'模型截距:{model_train.intercept_}') # 線性迴歸y=kx+b的b就是截距
直接可以查看準確率,當然這個模型訓練的不咋樣,準確率只有0.22.這不是重點,重點是這裡直接可以通過模型.score(輸入,輸出)進行查詢看看模型到底效果如何,同時還可以看到引數和截距
y=kx+b
k找到了,b找到了,還有模型的正確率,基本上這次建模就已經完成了
6 畫圖
當然如果想更直觀的看到模型效果如何,可以直接把他們畫圖畫出來
fig = plt.figure()
# 定義x軸
x = np.arange(len(X_test))
plt.plot(x, Y_test, 'r-', label='真實值')
plt.plot(x, Y_predict, 'g-', label='預測值')
plt.legend(loc='upper left')
plt.title('以時間的自變數預測功率與時間的關係模型')
plt.show()
7 模型儲存
模型已經訓練完成,如果效果還不錯的話則需要把模型儲存起來方便以後有了新的資料直接呼叫這個模型進行預測
這裡就需要想清楚需要儲存哪些東西,首先ss肯定要存,你來的新資料難道就不歸一化麼?不歸一化的話帶進去的值就不滿足預測模型需要的值,所以ss模型一定要存
用腳指頭想也應該想到需要存預測模型,費了半天勁就是為了搞到預測模型,不存它存誰
from sklearn.externals import joblib
joblib.dump(ss, './standard_to_lr_data_model.model') # ss = StandardScaler() 標準化類:已經訓練過的
joblib.dump(model_train, './lr_predict_model.model') # 就是匯入的訓練物件model_train = LinearRegression(fit_intercept=True):已經訓練過的
最後給總程式碼
# coding=gbk
# 因為裡面存在中文,所以這裡直接新增後就不報錯了
# 匯入對應的庫
import numpy as np # 資料計算
import matplotlib.pyplot as plt # 畫圖
import pandas as pd # 資料處理
from sklearn.linear_model import LinearRegression # 線性迴歸模型
from sklearn.model_selection import train_test_split # 構建資料分割器
from sklearn.preprocessing import StandardScaler # 歸一化資料的模型
def date_reduce(dt):
# 傳進來的是兩列資料,這兩列資料是時間,實際是一個列表兩個元素,用空格拼接稱為字串
import time
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)
# 配置對應的中文顯示問題
plt.rcParams['font.sans-serif'] = [u'simHei']
plt.rcParams['axes.unicode_minus'] = False
# 匯入檔案
data = pd.read_csv('./data/household_power_consumption_1000.txt', sep=';', low_memory=False)
# 因為需要處理時間這個變數,其對應的是前兩列做對應的變化
X = data.iloc[:, :2].apply(lambda x: pd.Series(date_reduce(x)), axis=1)
# 值得一提的是apply函式的使用,apply函式的引數,第一個引數是一個函式,這個引數是需要傳遞一個函式,可以是一個隱函式,這裡如果不是隱函式則傳遞的函式自動的會根據後面引數的決定傳遞的物件,第二個引數是axis,表示軸,如果為1,則表示一行一行的傳進去,那麼誰呼叫就傳遞誰,後面的args和kwargs都是設定傳遞給第一個引數的,因為第一個引數是寫函式名,寫函式名沒有辦法直接執行,所以設計args和kwargs進行傳參
# pd.Series()是將裡面的內容轉變成為Series資料型別
# 處理好時間後繼續獲取對應的Y
Y = data.iloc[:, 2] # 就當做是numpy一樣直接取
# 分成訓練集和測試集
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0)
# print(X_train.describe())
# 查詢的裡面存在一個std,這個是標準差,平方後就是方差,查詢後發現裡面存在std為0的資料,標準差為0,表示資料根本就沒有傳送變化
# std大了也表示波動很大,對於方差大的,需要規範化到0-1之間,最好將所有的資料都歸一化到0-1之間,這樣才不會導致造成引數非常大
# 標準化資料
ss = StandardScaler() # 建立標準化物件,使用標準差為1的
# print(X_train)
X_train = ss.fit_transform(X_train) # 轉變為了可以訓練用的資料
# print(X_train)
# 寫出對應的測試資料
X_test = ss.transform(X_test) # 媽蛋,感覺跟fit_transform沒有啥子區別
# print(X_test)
# 使用DataFrame將資料處理稱為列表形式
# print(X_train.describe()) # 這裡會報錯,說的是numpy.ndarray沒有對應的屬性
# 這就是說明X_train通過StandardScaler的fit_transform之後稱為了numpy陣列物件
# print(pd.DataFrame(X_train).describe())
# print(data.info())
# 直接訓練模型
model_train = LinearRegression(fit_intercept=True) # 建立模型物件
# 訓練模型
model_train.fit(X_train, Y_train)
# 獲得預測值
Y_predict = model_train.predict(X_test)
# 檢視模型的r方 也叫做準確率
print(f'訓練樣本的R**2(準確率):{model_train.score(X_train,Y_train)}')
print(f'測試樣本的R**2(準確率):{model_train.score(X_test,Y_test)}')
# 預測值減去實際值的平方後再求均值 就是均方差
mse = np.average((Y_predict - Y_test) ** 2)
# 如果再開根號 就是標準差
rmse = np.sqrt(mse)
print(f'模型引數θ:{model_train.coef_}') # coef_可以直接獲取到訓練後的引數
print(f'模型截距:{model_train.intercept_}') # 線性迴歸y=kx+b的b就是截距
# 兩個模型實際上已經完成了,現在可以進行畫圖操作
fig = plt.figure()
# 定義x軸
x = np.arange(len(X_test))
plt.plot(x, Y_test, 'r-', label='真實值')
plt.plot(x, Y_predict, 'g-', label='預測值')
plt.legend(loc='upper left')
plt.title('以時間的自變數預測功率與時間的關係模型')
plt.show()
# 到這裡為止,已經訓練好了兩個模型
# 1 標準化模型,這個模型是針對資料進行歸一化的模型,儲存它的意義在於如果我們最後的模型成功了,那麼新的資料進來後也需要執行歸一化,這就要求每次歸一化的模型必須使用的一樣,否則帶入計算模型中則無法得到預期的值
# 2 計算模型,這個模型就是我們的線性迴歸模型,這個模型裡面存放著針對歸一化之後的資料計算操作的封裝
# 使用固定的庫進行儲存模型
from sklearn.externals import joblib
joblib.dump(ss, './standard_to_lr_data_model.model') # ss = StandardScaler() 標準化類:已經訓練過的
joblib.dump(model_train, './lr_predict_model.model') # 就是匯入的訓練物件model_train = LinearRegression(fit_intercept=True):已經訓練過的