Python Logistic Regression

Posted by kifish on April 27, 2018

见http://blog.csdn.net/zouxy09/article/details/20319673:
getA()函数是在numpy模块中的,功能是将矩阵转化成数组

主要bug原因是:有的地方把python当matlab了,比如对matrix的索引,两者其实有一些区别的。

原文中的代码可能随机化的地方有bug。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import time
import random
#sigmoid
def sigmoid(z):
    return 1.0/(1+np.exp(-z))

#opts = {alpha:a_val,max_epoch_times:m_val}
def lr_train(X,Y,opts):
    start = time.time()
    m,n = np.shape(X)
    # X = [np.ones(m,1),X] if matlab
    X = np.concatenate((np.ones((m,1)),X),axis=1)

    n += 1
    alpha = opts['alpha']
    max_epoch_times = opts['max_epoch_times']
    w = np.ones((1,n))
    w = np.mat(w)
    for epoch in range(max_epoch_times):
        if opts['optimizer'] == 'GradientDescent':
            pred = sigmoid(X*np.transpose(w))
            err = Y-pred
            #maxmize the likelihood
            w += alpha*np.transpose(err)*X

        elif opts['optimizer'] == 'StoGradDescent':
            temp_m_range = range(m)
            temp_m_range = list(temp_m_range)
            #we can also use permutation to simplify
            #actually the result is equal to direct loop
            for iter in range(m):
                sample_idx = int(random.uniform(0,len(temp_m_range)))
                pred = sigmoid(X[sample_idx,:]*np.transpose(w))
                err = Y[sample_idx]-pred
                w += alpha*err*X[sample_idx,:]
                del_idx = 0
                for idx,val in enumerate(temp_m_range):
                    if val == sample_idx:
                        del_idx = idx
                        break
                del temp_m_range[del_idx]

        #decrease alpha with iteration
        elif opts['optimizer'] == 'SmoothStoGradDescent':
            temp_m_range = range(m)
            temp_m_range = list(temp_m_range)
            #we can also use permutation to simplify
            for iter in range(m):
                alpha = 4.0 / (1.0 + epoch + iter) + 0.01
                sample_idx = int(random.uniform(0, len(temp_m_range)))
                pred = sigmoid(X[sample_idx,:] * np.transpose(w))
                err = Y[sample_idx] - pred
                w += alpha * err * X[sample_idx,:]
                #del temp_m_range[sample_idx] this code means one sample maybe selected for some times
                """
                The code below means one sample should not be selected for over one time.
                However, the code above can bring better effects which confuses me.
                Maybe SmoothStoGradDescent is deeply affected by randomness.
                """
                del_idx = 0
                for idx,val in enumerate(temp_m_range):
                    if val == sample_idx:
                        del_idx = idx
                        break
                del temp_m_range[del_idx]
        else:
            raise NameError('Invalid optimizer type!')

        print("Training is completed. It took %ss" %(time.time()-start))
        return w

def lr_test(X,Y,w):
    m,n = np.shape(X)
    X = np.concatenate((np.ones((m, 1)), X), axis=1)
    acc_cnt = 0
    for idx in range(m):
        pred = 0
        if sigmoid(X[idx,:] * np.transpose(w)) > 0.5:
            pred = 1
        if Y[idx] == pred:
            acc_cnt += 1
        acc = float(acc_cnt)/m
    return acc

def lr_plot(X,Y,w):
    m,n = np.shape(X)
    X= np.concatenate((np.ones((m, 1)), X), axis=1)
    n += 1
    #print(Y)
    #print("n = %d"%n)
    try:
        if n ==3 :
            for i in range(m):
                #note: here numpy is different from matlab
                #Y[i] != Y[i,0] because type(Y) is np.matrix
                if int(Y[i,0]) == 0:
                    plt.plot(X[i, 1], X[i, 2], 'ob')

                elif int(Y[i,0]) == 1:
                    plt.plot(X[i, 1], X[i, 2], 'or')

    except Exception as e:
        print("The number of features in invalid")

    #X is np.mat
    min_x = min(X[:,1])[0,0]
    max_x = max(X[:,1])[0,0]
    """
    print(w)
    print((type(w)))
    """
    w_array = w.getA() # conver mat to array
    """
    print(w_array)
    print(type(w_array))
    """
    y_min_x = float(-w_array[0][0]-w_array[0][1]*min_x)/w_array[0][2]
    y_max_x = float(-w_array[0][0]-w_array[0][1]*max_x)/w_array[0][2]
    plt.plot([min_x,max_x],[y_min_x,y_max_x],'-g')
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

def load_data():
    train_x = []
    train_y = []
    f = open('/Users/k/Downloads/data.txt')
    for line in f.readlines():
        line_arr = line.strip().split()
        train_x.append([float(line_arr[0]),float(line_arr[1])])
        train_y.append([float(line_arr[2])])
        #note:train_y.append([float(line_arr[2])]) should not be train_y.append(float(line_arr[2]))

    return np.mat(train_x),np.mat(train_y)

#load data
print("loading data")
train_x,train_y = load_data()
#maybe we should use other data to test
test_x = train_x
test_y = train_y

#train
print("training")
opts = {
    'alpha' :0.01,
    'max_epoch_times' : 100,
    'optimizer' : 'SmoothStoGradDescent'
    #'optimizer' : 'StoGradDescent'
    #'optimizer' : 'GradientDescent'
}

w = lr_train(train_x,train_y,opts)
#print(w)
#fit ;actually  we should use new data
acc = lr_test(test_x,test_y,w)

print("result:")
print("accuracy : %.3f%%" % (acc*100))
# (acc*100) should not be acc*100
lr_plot(train_x,train_y,w)