Lesson 7: 單變量函數的根與最小值

Objective:

以數值與符號的方式計算單變量函數的根 $f(x) = 0$ 。
以數值與符號的方式計算單變量函數的最小值 $min_x f(x)$ 。
熟悉在函數圖形中標示根或極值位置。

Prerequisite:

一般函數圖形的繪製，包括數值方式（Numerical）與符號方式（Symbolic）。
以數值與符號的方式計算積分。

範例 1：以數值方法計算多項式函數的根

Solve $f(x) = x^2-3x+2=0$ for $x$
Solve $f(x) = x^4-8x^3+16x^2-2x+8=0$ for $x$
繪製函數圖形並標示所有根的位置。

注意事項：

多項式函數是特別的函數，有專屬的解根的指令，譬如 numpy.polynomial.polynomial.polyroots。
本範例只針對實數根，因此必須將虛根排除。
計算一般函數的根，並非如對多項式函數一樣能以一個指令全部找到。
下圖的根較為集中，可以再畫一張圖，聚焦在根的附近，放大函數通過 $y=0$ 的視野。

import numpy as np
import matplotlib.pyplot as plt 
import numpy.polynomial.polynomial as poly

coef = [2, -3, 1] # from low to high
r = poly.polyroots(coef)
print('The roots are {}'.format(r))

coef = [8, -2, 16, -8, 1]
r = poly.polyroots(coef)
r_x =np.real(r[np.isreal(r)]) 
x = np.linspace(-1, 5, 100)
y = poly.polyval(x, coef)

plt.style.use('ggplot')
plt.plot(x, y, color = 'g', linewidth = 2)
plt.scatter(r_x, [0, 0], s = 20, c = 'r')
plt.grid(True)
plt.title('Root finding of a polynomial function')
plt.show()

範例 2：以數值方法計算一般函數（非多項式函數）的根

Solve $f(x) = x - e^{-x/2} = 0$ for $x$ 。
繪製函數圖形並標示所有根的位置。

注意事項：

針對非多項式函數，通常一次僅計算一個根（在設定的範圍內，如 bracket=[0.5, 1]），本範例採 sicpy.optimize.root_scalar 指令。而 sicpy.optimize.root_scalar 的計算結果並「不單純」，到底如何找到對應的根, 讀者可以嘗試列印 print(dir(sol)) 出來看看，有哪些 attributes 或 method 可用。
sicpy.optimize.root_scalar 所針對的函數 f，可以副程式 def 的方式表達或以匿名函數 f = lambda x:… 的方式。一般而言，簡單的函數以單行的匿名函數表示，較複雜的函數需要多行指令的，便以 def 方式。
下列程式碼使用 “plt.text” 來標示根的位置，為了準確地將符號落在根的位置，水平（horizontal）與垂直（vertical）的微調是必要的。

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def f(x) :
    return x - np.exp(-x/2)
# f = lambda x: x - np.exp(-x/2)
# ---- First, draw the graph ------
x = np.linspace(0, 3, 1000)
plt.plot(x, f(x), color = 'g', linewidth = 1)
# ---- Second, compute the root ------
sol = opt.root_scalar(f, bracket=[0.5, 1], method='brentq')  
print('The root is at x = {:.4f}'.format(sol.root))  
# mark the root on top of function line
plt.text(sol.root, 0, 'X', color = 'r',
    horizontalalignment='center',
    verticalalignment='center')  
plt.title('$f(x) = x - e^{-x/2}$ and its root')
plt.xlabel('X'), plt.grid(True)
plt.show()

另一個計算一般函數根的指令 sicpy.optimize.root。這個指令以初始值，如下方程式碼內的 x0 = 1，做為演算法尋找根的開端。當函數有多個根時，必須藉改變初始值的設定，才能找到不同的根。當然，如果能先畫出函數圖，較能準確地給定初始值，否則必須盲目以亂數方式做多次的嘗試。

import scipy.optimize as opt

f = lambda x : x - np.exp(-x/2)
sol = opt.root(f, x0 = 1) 
# print(dir(sol))
print('The root is x = {:.4f}'.format(sol.x[0]))
print('The function value is f = {:.6f}'.format(sol.fun[0]))

練習：計算一般多根函數的根

分別以指令 sicpy.optimize.root 與 sicpy.optimize.root_scalar 計算函數 $f(x) = x^4-8x^3+16x^2-2x+8=0$ 的兩個根。

注意事項：

先以目視方式給予適當的初始值（對 sicpy.optimize.root）或範圍（對 sicpy.optimize.root_scalar）。
以迴圈方式透過均勻分布的隨機亂數指定初始值，以散彈槍打鳥的方式，計算出相當數量的根，最後再整理出幾其中 unique 的根（numpy.unique）。

範例 3：以數值方法計算較複雜函數的根

計算函數的根 $F(x) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}}e^{-t^2/2} \;dt = 0.95$

注意事項：

這個問題可以改寫為較明顯的方式： $f(x) = 0$ , where $f(x) = F(x) - 0.95$
$F(x)$ 正是標準常態分配的累積密度函數。
在計算複雜函數的根之前，可以試著畫出函數來觀察根的大略位置，以便與計算出來的根做比對。這裡採用向量函數的方式來計算帶著積分的函數。
對向量函數不熟悉的讀者，可以先試著用迴圈來畫圖。先畫出來，再來仔細端詳向量函數的精簡。
這個函數 $f(x)$ 比較複雜些，因此採用 def 副程式的方式來包裝，並且因附加一個參數，在呼叫 sicpy.optimize.root_scalre 時，必須加入 args 參數，作為函數額外的常數資料，在此就是機率值 prob=0.95。
右下圖順便畫出 $F(x)$ 的意義。

import numpy as np
import scipy.optimize as opt
from scipy.stats import norm
import matplotlib.pyplot as plt
import scipy.integrate as integral


def f(x, prob) :
    g = lambda x : np.exp(-x**2/2) / np.sqrt(2*np.pi)
    return integral.quad(g, -np.inf, x)[0] - prob

# Draw function graph
prob = 0.95
x = np.linspace(-5, 5, 100)
vec_f = np.vectorize(f)
F = vec_f(x, prob)
plt.plot(x, F, linewidth = 2)

sol = opt.root_scalar(f, args = prob, bracket = [0, 3])   # pass prob as an extra argument
print('The root is x = {:.3f}'.format(sol.root)) 

plt.hlines(0, x.min(), x.max(), color = 'g')
plt.vlines(sol.root, F.min(), F.max(), color = 'g')
plt.grid(True)
plt.title('$F(x) - {:.2f}$'.format(prob))
plt.show()

# Draw the area in pdf 
y = norm.pdf(x)
plt.plot(x, y)
x_area = np.linspace(x.min(), sol.root, 100 )
y_area = norm.pdf(x_area)
plt.fill_between(x_area, y_area, 0, color = 'b', alpha = 0.3)
plt.text(0, 0.1, prob)
plt.title('$F(x)={:.2f}, x = {:.3f}$'.format(prob, sol.root))
plt.show()

練習：以數值方法計算較複雜函數的根

令 $\beta(x;a,b)$ 代表 $\beta$ 機率密度函數（PDF)， $a, b$ 為分配參數，

求解 $\beta(x;2, 6) = \beta(x;4, 2)$ 。即求解 $f(x) = \beta(x;2, 6) - \beta(x;4, 2) = 0$ 的解。
繪製函數圖 $f(x)$ 並標示答案的位置。

請注意： $\beta$ 機率密度函數（PDF) 可以直接採用 scipy.beta.pdf(…)

範例 4：以符號的方式計算函數的根（Symbolic root-finding）

利用符號運算套件 sympy 求解下列方程式，並以符號繪圖的方式繪製函數圖形。

$f(x) = x^2-3x+2 = 0$
$f(x) = x - e^{-x/2} = 0$

注意事項：

sympy.solve 與 sympy.solveset 都是以符號的方式計算函數的根。
以符號方式計算的根，仍是符號，必須經過轉換才能當數值使用。
程式碼最下方被註解的部分，用來列印符號的結果，讀者可以打開來執行看看。

import sympy as sym
from sympy import plot

x = sym.Symbol('x')
coef = [1, -3, 2] # from high to low
f = np.polyval(coef, x)
# f = x**2 - 3*x + 2
r = sym.solve(f)
# convert r to float numbers
r_float = [float(x) for x in r]

fig1 = plot(f, xlim = [-1, 4], ylim = [-1, 8], \
    show = False, title = 'Symbolic root-finding {}'.format(r_float))
fig2 = plot(0*x, xlim = [-1, 4], line_color = 'g', show = False)
fig1.append(fig2[0])
fig1.show()

# Alternative: sympy claims that solveset will replace solve in the future
rset = sym.solveset(f, x, domain = sym.S.Reals) #sym.Interval(0, 1)
rset_float = [float(x) for x in rset] # convert to float
print('The roots are {}'.format(rset_float))
#-------------------------------------------
# sym.pprint(f)
# exp = sym.sin(x) / sym.cos(x)
# sym.pprint(exp)
# sym.pprint(sym.simplify(exp))
#--------------------------------------------

範例 5：數值最小值計算 Numerical minimization

計算函數 $f(x) = x^4 - 8x^3+16x^2-2x+8$ 的最小值。

注意事項：

本範例使用 scipy.optimize.minimize_scalar 計算函數的區域最小值（Local minimum）。
指令 scipy.optimize.minimize_scalar 一次只能計算一個區域最小值，決定區域的參數是 bracket。決定的方式可以用兩個數字或三個數字，請參閱手冊的詳細解釋，譬如，本範例採 bracket = [2, 3, 5], 意思是 $f(3) < f(2)$ 及 $f(3) < f(5)$ 。
指令 scipy.optimize.minimize_scalar 可以選擇演算法，內設的演算法是 Brent。其他演算法請詳見手冊。
下列程式碼只找出一個區域最小值，請自行找出另一個，並標註記號。
本範例使用的 scipy.optimize.minimize_scalar 只能找到一個最小值，而且是區域最小值（指定了範圍）。如果函數有絕對最大值，程式該怎麼寫，才能找到絕對最大值。

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

# draw the graph
coef = [1, -8, 16, -2, 8] # high to low
f= lambda x : np.polyval(coef, x)
x = np.linspace(-1, 5, 1000)
plt.plot(x, f(x))
plt.xlabel('X'), plt.grid(True)

# compute the local minimum
res = opt.minimize_scalar(f, bracket=[2, 3, 5])
print(res) # find out the return values
print('The function has a local minimum at x = {:.4f}'.format(res.x))
print('The corresponding function value is {:.4f}'.format(res.fun))

plt.text(res.x, res.fun, 'X', color = 'r',
    horizontalalignment='center',
    verticalalignment='center')
plt.title('The local minimum of $f(x)$')
plt.show()

練習：數值最小值計算 Numerical minimization

$\displaystyle \min_{x \in (0,3)} \tan^{-1}\frac{5}{x}+\tan^{-1}\frac{2}{3-x}$

注意事項：

函數 $\tan^{-1}()$ 可以用 numpy.arctan()。
計算最最小值之前，最好先畫出函數圖，對最小值的位置有大概的掌握。
建議讀者嘗試使用其他演算法，譬如 res = opt.minimize_scalar(f, bounds = [0.1, 2], method = ‘bounded’), 其中演算法 method 採 bounded 搭配的參數 bounds 用來設定搜尋範圍。

範例 6：數值最大值計算 Numerical maximization

計算函數的區域最大值

$\displaystyle \max_{0 \leq \theta \leq 3\pi} V(\theta)$ ,
where $V(\theta)=\frac{\pi}{3}(\frac{10\pi - 5\theta}{2\pi})^2\sqrt{25 - (\frac{10\pi - 5\theta}{2\pi})^2}$

注意事項：

函數的區域最大值計算被歸類於最小值的計算，換句話說，沒有為最大值準備的指令，只需將函數加上負號即可。即 $\max_x f(x) \equiv \min_x -f(x)$ .
本範例的函數較為複雜，不適合處裡單行函數設定的 lambda，而是採 def 函數設定方式。

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

# First, draw function graph
# define a function for a complicated funcion
def V(t) :
    tmp = ((10*np.pi - 5*t) / (2*np.pi))**2
    V = np.pi / 3 * tmp * np.sqrt(25 - tmp)
    return V

t = np.linspace(0, 3*np.pi, 1000)
plt.plot(t, V(t))

f = lambda x : -V(x)
res = opt.minimize_scalar(f, bounds = [0, 2])

plt.text(res.x, -res.fun, '+', color ='r', fontsize = 16,
    horizontalalignment='center',
    verticalalignment='center')
plt.xlabel('X'), plt.grid(True)
plt.title('The fnction maximum is {:.4f} at x = {:.4f}'.format(-res.fun, res.x))
plt.show()

習題： 每題都需要畫出目標函數的圖形並標示出關鍵值的位置（最小值、最大值或實數根）。

$\displaystyle \min_x \sqrt{\frac{x^2+1}{x+1}}$
$\displaystyle \min_{-4 \leq x \leq 3} (x+1)^5 sin(x-3)$
計算 $L(x) = 10$ 的解 $x$ , 其中
$L(x) = \int_a^x \sqrt{1 + (f'(t))^2}\; dt, \;\;$ for $f(t) = t^2/2$ and $a = 0$ .
最大概似函數估計（MLE）：計算 $\displaystyle \max_{\lambda} \ln \Pi_{i=1}^N f(x_i; \lambda)$ ,

其中 $f(x_i; \lambda)$ 代表指數分配（參數 $\lambda$ ）的概似函數，即 $f(x_i; \lambda) = \lambda e^{-\lambda x_i}$ 。令樣本數 $N= 10, 20 ,30, 50, 100, 300, 500$ ，分別生成樣本 $x_i$ （令真實 $\lambda = 2$ ，或自己設定），並採最大概似估計法(log MLE)估計 $\lambda$ 。
- 任取一組樣本，繪製目標函數 $\displaystyle \ln \Pi_{i=1}^N f(x_i; \lambda)$ ，並標示出最大值的位置。
- 畫一張折線圖，呈現每個樣本數的 $\lambda$ 估計值。為得到在每一個樣本數下較穩定的估計值，對每個樣本數皆執行 10,000 次，最後取其平均數作為估計值，如下圖左。也因為執行了 10,000 次，因此可以得出估計值的標準差，同樣也為每個樣本數的標準差畫一張折線圖，如下圖右。
  請注意：本題雖然可以直接以紙筆推演出最後的封閉解（樣本平均值），不過為配合本章的主題，仍採計算對數概似函數最大值的方式進行。另外，也可以嘗試不取對數的做法，即 $\displaystyle \max_{\lambda} \Pi_{i=1}^N f(x_i; \lambda)$