1. 程式人生 > >機器學習2---線性模型

機器學習2---線性模型

 

LDA的程式碼可參見:https://blog.csdn.net/yt71656/article/details/45199603

 

來補充程式碼,萌新寫的,也沒有經過整理,只是為了熟悉思路,大神輕噴。

線性迴歸,python3。使用的資料是吳恩達的機器學習資料。

import numpy as np
import matplotlib.pyplot as plt

A = np.zeros((97,2),dtype=float)    #先建立一個 3x3的全零方陣A,並且資料的型別設定為float浮點型
 
f = open('\ex1data1.txt')               #開啟資料檔案檔案
lines = f.readlines()           #把全部資料檔案讀到一個列表lines中
A_row = 0                       #表示矩陣的行,從0行開始
for line in lines:              #把lines中的資料逐行讀取出來
    list = line.strip('\n').split(' ')      #處理逐行資料:strip表示把頭尾的'\n'去掉,split表示以空格來分割行資料,然後把處理後的行資料返回到list列表中
    A[A_row:] = list[0:2]                    #把處理後的資料放到方陣A中。list[0:3]表示列表的0,1,2列資料放到矩陣A中的A_row行
    A_row+=1                                #然後方陣A的下一行接著讀

x=A[:,0]
y=A[:,1]



# y = th1 + th2*x
th1=0
th2=1
N=20

m=len(A)
rl=0.02
sumt=[]  # 記錄代價函式取值
for i in range(N):    ## 優化10次
    sum=0
    sum1=0
    sum2=0
    for j in range(m):
        sum1=sum1+th1+th2*x[j]-y[j]   # th1
        sum2=sum2+(th1+th2*x[j]-y[j])*x[j]   # th2
        sum=sum+(th1+th2*x[j]-y[j])**2
    if (sum/m)<0.1:
        break
    sumt.append(sum/m)
        
    sum1=sum1/m
    sum2=sum2/m
    th1=th1-rl*sum1
    th2=th2-rl*sum2
    
print(th1,'+',th2,'x')

#繪製散點圖,examDf.jt為X軸,examDf.hk為Y軸
plt.scatter(x,y,color = 'darkgreen',label = "Exam Data")
 
#新增圖的標籤(x軸,y軸)
plt.xlabel("x")#設定X軸標籤
plt.ylabel("y")#設定Y軸標籤

x1 = np.linspace(0, 25, 100)
y1 = th1+th2*x1
plt.plot(x1, y1, color="blue")
plt.show()#顯示影象

plt.xlabel("n")#設定X軸標籤  次數
plt.ylabel("J")#設定Y軸標籤  代價函式
N1 = np.linspace(0, N, N)
plt.plot(N1, sumt, color="blue")
plt.show()#顯示影象

執行結果:

中間的藍色直線是擬合直線,其實看起來不太好,畢竟也只跑了20次。

代價函式的一個影象,橫座標是學習次數,可以看到是下降趨勢。

 

小tips:

(1)一定要動手啊~

(2)有時候其實不是演算法思路的問題,也不是程式碼的問題(簡潔性、優美之類的除外),就是取值的問題,多看看大神的程式碼真的有幫助。

多元線性迴歸分析的程式碼,資料是自己做的(3個輸入變數,一個輸出變數),這次是基於向量的思想在寫程式,比之前用迴圈的程式能簡潔些,哈哈哈~

# -*- coding: utf-8 -*-
"""
Created on Thu Nov  1 15:10:44 2018

@author: user
"""

import numpy as np
import matplotlib.pyplot as plt

# 生成隨機數
x=np.random.random([1,100])*5

x0=np.ones([1,100])
x1=x
x2=np.multiply(x,x)    # 矩陣元素做點乘
x3=np.multiply(x2,x)
x1=x1+np.random.random([1,100])-0.5
x2=x2+np.random.random([1,100])-0.5
x3=x3+np.random.random([1,100])-0.5
X=np.vstack((x0,x1,x2,x3))   # 按列合併,即將矩陣直接放在下面;np.hstack是按行合併
y=2*x0+3*x1+4*x2+5*x3    # m=100, n=3


# 求解 y=th0+th1*x1+th2*x2+th3*x3 中的引數th

# 初始化 
m=100 # 100個樣本
n=3  # 3個輸入變數(特徵)
ths=np.array([0,1,0,1]).reshape([4,1])
lr=0.00065
N=50  # 迭代次數
Co=[]
T=np.zeros([n+1,1])

# 梯度下降法訓練引數
for iter in range(N):
    cost=0
    for j in range(n+1):
        mm=X[j,:]
        mm=np.mat(mm)
        a=np.dot(ths.T,X)   # 矩陣做乘法
        b=np.dot(a,mm.T)
        c=np.dot(y,mm.T)
        T[j,0]=b-c
    ths=ths-lr*T/m
    for i in range(m):
        a=np.dot(ths.T,np.mat(X[:,i]).T)
        cost=cost+(a[0,0]-y[0,i])**2 
    Co.append(cost)
    
plt.xlabel("n")#設定X軸標籤  次數
plt.ylabel("J")#設定Y軸標籤  代價函式
N1 = np.linspace(0, N, N)
plt.plot(N1, Co, color="blue")
plt.show()#顯示影象
    
    

執行的代價函式結果:

感覺這個結果很佛性。。。其實引數結果還可以,設定初始值是個大學問。。。

 

小tip:

(1)矩陣的乘法要分清楚:*和dot() 不一樣。

(2)變數命名要看清楚,這個當時調了好久,後來發現是名字重複了。。。。。當時的心情,本來想去吃飯的。。。

(3)正則化,下次可以用上。

羅傑斯特迴歸

兩個程式碼,後期會修改比對:

import numpy as np
import math
import matplotlib.pyplot as plt

# 讀入資料
A = np.zeros((100,3),dtype=float)    
 
f = open('。。。')               #開啟資料檔案檔案
lines = f.readlines()           #把全部資料檔案讀到一個列表lines中
A_row = 0                       #表示矩陣的行,從0行開始
for line in lines:              #把lines中的資料逐行讀取出來
    list = line.strip('\n').split(' ')      #處理逐行資料:strip表示把頭尾的'\n'去掉,split表示以空格來分割行資料,然後把處理後的行資料返回到list列表中
    A[A_row:] = list[0:3]                    #把處理後的資料放到方陣A中。list[0:3]表示列表的0,1,2列資料放到矩陣A中的A_row行
    A_row+=1                                #然後方陣A的下一行接著讀

x=A[:,0:2]
y=A[:,2]
x=np.mat(x)
x=x.T
y=np.mat(y)
y=y.T


# 初始化引數
ths=np.array([1,1],dtype=float)
ths=np.mat(ths).T
m=len(y)
N=10
lr=0.001
J=[]

def mysigmoid(x):      
    h=1/(1+np.vectorize(math.exp)(-x))
    return h

# 訓練引數
for it in range(N):
    sum1= np.vectorize(math.log)(mysigmoid(-ths.T*x))*y+np.vectorize(math.log)(1-mysigmoid(-ths.T*x))*(1-y)
    temp= mysigmoid(-ths.T*x)
    sum2= x* (temp.T-y)
    J.append(-sum1[0,0]/m)
    dths= sum2/m
    ths=ths-lr*dths
    
plt.xlabel("n")#設定X軸標籤  次數
plt.ylabel("J")#設定Y軸標籤  代價函式
N1 = np.linspace(0, N, N)
plt.plot(N1, J, color="blue")
plt.show()#顯示影象


H=mysigmoid(-ths.T*x)

# -*- coding: utf-8 -*-
"""
Created on Mon Nov  5 17:16:00 2018

@author: user
"""

import numpy as np
 
 
class LogisticRegression:
    """
    Logistic Regression:邏輯迴歸
    Author:CommissarMa
    """
 
    def __init__(self, m, n, X, y, alpha, iterThreshold):
        """
        建構函式:初始化成員變數
        :param m: 記錄數量(int)
        :param n: 特徵數量(int)
        :param X: 記錄矩陣(n*m)(float)
        :param y: 類別向量(1*m)(取值範圍:0或1)
        :param alpha: 更新速率(float)
        :param iterThreshold: 梯度下降迭代停止條件(float)
        :var w: 引數向量(n*1)(float)
        :var b: 引數(float)
        """
        self.m = m
        self.n = n
        self.X = X
        self.y = y
        self.alpha = alpha
        self.iterThreshold = iterThreshold
        self.w = np.ones((n, 1))
        self.b = 0
        return
 
    def train(self):
        """
        訓練:使用資料進行訓練,使用梯度下降進行迭代使得損失值不斷下降直到小於設定的迭代停止條件
        :return:訓練完成後得到最優的w和b
        """
        JLast = -1  # 用來存放上一次迭代的損失值。用-1是因為損失值>=0
        count = 0  # 迭代次數
        while True:
            count += 1
            J = 0  # 損失值
            dw = np.zeros((self.n, 1))  # a對w的導數(n*1)
            db = 0  # a對b的導數
            Z = np.dot(self.w.T, self.X) + self.b  # Z=wT*X+b
            a = 1 / (1 + np.exp(-Z))  # Sigmoid函式
            J += -(np.dot(self.y, np.log(a).T) + np.dot(1 - self.y, np.log(1 - a).T))  # 損失函式的計算
            dz = a - self.y  # a對z的導數(m*1)
            dw += np.dot(self.X, dz.T)
            db += np.sum(dz, axis=1)
            J /= self.m  # 平均損失
            dw /= self.m
            db /= self.m
            self.w -= self.alpha * dw
            self.b -= self.alpha * db
            print("第" + str(count) + "次梯度下降的損失值J:" + str(J))
            if (np.abs(J - JLast) < self.iterThreshold and JLast > 0) or count>200:
                break
            JLast = J
        return self.w, self.b
 
    def predict(self, x):
        """
        預測:對新的記錄進行預測,給出預測的類別
        :param x:需要進行預測的一條記錄(n*1)
        :return:如果預測出的概率大於0.5就返回類別1,小於等於0.5就返回類別0
        """
        result = np.dot(self.w.T, x) + self.b
        result = 1 / (1 + np.exp(-result))
        if result > 0.5:
            return 1
        else:
            return 0
 
 
if __name__ == '__main__':
    
    A = np.zeros((100,3),dtype=float)    
    f = open('。。。')               #開啟資料檔案檔案
    lines = f.readlines()           #把全部資料檔案讀到一個列表lines中
    A_row = 0                       #表示矩陣的行,從0行開始
    for line in lines:              #把lines中的資料逐行讀取出來
        list = line.strip('\n').split(' ')      #處理逐行資料:strip表示把頭尾的'\n'去掉,split表示以空格來分割行資料,然後把處理後的行資料返回到list列表中
        A[A_row:] = list[0:3]                    #把處理後的資料放到方陣A中。list[0:3]表示列表的0,1,2列資料放到矩陣A中的A_row行
        A_row+=1
    x=A[:,:2]
    m,n=np.shape(x)
    X=x.reshape([n,m])
    y=A[:,2]
    
    alpha = 0.001  # 設定更新速率
    iterThreshold = 0.00001  # 設定迭代停止條件
    lr = LogisticRegression(m, n, X, y, alpha, iterThreshold)
    lr.train()
    print(lr.predict(np.array([np.random.rand(1), np.random.rand(1)])))
    print(help(LogisticRegression))