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

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

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

∑imod2=0xi\sum_{i \mod 2=0}x_iimod2=0∑xi

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

将向量x\mathbf{x}x中等于1的分量累加.
∑xi=1xi\sum_{x_i=1}x_ixi=1∑xi
计算10!10!10!
∏i=110i\prod_{i=1}^{10}ii=1∏10i
表示函数f(x)=3x2+x+2f(x)=3x^2+x+2f(x)=3x2+x+2在区间[1,5][1,5][1,5]上与x轴之间的面积.
∫15(3x2+x+2)dx\int_1^5(3x^2+x+2)\mathrm{d} x∫15(3x2+x+2)dx

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

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

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

例如∫12(x2+1)dx\int_1^2(x^2+1)\mathrm d x∫12(x2+1)dx,通过手算结果为103\frac{10}{3}310

integration = 0.0
delta = 0.001
x = 1.0
while x <= 2:integration += (x ** 2 + 1) * deltax += 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)(n=3, m=1)(n=3,m=1)来验证最小二乘法.

例如

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

我们假设存在f(xi)=wxi+bf(x_i)=wx_i+bf(xi)=wxi+b，使得f(xi)≈yif(x_i) \approx y_if(xi)≈yi.
很显然, 这是一个一元线性回归.
那么我们利用均方误差定义其损失函数：
L(w,b)=∑i=1n(f(xi)−yi)2L(w, b)=\sum_{i=1}^n(f(x_i)-y_i)^2L(w,b)=i=1∑n(f(xi)−yi)2
（ps:有人在定义时乘上了一个12\frac {1}{2}21, 可与求导后的222约掉.）
确定www和bbb的方法就是使得我们的损失最小, 即：
(w,b)=arg min⁡(w,b)∑i=1n(f(xi)−yi)2(w, b)= \argmin_{(w, b)}\sum_{i=1}^n(f(x_i)-y_i)^2(w,b)=(w,b)argmini=1∑n(f(xi)−yi)2
以上的求解方法就可称为“最小二乘法”.
对于损失函数L(w,b)L(w, b)L(w,b), 它是一个凸函数, 我们可以利用求导找到它的最小值.
对www求导：
∂L(w,b)∂w=2(w∑i=1n−∑i=1n(yi−b)xi)\frac {\partial L(w, b)}{\partial w} =2(w\sum_{i=1}^n-\sum_{i=1}^n(y_i-b)x_i)∂w∂L(w,b)=2(wi=1∑n−i=1∑n(yi−b)xi)
对bbb求导:
∂L(w,b)∂b=2(nb−∑i=1n(yi−wxi))\frac {\partial L(w, b)}{\partial b}=2(nb-\sum_{i=1}^n(y_i-wx_i))∂b∂L(w,b)=2(nb−i=1∑n(yi−wxi))
令上式都为零,解得：
w=∑i=1nyi(xi−x‾)∑i=1nxi2−1n(∑i=1nxi)2w= \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}w=∑i=1nxi2−n1(∑i=1nxi)2∑i=1nyi(xi−x)
b=1n∑i=1n(yi−wxi)b= \frac{1}{n}\sum_{i=1}^n(y_i-wx_i)b=n1i=1∑n(yi−wxi)
其中x‾\overline xx为x\mathbf{x}x分量的均值.
带入表格中的数据, 可以算出：
w=1.989w =1.989w=1.989.
b=−23.702b = -23.702b=−23.702.
L(w,b)=1.866L(w, b) = 1.866L(w,b)=1.866.
f(x)=1.989x−23.702f(\mathbf{x})=1.989\mathbf{x}-23.702f(x)=1.989x−23.702
相关代码如下:

import numpy as np
# import matplotlib.pyplot as plt# 最小二乘法def lost_function(w, b, data):total_cost = 0m = len(data)# 逐点计算平方损失误差，然后求平均数for i in range(m):x = data[i, 0]y = data[i, 1]total_cost += (y - w * x - b) ** 2return total_cost / m# 定义一个求均值的函数,方便使用
def average(data):sum = 0m = len(data)for i in range(m):sum += data[i]return sum / mdef fit(data):m = len(data)x_avg = average(data[:, 0])sum1 = 0sum2 = 0sum_delta = 0for i in range(m):x = data[i, 0]y = data[i, 1]sum1 += y * (x - x_avg)sum2+= x ** 2# 根据公式计算ww = 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 / mreturn w, bdata = 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][0,1][0,1]之间, 具有概率意义.
特点2: 此“回归”非彼“回归”, 实际上它是一种分类方法.
特点3: 利用极大化似然函数方法.
特点4: 处理非线性数据较麻烦.
特点5: 很难处理数据不平衡的问题.

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

假定X=[147.4154.3160.3]\mathbf{X}=\begin{bmatrix} 1 &47.4 \\ 1 &54.3 \\1&60.3 \end{bmatrix}X=⎣⎡11147.454.360.3⎦⎤,Y=[71.582.497.3]\mathbf{Y}= \begin{bmatrix} 71.5 \\ 82.4 \\97.3\end{bmatrix}Y=⎣⎡71.582.497.3⎦⎤.

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

解得w=[−23.7021.989]\mathbf{w}=\begin{bmatrix}-23.702 \\ 1.989\end{bmatrix}w=[−23.7021.989],Y^=[70.60284.33096.268]\hat\mathbf{Y}=\begin{bmatrix} 70.602 \\84.330 \\ 96.268\end{bmatrix}Y^=⎣⎡70.60284.33096.268⎦⎤, L2\rm{L}2L2损失∣∣Xw−Y∣∣22=1.866||\mathbf{X}\mathbf{w}-\mathbf{Y}||_2^2=1.866∣∣Xw−Y∣∣22=1.866

import numpy as npX = 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]]