Chap.I 理论基础

Part 2：微积分

4. Critical Points and Optimization 临界点与优化

Optimization：找到最佳解的方式。

4-1. Function Direction at a Point 某一点的函数方向

若 f'(x) = 0，代表其 x 解点位为 f(x) 极值 (二次函数为 max or min)。

def k(x):    return -10*(x**2) + (100*x)  + 3def kd(x):    return -20*x + 100from matplotlib import pyplot as pltx = list(range(0, 11))y = [k(i) for i in x]yd = [kd(i) for i in x]     # y 微分后的线段# 作图plt.xlabel('x (time in seconds)')plt.ylabel('k(x) (height in feet)')plt.xticks(range(0,15, 1))plt.yticks(range(-200, 500, 20))plt.grid()plt.plot(x, y, 'g', x, yd, 'purple')    # 原有曲线/微分后线plt.show()

其中要注意的地方是，当 yd=0，即为原曲线的"极值"。

4-2. Critical Points 临界点

有时 f'(x) = 0 时，并不一定会有 max or min。

def v(x):    return (x**3) - (2*x) + 100def vd(x):    return 3*(x**2) - 2from matplotlib import pyplot as pltx = list(range(-10, 11))y = [v(i) for i in x]yd = [vd(i) for i in x]# 作图plt.xlabel('x')plt.ylabel('v(x)')plt.xticks(range(-10,15, 1))plt.yticks(range(-1000, 2000, 100))plt.grid()plt.plot(x, y, 'g', x, yd, 'purple')plt.scatter(0, (2/3)**0.5, c='r') # 极值点plt.show()

4-3. Second Order Derivatives

二阶导数用于判断一阶导数是 max, min or inflexion
a. 若 f''(x) > 0，表示 f'(x)=0 过后会上升，此函数有最小值。
b. 若 f''(x) < 0，表示 f'(x)=0 过后会下降，此函数有最大值。
c. 若 f''(x) = 0，表示此函数仅有转折点。

def v(x):    return (x**3) - (6*(x**2)) + (12*x) + 2def vd(x):    return (3*(x**2)) - (12*x) + 12def v2d(x):    return (3*(2*x)) - 12from matplotlib import pyplot as pltx = list(range(-5, 11))y = [v(i) for i in x]yd = [vd(i) for i in x]y2d = [v2d(i) for i in x]# 作图plt.xlabel('x')plt.ylabel('v(x)')plt.xticks(range(-10,15, 1))plt.yticks(range(-2000, 2000, 50))plt.grid()plt.plot(x, y, 'g', x, yd, 'purple', x, y2d, 'b')plt.scatter(2, v(2), c='r')plt.show()

5. Partial Derivatives 偏导数（偏微分）

5-1. Computing a Gradient 计算梯度

自变数 > 1，就叫做梯度非斜率
f(x, y) = x^2 + y^2
f(x, y)/dx = 2x
f(x, y)/dy = 2y
梯度逼近示意图：

5-2. Using the gradient 利用梯度找最小值

import numpy as npimport matplotlib.pyplot as plt# 目标函数（损失函数）：f(x)def func(x): return x ** 2 # np.square(x)# 目标函数的一阶导数：f'(x)def dfunc(x): return 2 * xdef GD(x_start, df, epochs, lr):    """ 梯度下降法。给定起始点与目标函数的一阶导函数，求在epochs次反覆运算中x的更新值        其原理为让起始 x = x - dfunc(x)*lr，直到 dfunc(x) = 0，x 就不会再变，求得极值。        :param x_start: x的起始点        :param df: 目标函数的一阶导函数        :param epochs: 反覆运算週期        :param lr: 学习率        :return: x在每次反覆运算后的位置（包括起始点），长度为epochs+1     """        # 此迴圈用作记录所有起始点用    xs = np.zeros(epochs+1)     # array[0, 0, ...0] 有1001个    x = x_start                 # 起始点    xs[0] = x    for i in range(epochs):     # i = 0~999        x += -df(x) * lr        # x = x - dfunc(x) * learning_rate        xs[i+1] = x    return xs

情境1. 超参数如下：

# 超参数（Hyperparameters）：在模型优化（训练）前，可调整的参数。x_start = 2                 # 起始点 (可由任意点开始)epochs = 1000               # 执行週期数 (程式走多少迴圈停止)learning_rate = 0.01        # 学习率 (每次前进的步伐多大)# 梯度下降法x = GD(x_start, dfunc, epochs, learning_rate)print(x)>> [-5. -4. -3.2 -2.56 -2.05 -1.64 -1.31 -1.05 -0.84 -0.67 -0.54 -0.43 -0.34 -0.27 -0.22 -0.18]from numpy import aranget = arange(-3, 3, 0.01)plt.plot(t, func(t), c='b')     # 画出所求原函数plt.plot(x, func(x), c='r', label='lr={}'.format(learning_rate)) # 画出红连线plt.scatter(x, func(x), c='r')  # 画出红点plt.legend()plt.show()

会发现其慢慢逼近 min。

情境2. 若超参数如下：

x_start = 5epochs = 15learning_rate = 0.01

会发现 Learning rate 太小，根本还没求到 min 就停止前进。
(可以透过增加 epochs 来达到同样效果)

情境3. 若超参数如下：

x_start = 5epochs = 15learning_rate = 0.9

会发现 Learning rate 太大，逼近速度较快但有跳过最佳解的风险。

设计迴归方式时，超参数的选定会显得至关重要。

.

6. Integration 积分

6-1. Basic

积分概念为「函数区域面积总和」。

import matplotlib.pyplot as pltimport numpy as np# Define function fdef f(x):    return xx = range(0, 11)y = [f(a) for a in x]# 作图plt.xlabel('x')plt.ylabel('f(x)')plt.grid()plt.plot(x, y, color='purple')p = np.arange(0, 2, 1/20)plt.fill_between(p, f(p), color='orange')plt.show()

下图橘色区块，即是"函数 f(x)=x 在 x=0~2 间面积。

6-2. 利用 scipy 计算定积分

import scipy.integrate as integrate# Define function fdef f(x):    return x + 1# 计算 x=a→b 定积分 intergrate.quad(func, a, b)i, e = integrate.quad(lambda x: f(x), 0, 2)i, e    # i 为定积分结果，e 为计算误差值>>  4.0 ,  4.440892098500626e-14

6-3. Infinite limits 无穷极限

利用 numpy 的 inf 计算无穷大，如下：

import scipy.integrate as inteimport numpy as npi, e = inte.quad(lambda x: np.e**(-5*x), 0, np.inf)print(i)>>  0.20000000000000007

6-4. Normal distribution 常态分配

A. 常态分配机率积分

nor = lambda x: np.exp(-x**2 / 2.0)/np.sqrt(2.0 * np.pi)i, e = inte.quad(nor, -np.inf, np.inf)print('Total Probability:', i)>>  Total Probability: 0.9999999999999998

B. 一倍标準差

i, e = inte.quad(nor, -1, 1)print('1-sigma:', i)>>  1-sigma: 0.682689492137086

C. 1.96 倍标準差（信赖水準 5%）

i, e = inte.quad(nor, -1.96, 1.96)print('1.96-sigma:', i)>>  2-sigma: 0.9500042097035591 

D. 三倍标準差（可达 99.7%）

i, e = inte.quad(nor, -3, 3)print('3-sigma:', i)>> 3-sigma: 0.9973002039367399# 作图x = np.linspace(-4, 4, 1000)y = nor(x)plt.xlabel('x')plt.ylabel('Normal Distribution')plt.grid()plt.plot(x, y, color='purple')sig1 = np.linspace(-1, 1, 500)   # 68.3%sig2 = np.linspace(-2, 2, 500)   # 95.4% (假设检定常用 2 或 1.96-sigma)sig3 = np.linspace(-3, 3, 500)   # 99.7%plt.fill_between(sig3, nor(sig3), color='0.5')plt.fill_between(sig2, nor(sig2), color='c')plt.fill_between(sig1, nor(sig1), color='orange')plt.show()

下图橘色为 1-sigma，青色为 2-sigma，灰色为 3-sigma。

.
.
.
.
.

Homework Answer:

Ex.1 (2i)(https://reurl.cc/DgbDqe)：
(x+2)(x-3) = (x-5)(x-6)

# 2i. import sympy as spx = sp.Symbol('x')sp.solvers.solve((x+2)*(x-3)-(x-5)*(x-6), x)>> [18/5]# 或：improt sympy as spexp1='(x+2)*(x-3), (x-5)*(x-6)'sp.solvers.solve(sp.core.sympify(f"Eq({exp1})"))>> [18/5]

Ex.2 (3c)(https://reurl.cc/DgbDqe)：
|10-2x| = 6

# 若要计算绝对值，须加上 real=Truex = sp.Symbol('x', real=True) # real=True 表 x 是实数sp.solvers.solve(abs(10-2*x)-6, x))>> [2 8]

Ex.3 (3c)(https://reurl.cc/EnbDQk)：
4x + 3y = 4
2x + 2y - 2z = 0
5x + 3y + z = -2

# 1. numpy 解法import numpy as npleft = str(input("请输入等号左测导数 (以','逗号隔开): "))left = left.split(',')left = list(map(int, left))right = str(input("请输入等号右测导数 (以','逗号隔开): "))right = right.split(',')right = list(map(int, right))a_shape = int(len(left)**0.5)  # 得到矩阵形状a = np.array(left).reshape(a_shape, a_shape) # 重整矩阵为正方形b = np.array(right)np.linalg.solve(a, b)>> 请输入等号左测导数 (以','逗号隔开): 4,3,0,2,2,-2,5,3,1>> 请输入等号右测导数 (以','逗号隔开): 4,0,-2>> [-11.  16.   5.]# 2. sympy 解法from sympy.solvers import solveeq = []while True:    inp = input('请输入方程式: ')    if inp =='':        break    eq.append(inp)if len(eq)>0:    for i in range(len(eq)):        arr = eq[i].split('=')        if len(arr) ==2:            eq[i] = arr[0] + '-(' + arr[1] + ')'    solve(eq)>> 请输入方程式: 4x + 3y = 4>> 请输入方程式: 2x + 2y - 2z = 0>> 请输入方程式: 5x + 3y + z = -2>> 请输入方程式: >> {x: -11, y: 16, z: 5}

from sympy.solvers import solvedef run():    eq_clean = []    exp = text.get('1.0', 'end')    # 把输入区块的东西丢进 exp    eq = exp.split('\n')            # 用 '\n' 拆分    if len(eq) > 0:                 # 若输入的东西不是空白的        for i in range(len(eq)):    # i = 输入几个多项式            if eq[i] == '':                continue            arr = eq[i].split('=')  # arr = 每个多项式用 '=' 拆分            if len(arr) == 2:       # 若只有两项 (也就是只有等号左右侧)                eq[i] = arr[0] + '-(' + arr[1] + ')'            eq_clean.append(eq[i])  # 把 eq[0] 加入 eq_clean        ans.set(solve(eq_clean))    # 把 eq_clean 丢进 sympy 求解，放入变数 ans

接着我们设计弹出视窗 (tkinter)：

import tkinter as tk# 1. 宣告一个根视窗root = tk.Tk()# 2. 建立显示字元区块# text= : 显示字元# height=1 : 高度为1# font= : 字体 大小 粗体tk.Label(root, text='请输入方程式: (2x 须输入 2*x)', height=1, font='Helvetic 18 bold').pack(fill=tk.X)# 3. 建立可输入区块text = tk.Text(root, height=5, width=30, font='Helvetic 18 bold')text.pack()# 4. 建立按纽触发 run() 函数# commend= : 事件处理，通常会接一个函数tk.Button(root, text='求解', font='Helvetic 36 bold', command=run).pack(fill=tk.X)# 5. 建立显示答案区块# textvariable= : 需求显示的字元ans = tk.StringVar()tk.Label(root, text='', height=1, font='Helvetic 18 bold', textvariable=ans).pack(fill=tk.X)# 6. 监听事件root.mainloop()

实际画面如下：

当然，可以用网页设计输入画面，会更简洁：

import streamlit as st# 要用 cmd : streamlit run f:/LessonOnClass/Part2/20210417/00.Practice.pyimport sympy as spx, y, z = sp.symbols('x y z')exp = st.text_area('请输入联立方程式: ', 'x+y=5\nx-y=1') # title1st.markdown('# 联立方程式求解') # title2if st.button('求解'):   # 若按下按钮则执行以下    eq_c = []    eq = exp.split('\n')    if len(eq)>0:        for i in range (len(eq)):            if eq[i] == '':                continue            arr = eq[i].split('=')            if len(arr) ==2:                eq[i] = arr[0] + '-(' + arr[1] + ')'            eq_c.append(eq[i])    ans = sp.solvers.solve(eq_c)    st.write(f'结果: {ans}')   # 将 ans 写入网页中

实际画面如下：

Homework:

利用梯度下降法，求得目标函数 f(x) 的区域内极值。

f(x) = x^3 - 2x + 100

f(x) = -5x^2 + 3x + 6

f(x) = 2x^4 - 3x^2 + 2x - 20

f(x) = sin(x)E*(-0.1*(x-0.6)**2)
*提示：用 y_prime = y.diff(x) 求微分后方程式。