Python 演算法 Day 4 - 理论基础 微积分

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)  + 3

def kd(x):
    return -20*x + 100

from matplotlib import pyplot as plt
x = 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,即为原曲线的"极值"。
https://ithelp.ithome.com.tw/upload/images/20210707/20138527dR80OUeiPY.png

4-2. Critical Points 临界点

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

def v(x):
    return (x**3) - (2*x) + 100
def vd(x):
    return 3*(x**2) - 2

from matplotlib import pyplot as plt
x = 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()

https://ithelp.ithome.com.tw/upload/images/20210707/20138527GQ0sqSC38u.png

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) + 2
def vd(x):
    return (3*(x**2)) - (12*x) + 12
def v2d(x):
    return (3*(2*x)) - 12

from matplotlib import pyplot as plt
x = 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()

https://ithelp.ithome.com.tw/upload/images/20210707/20138527koUmGEEFuQ.png

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
梯度逼近示意图:
https://ithelp.ithome.com.tw/upload/images/20210707/20138527oyxVDwN1O6.png

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

import numpy as np
import matplotlib.pyplot as plt

# 目标函数(损失函数):f(x)
def func(x): return x ** 2 # np.square(x)
# 目标函数的一阶导数:f'(x)
def dfunc(x): return 2 * x

def 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 arange
t = 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。
https://ithelp.ithome.com.tw/upload/images/20210719/20138527Y1xHQBAGe4.png

情境2. 若超参数如下:

x_start = 5
epochs = 15
learning_rate = 0.01

会发现 Learning rate 太小,根本还没求到 min 就停止前进。
(可以透过增加 epochs 来达到同样效果)
https://ithelp.ithome.com.tw/upload/images/20210719/20138527DLn3WKJcrK.png

情境3. 若超参数如下:

x_start = 5
epochs = 15
learning_rate = 0.9

会发现 Learning rate 太大,逼近速度较快但有跳过最佳解的风险。
https://ithelp.ithome.com.tw/upload/images/20210719/20138527Nhiljy0lbu.png

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

.

6. Integration 积分

6-1. Basic

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

import matplotlib.pyplot as plt
import numpy as np

# Define function f
def f(x):
    return x
x = 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 间面积。
https://ithelp.ithome.com.tw/upload/images/20210710/20138527YUaHnuToHx.png

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

import scipy.integrate as integrate
# Define function f
def 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 inte
import numpy as np

i, 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。
https://ithelp.ithome.com.tw/upload/images/20210712/20138527gG94osiDaL.png
.
.
.
.
.

Homework Answer:

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

# 2i. 
import sympy as sp
x = sp.Symbol('x')
sp.solvers.solve((x+2)*(x-3)-(x-5)*(x-6), x)
>> [18/5]

# 或:
improt sympy as sp
exp1='(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=True
x = 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 np
left = 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 solve
eq = []
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}
  1. 使用 tkinter,让使用者输入方程序,自动计算答案。
from sympy.solvers import solve

def 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()

实际画面如下:
https://ithelp.ithome.com.tw/upload/images/20210630/2013852787e6dqbOfC.png

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

import streamlit as st
# 要用 cmd : streamlit run f:/LessonOnClass/Part2/20210417/00.Practice.py
import sympy as sp

x, y, z = sp.symbols('x y z')

exp = st.text_area('请输入联立方程序: ', 'x+y=5\nx-y=1') # title1
st.markdown('# 联立方程序求解') # title2

if 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 写入网页中

实际画面如下:
https://ithelp.ithome.com.tw/upload/images/20210707/20138527yuyEI2b689.png

Homework:

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

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

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

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

  4. f(x) = sin(x)E*(-0.1*(x-0.6)**2)
    *提示:用 y_prime = y.diff(x) 求微分後方程序。


<<:  Day5 参加职训(机器学习与资料分析工程师培训班),Python程序设计

>>:  Exactly how To Come To Be A Salesforce Omni Programmer Qualification?

Day29 NiFi 与其他工具的比较

这边我想特别写出这一篇的原因是当初我在学习与操作 NiFi 的过程中,我曾感到一些疑惑,会觉得感觉...

[13th][Day23] http response header(下)

接续昨天 response headers 的部分 一样是看 Julia Evans 大大的可爱的图...

TypeScript | namespace 心得纪录

这边一样是在研究 D14 介面的函式超载运用,又提到 namespace 的知识点,所以回头来补相关...

连续 30 天 玩玩看 ProtoPie - Day 15

昨天顺利让使用者讲的话能够显示在萤幕上後,来做点有用的语音控制吧。 Trigger 选择一个 Voi...

见习村27 - First non-repeating character

27 - First non-repeating character Don't say so mu...