smale学习之数学表达式(day3)

闵老板帖子：数学表达式魔训
文章中的代码均由python实现.

---

Q1.将向量下标为偶数的分量 (x2, x4, …) 累加, 写出相应表达式.

$$\sum _{i\mathrm{m}\mathrm{o}\mathrm{d}2=0}{x}_i$$

Q2.各出一道累加、累乘、积分表达式的习题, 并给出标准答案.

• 将向量$\mathbf{x}$中等于1的分量累加.
  $$\sum _{x_i=1}{x}_i$$
• 计算$10!$
  $$\prod _{i=1}^{10}i$$
• 表示函数$f(x)=3x^2+x+2$在区间$[1,5]$上与x轴之间的面积.
  $${\int }_1^5(3x^2+x+2)\mathrm{d}x$$

Q3.你使用过三重累加吗? 描述一下其应用.

没有, 但使用过三重循环.

Q4.给一个常用的定积分, 将手算结果与程序结果对比.

例如${\int }_1^2(x^2+1)\mathrm{d}x$,通过手算结果为$\frac{10}{3}$

integration = 0.0
delta = 0.001
x = 1.0
while x <= 2:
    integration += (x ** 2 + 1) * delta
    x += delta
print(integration)

输出为:

3.336833499999815

可以发现是很接近的.
python中也有专门计算定积分的库.

from sympy import *
x = symbols('x')
print(integrate(x**2+1, (x, 1, 2)))

输出为:

10/3

Q5.自己写一个小例子$(n=3,m=1)$来验证最小二乘法.

例如

$\mathbf{x}$	$\mathbf{y}$
47.4	71.5
54.3	82.4
60.3	97.3

我们假设存在$f(x_i)=w{x}_i+b$，使得$f(x_i)\approx y_i$.
很显然, 这是一个一元线性回归.
那么我们利用均方误差定义其损失函数：
$$L(w,b)=\sum _{i=1}^n(f(x_i)-y_i{)}^2$$
（ps:有人在定义时乘上了一个$\frac{1}{2}$, 可与求导后的$2$约掉.）
确定$w$和$b$的方法就是使得我们的损失最小, 即：
$$(w,b)=\underset{(w,b)}{\mathrm{argmin}}\sum _{i=1}^n(f(x_i)-y_i{)}^2$$
以上的求解方法就可称为“最小二乘法”.
对于损失函数$L(w,b)$, 它是一个凸函数, 我们可以利用求导找到它的最小值.
对$w$求导：
$$\frac{\partial L(w,b)}{\partial w}=2(w\sum _{i=1}^n-\sum _{i=1}^n(y_i-b)x_i)$$
对$b$求导:
$$\frac{\partial L(w,b)}{\partial b}=2(nb-\sum _{i=1}^n(y_i-w{x}_i))$$
令上式都为零,解得：
$$w=\frac{{\sum }_{i=1}^n{y}_i(x_i-\overline{x})}{{\sum }_{i=1}^n{x_i}^2-\frac{1}{n}({\sum }_{i=1}^n{x}_i{)}^2}$$
$$b=\frac{1}{n}\sum _{i=1}^n(y_i-w{x}_i)$$
其中$\overline{x}$为$\mathbf{x}$分量的均值.
带入表格中的数据, 可以算出：
$w=1.989$.
$b=-23.702$.
$L(w,b)=1.866$.
$f(\mathbf{x})=1.989\mathbf{x}-23.702$
相关代码如下:

import numpy as np
# import matplotlib.pyplot as plt


# 最小二乘法

def lost_function(w, b, data):
    total_cost = 0
    m = len(data)

    # 逐点计算平方损失误差，然后求平均数
    for i in range(m):
        x = data[i, 0]
        y = data[i, 1]
        total_cost += (y - w * x - b) ** 2

    return total_cost / m


# 定义一个求均值的函数,方便使用
def average(data):
    sum = 0
    m = len(data)
    for i in range(m):
        sum += data[i]
    return sum / m


def fit(data):
    m = len(data)
    x_avg = average(data[:, 0])

    sum1 = 0
    sum2 = 0
    sum_delta = 0

    for i in range(m):
        x = data[i, 0]
        y = data[i, 1]
        sum1 += y * (x - x_avg)
        sum2+= x ** 2
    # 根据公式计算w
    w = sum1 / (sum2 - m * (x_avg ** 2))

    for i in range(m):
        x = data[i, 0]
        y = data[i, 1]
        sum_delta += (y - w * x)
    b = sum_delta / m

    return w, b


data = np.genfromtxt('data1.csv', delimiter=',')

# 获取x,y两列数据
x = data[:, 0]
y = data[:, 1]
w, b = fit(data)

print("w is: ", w)
print("b is: ", b)

cost = lost_function(w, b, data)

print("cost is: ", cost)

输出为:

w is:  1.9895608351331777
b is:  -23.702951763858255
cost is:  1.8658615310775

data1.csv如下:

Q6.逻辑回归推导, 并说出这个方法的特点.

• 推导过程省略.
• 特点1: 使用Sigmoid函数将值压缩到$[0,1]$之间, 具有概率意义.
• 特点2: 此“回归”非彼“回归”, 实际上它是一种分类方法.
• 特点3: 利用极大化似然函数方法.
• 特点4: 处理非线性数据较麻烦.
• 特点5: 很难处理数据不平衡的问题.

---

补充Q5, $\mathbf{X}$应该是一个$3\times 2$的矩阵.

假定$\mathbf{X}=\left[\begin{array}{cc}1& 47.4\\ 1& 54.3\\ 1& 60.3\end{array}\right]$,$\mathbf{Y}=\left[\begin{array}{c}71.5\\ 82.4\\ 97.3\end{array}\right]$.

由$\mathbf{w}=({\mathbf{X}}^{\mathrm{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{Y}$

解得$\mathbf{w}=\left[\begin{array}{c}-23.702\\ 1.989\end{array}\right]$,$\hat{\mathbf{Y}}=\left[\begin{array}{c}70.602\\ 84.330\\ 96.268\end{array}\right]$, $\mathrm{L}2$损失$||\mathbf{X}\mathbf{w}-\mathbf{Y}|{|}_2^2=1.866$

import numpy as np

X = np.array([[1, 47.4], [1, 54.3], [1, 60.3]])
Y = np.array([71.5, 82.4, 97.3])
Y = np.matrix(Y)
Y = Y.T  # 1行3列转变为3行1列
print("矩阵Y为:")
print(Y)
print('--------')  # 分割符
X = np.matrix(X)
print("矩阵X为:")
print(X)
print('--------')
X_T = np.transpose(X)  # X的转置矩阵
print(X_T)
print('--------')
temp = np.matmul(X_T, X)  # 矩阵乘法
print(temp)
print('--------')
temp = temp.I  # 求得其逆矩阵
print(temp)
print('--------')
temp = np.matmul(temp, X_T)
print(temp)
print('--------')
w = np.matmul(temp, Y)  # 得到w的矩阵
print("矩阵w为:")
print(w)
print('--------')
Y_hat = np.matmul(X, w)
print("预测Y为:")
print(Y_hat)

输出为:

矩阵Y为:
[[71.5]
 [82.4]
 [97.3]]
--------
矩阵X为:
[[ 1.  47.4]
 [ 1.  54.3]
 [ 1.  60.3]]
--------
[[ 1.   1.   1. ]
 [47.4 54.3 60.3]]
--------
[[3.00000e+00 1.62000e+02]
 [1.62000e+02 8.83134e+03]]
--------
[[ 3.53225342e+01 -6.47948164e-01]
 [-6.47948164e-01  1.19990401e-02]]
--------
[[ 4.60979122e+00  1.38948884e-01 -3.74874010e+00]
 [-7.91936645e-02  3.59971202e-03  7.55939525e-02]]
--------
矩阵w为:
[[-23.70295176]
 [  1.98956084]]
--------
预测Y为:
[[70.60223182]
 [84.33020158]
 [96.26756659]]