# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # 练习2：Logistic Regression 逻辑回归

# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# ## 查看数据
#
# Suppose that you are the administrator of a university department and
# you want to determine each applicant’s chance of admission based on their
# results on two exams. You have historical data from previous applicants
# that you can use as a training set for logistic regression. For each training
# example, you have the applicant’s scores on two exams and the admissions
# decision.
#
# 也就是数据集有三个变量，分别是 exam1和exam2的分数，以及目标 admissions 

# 换用pandas来读数据
data1 = pd.read_csv('./data/ex2data1.txt',names=['exam1','exam2','admission'])
data1.head()

# ### 绘图展示下

# +
positive = data1.loc[data1['admission']==1] # 正样本
negative = data1.loc[data1['admission']==0] # 负样本

plt.figure(figsize=(10,6))
plt.plot(positive.iloc[:,0], positive.iloc[:,1], 'k+', label='Admitted')
plt.plot(negative.iloc[:,0], negative.iloc[:,1], 'yo', label='Not admitted')
plt.xlabel('Exam 1 score')
plt.ylabel('Exam 2 score')
plt.legend()
# -

# ## 初始化变量 X ，y，参数矩阵 $\theta$

# +
X = data1.iloc[:,:-1]
y = data1.iloc[:,-1:] 

X.insert(0,'x0',1) # 添加 1 列

# 转为矩阵形式
X = np.matrix(X)
y = np.matrix(y)

# 初始化参数矩阵
theta = np.matrix(np.zeros((1,3)))
X.shape, y.shape ,theta.shape


# -

# ## Sigmod函数
# \\[g\left( z \right)=\frac{1}{1+{{e}^{-z}}}\\]

def sigmoid(z):
    return 1 / (1 + np.exp(-z))


myx = np.arange(-10,10,.1)
plt.plot(myx,sigmoid(myx))
plt.title("Sigmoid Function")
plt.grid(True)


# ## 代价函数
#
# $$J(\theta)=-\frac{1}{m}\sum\limits_{i=1}^m [y^{(i)}log(h_{\theta}(x))+(1-y^{(i)})log(1-h_{\theta}(x))]$$

def computeCost(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    cost = np.multiply(-y, np.log(sigmoid(np.dot(X , theta.T)))) - np.multiply((1-y), np.log(1-(sigmoid(np.dot(X , theta.T)))))
    return np.sum(cost)/((len(X)))


print('初始代价：', computeCost(theta,X, y) ) 


# ## 梯度下降
#
# $$ \theta_j = \theta_j - \alpha\dfrac{\partial}{\partial\theta_j}J(\theta_j)  $$
# $$\dfrac{\partial}{\partial\theta}J(\theta)=\frac{1}{m}\sum\limits_{i=1}^m (h_{\theta}(x^{(i)})-y^{(i)})x^{(i)} $$

# 用于计算一次迭代的梯度
def computeGradient(theta,X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    m = len(X)
    n = theta.shape[1]  
    temp= np.zeros(X.shape[1])
    
    error = sigmoid(np.dot(X,theta.T)) - y
    for j in range(n):
        temp[j] =  ((1/m) * np.sum(np.multiply(error, X[:,j])))    
        
    return temp


computeGradient(theta, X, y)

# ## 此处要求使用fminuc函数来最小化代价，python中与该函数相近的是 fmin_tnc

import scipy.optimize as opt
result = opt.fmin_tnc(func=computeCost, x0=theta, fprime=computeGradient, args=(X, y))
result

# ### 这是另一位大佬kaleko的解法，使用的是fmin函数

# +
from scipy import optimize

def optimizeTheta(mytheta,myX,myy,mylambda=0.):
    result = optimize.fmin(computeCost, x0=mytheta, args=(myX, myy), maxiter=400, full_output=True)
    return result[0], result[1]

theta_new, mincost = optimizeTheta(theta,X,y)
# -

theta_new, mincost

# ## 绘制决策边界 Decision Boundary
#
# 我们知道假设函数（预测函数）$h_{\theta}(z)=g(\theta^{T}x)=\dfrac{1}{1+e^{-\theta^{T}x}}$，当我们得到训练好的参数后，
#
# 可知 $h_{\theta}(z)=g(\theta_0+\theta_1x_1+\theta_2x_2)$
#
# 由于sigmoid函数： 当$z>0,g(z)>0.5$，当$z<0,g(z)<0.5$，所以决策边界就是当$h_{\theta}(z)=g(\theta_0+\theta_1x_1+\theta_2x_2)=0.5$时才发生
#
# 令 $z=\theta_0+\theta_1x_1+\theta_2x_2=0$，则 $x_2=(-\frac{1}{\theta_2})(\theta_0+\theta_1x_1)$

# 预测函数 
# 假设现在有一学生 exam1 考了45分，exam2考了85分，预测其通过率
X_test = np.array([1.,45.,85.])
sigmoid(np.dot(X_test , theta.T))

# +
# 先绘制原先的散点图
positive = data1.loc[data1['admission']==1] # 正样本
negative = data1.loc[data1['admission']==0] # 负样本

plt.figure(figsize=(10,6))
plt.plot(positive.iloc[:,0], positive.iloc[:,1], 'k+', label='Admitted')
plt.plot(negative.iloc[:,0], negative.iloc[:,1], 'yo', label='Not admitted')

# 绘制决策边界
boundary_x = X[:,1]
boundary_y = (-1./theta_new[2]) * (X[:,1] * theta_new[1] + theta_new[0])

plt.plot(boundary_x,boundary_y, 'b', label='decision boundary')
plt.xlabel('Exam 1 score')
plt.ylabel('Exam 2 score')

plt.legend()
# -

# # 加入正则化的逻辑回归

# ## 查看数据
#
# In this part of the exercise, you will implement regularized logistic regression to predict whether microchips from a fabrication plant passes quality assurance (QA). During QA, each microchip goes through various tests to ensure it is functioning correctly.
#
# Suppose you are the product manager of the factory and you have the test results for some microchips on two different tests. From these two tests,you would like to determine whether the microchips should be accepted or rejected. To help you make the decision, you have a dataset of test results on past microchips, from which you can build a logistic regression model
#
# 数据集有三个变量，前两个是微芯片的两个不同的历史测试得分，后一个是目标值，是历史测试结果（通过或不通过QA）

data2 = pd.read_csv('./data/ex2data2.txt',names=['test1','test2','accepted'])
data2.head()

# +
positive = data2.loc[data2['accepted']==1] # 正样本
negative = data2.loc[data2['accepted']==0] # 负样本

plt.figure(figsize=(10,6))
plt.plot(positive.iloc[:,0], positive.iloc[:,1], 'k+', label='Accepted')
plt.plot(negative.iloc[:,0], negative.iloc[:,1], 'yo', label='Not Accepted')
plt.xlabel('Test 1 score')
plt.ylabel('Test 2 score')
plt.legend()
# -

# ## 通过Feature mapping增加一些特征
#
# $mapFeature(x) = \begin{bmatrix}1\\ x_1 \\x_2\\x_1^2\\x_1x_2\\x_2^2\\x_1^3 \\...\\x_1x_2^5\\x_2^6\end{bmatrix}$ 
#
# 最后共28维特征

# +
x1 = data2['test1']
x2 = data2['test2']

data2.insert(3, 'x0', 1)

degree = 6
for i in range(1, degree+1):
    for j in range(0, i+1):
        colsname = ('' if (i-j)==0  else ('x_1'+str(i-j))) + ('' if (j)==0 else ('x_2'+str(j)))
        data2[colsname] = np.power(x1, i-j) * np.power(x2, j)

data2.drop('test1', axis=1, inplace=True)
data2.drop('test2', axis=1, inplace=True)

data2.head()
# -

# ## 接下来同样初始化变量 X2 ，y2，参数矩阵 $\theta_2$

# +
X_2 = data2.iloc[:,1:]
y_2 = data2.iloc[:,:1]  

# 转为矩阵形式
X_2 = np.matrix(X_2)
y_2 = np.matrix(y_2)

# 初始化参数矩阵
theta_2 = np.matrix(np.zeros((1,X_2.shape[1])))
X_2.shape, y_2.shape ,theta_2.shape


# -

# ## 正则化代价函数
#
# $$J(\theta)=\frac{1}{m}\sum\limits_{i=1}^m [-y^{(i)}log(h_{\theta}(x))-(1-y^{(i)})log(1-h_{\theta}(x))]+\frac{\lambda}{2m}\sum\limits_{j=1}^{n}\theta_j^2$$
#
# 注意：$\theta_0$不参与正则化

def computeRegCost(theta, X, y, mylambda):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    m = len(X)
    cost = np.multiply(-y, np.log(sigmoid(np.dot(X , theta.T)))) - np.multiply((1-y), np.log(1-(sigmoid(np.dot(X , theta.T)))))
    reg = (mylambda/(2*m)) * (np.sum(np.power(theta[:,1:theta.shape[1]], 2)))
    return (np.sum(cost)/((len(X))) + reg)


# +
mylambda = 1.

computeRegCost(theta_2, X_2, y_2, mylambda)


# -

# ## 梯度下降
#
# $ Repeat\ until\ convergence\{\\ \theta_0:=\theta_0-\alpha\frac{1}{m}\sum\limits_{i=1}^m [(h_{\theta}(x^{(i)})-y^{(i)})x_0^{(i)}] \\ \begin{aligned} \theta_j & :=\theta_j-\alpha\frac{1}{m}\sum\limits_{i=1}^m [(h_{\theta}(x^{(i)})-y^{(i)})x^{(i)} + \frac{\lambda}{m}\theta_j] \\& for \ j=1,2,..,n\\\ \end{aligned} $
#
#

# 计算迭代一次的梯度
def computeRegGradient(theta, X, y, mylambda):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    m = len(X)
    n = theta.shape[1]  
    temp= np.zeros(X.shape[1])
    cost = sigmoid(np.dot(X , theta.T)) - y 
    for j in range(n):
        if (j == 0):
            temp[j] =  (1/m) * np.sum(np.multiply(cost, X[:,j])) 
        else:
            temp[j] =  (1/m) * np.sum(np.multiply(cost, X[:,j]) + (mylambda/m) * theta[:,j]) 
    return temp


computeRegGradient(theta_2, X_2, y_2, mylambda)


# ## 最小化参数

# ### 这里参考的是kaleko大佬的代码，使用了minimize函数

# +
def optimizeRegularizedTheta(mytheta,myX,myy,mylambda):
    result = optimize.minimize(computeRegCost, mytheta, args=(myX, myy, mylambda),  method='Newton-CG', jac=computeRegGradient )
    return np.array([result.x]), result.fun
    
theta2_new, mincost = optimizeRegularizedTheta(theta_2,X_2,y_2,mylambda)
theta2_new, mincost
# -

# ### 这是另一种优化方式

# +
result2 = opt.fmin_tnc(func=computeRegCost, x0=theta_2, fprime=computeRegGradient, args=(X_2, y_2, mylambda))

computeRegCost(np.matrix(result2[0]), X_2, y_2, mylambda)


# -

# ## 绘制决策边界
#
# 我们知道决策边界就是当 $h_{\theta}(z)=g(\theta^{T}x)=0.5$ 也就是 $z=\theta^{T}x=0$ 时才发生，所以目的就是找出使等式成立的 $x_1，x_2$

# ### 这里使用的也是kaleko大佬的代码

def mapFeature( x1col, x2col ):
    """ 
    Function that takes in a column of n- x1's, a column of n- x2s, and builds
    a n- x 28-dim matrix of featuers as described in the homework assignment
    """
    degrees = 6
    out = np.ones( (x1col.shape[0], 1) )

    for i in range(1, degrees+1):
        for j in range(0, i+1):
            term1 = x1col ** (i-j)
            term2 = x2col ** (j)
            term  = (term1 * term2).reshape( term1.shape[0], 1 ) 
            out   = np.hstack(( out, term ))
    return out


def darw_data_and_boundary(x_plt, y_plt, theta, mylambda):
    z_plt = np.zeros((len(x_plt),len(y_plt)))
    for i in range(len(x_plt)):
        for j in range(len(y_plt)):
            myfeaturesij = mapFeature(np.array([x_plt[i]]),np.array([y_plt[j]]))
            z_plt[i][j] = np.dot(theta,myfeaturesij.T)
    z_plt = z_plt.transpose()

    plt.plot(positive.iloc[:,0], positive.iloc[:,1], 'k+', label='Accepted')
    plt.plot(negative.iloc[:,0], negative.iloc[:,1], 'yo', label='Not Accepted')
    plt.xlabel('Test 1 score')
    plt.ylabel('Test 2 score')
    plt.legend()
    u, v = np.meshgrid( x_plt, y_plt )
    mycontour = plt.contour( u, v, z_plt, [0])
    plt.title("Decision Boundary- lambda=%.2f" %mylambda)


# +
x1_min = data2['x_11'].min()
x1_max = data2['x_11'].max()
x2_min = data2['x_21'].min()
x2_max = data2['x_21'].max()
x_plt = np.linspace(x1_min,x1_max,50)
y_plt = np.linspace(x2_min,x2_max,50)

plt.figure(figsize=(10,6))
darw_data_and_boundary(x_plt, y_plt, theta2_new, mylambda)
# -

# ### 修改lambda值，观察不同情况下的决策边界

# +
plt.figure(figsize=(16,10))
plt.subplot(221)
mylambda1 = 0.
theta2_new1, mincost1 = optimizeRegularizedTheta(theta_2,X_2,y_2,mylambda1)
darw_data_and_boundary(x_plt, y_plt, theta2_new1,mylambda1)

plt.subplot(222)
mylambda2 = 1.
theta2_new2, mincost2 = optimizeRegularizedTheta(theta_2,X_2,y_2,mylambda2)
darw_data_and_boundary(x_plt, y_plt, theta2_new2,mylambda2)

plt.subplot(223)
mylambda3 = 10.
theta2_new3, mincost3 = optimizeRegularizedTheta(theta_2,X_2,y_2,mylambda3)
darw_data_and_boundary(x_plt, y_plt, theta2_new3,mylambda3)


plt.subplot(224)
mylambda4 = 100.
theta2_new4, mincost4 = optimizeRegularizedTheta(theta_2,X_2,y_2,mylambda4)
darw_data_and_boundary(x_plt, y_plt, theta2_new4,mylambda4)
# -

# ## 使用sklearn工具包

# +
from sklearn import linear_model#调用sklearn的线性回归包
model = linear_model.LogisticRegression(penalty='l2', C=1.0)
yy2 = np.array(y_2)
model.fit(X_2, yy2.ravel())

model.score(X_2, yy2) # 模型准确率
# -



StarcoderdataJupyterScriptsDedupFiltered