汪群超 Chun-Chao Wang

Dept. of Statistics, National Taipei University, Taiwan

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

Objective:

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

Prerequisite:

  1. 一般函數圖形的繪製,包括數值方式(Numerical)與符號方式(Symbolic)。
  2. 以數值與符號的方式計算積分。

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

  1. Solve f(x) = x^2-3x+2=0 for x

  2. Solve f(x) = x^4-8x^3+16x^2-2x+8=0 for x

  3. 繪製函數圖形並標示所有根的位置。


注意事項:

  1. 多項式函數是特別的函數,有專屬的解根的指令,譬如 numpy.polynomial.polynomial.polyroots

  2. 本範例只針對實數根,因此必須將虛根排除。

  3. 計算一般函數的根,並非如對多項式函數一樣能以一個指令全部找到。

  4. 下圖的根較為集中,可以再畫一張圖,聚焦在根的附近,放大函數通過 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:以數值方法計算一般函數(非多項式函數)的根

  1. Solve f(x) = x - e^{-x/2} = 0 for x

  2. 繪製函數圖形並標示所有根的位置。


注意事項:

  1. 針對非多項式函數,通常一次僅計算一個根(在設定的範圍內,如 bracket=[0.5, 1]),本範例採 sicpy.optimize.root_scalar 指令。 而 sicpy.optimize.root_scalar 的計算結果並「不單純」,到底如何找到對應的根, 讀者可以嘗試列印 print(dir(sol)) 出來看看,有哪些 attributes 或 method 可用。

  2. sicpy.optimize.root_scalar 所針對的函數 f,可以副程式 def 的方式表達或以匿名函數 f = lambda x:… 的方式。一般而言,簡單的函數以單行的匿名函數表示,較複雜的函數需要多行指令的,便以 def 方式。

  3. 下列程式碼使用 “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 的兩個根。

注意事項:

  1. 先以目視方式給予適當的初始值(對 sicpy.optimize.root)或範圍(對 sicpy.optimize.root_scalar)。

  2. 以迴圈方式透過均勻分布的隨機亂數指定初始值,以散彈槍打鳥的方式,計算出相當數量的根,最後再整理出幾其中 unique 的根(numpy.unique)。


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

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

注意事項:

  1. 這個問題可以改寫為較明顯的方式: f(x) = 0, where f(x) = F(x) - 0.95

  2. F(x) 正是標準常態分配的累積密度函數。

  3. 在計算複雜函數的根之前,可以試著畫出函數來觀察根的大略位置,以便與計算出來的根做比對。這裡採用向量函數的方式來計算帶著積分的函數。

  4. 對向量函數不熟悉的讀者,可以先試著用迴圈來畫圖。先畫出來,再來仔細端詳向量函數的精簡。

  5. 這個函數 f(x) 比較複雜些,因此採用 def 副程式的方式來包裝,並且因附加一個參數,在呼叫 sicpy.optimize.root_scalre 時,必須加入 args 參數,作為函數額外的常數資料,在此就是機率值 prob=0.95。

  6. 右下圖順便畫出 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 為分配參數,

  1. 求解 \beta(x;2, 6) = \beta(x;4, 2)。即求解 f(x) =  \beta(x;2, 6) - \beta(x;4, 2) = 0 的解。

  2. 繪製函數圖 f(x) 並標示答案的位置。

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


範例 4:以符號的方式計算函數的根(Symbolic root-finding)

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

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


注意事項:

  1. sympy.solvesympy.solveset 都是以符號的方式計算函數的根。

  2. 以符號方式計算的根,仍是符號,必須經過轉換才能當數值使用。

  3. 程式碼最下方被註解的部分,用來列印符號的結果,讀者可以打開來執行看看。

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 的最小值。

注意事項:

  1. 本範例使用 scipy.optimize.minimize_scalar 計算函數的區域最小值(Local minimum)。

  2. 指令 scipy.optimize.minimize_scalar 一次只能計算一個區域最小值,決定區域的參數是 bracket。決定的方式可以用兩個數字或三個數字,請參閱手冊的詳細解釋,譬如,本範例採 bracket = [2, 3, 5], 意思是 f(3) < f(2)f(3) < f(5)

  3. 指令 scipy.optimize.minimize_scalar 可以選擇演算法,內設的演算法是 Brent。其他演算法請詳見手冊。

  4. 下列程式碼只找出一個區域最小值,請自行找出另一個,並標註記號。

  5. 本範例使用的 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}

注意事項:

  1. 函數 \tan^{-1}() 可以用 numpy.arctan()

  2. 計算最最小值之前,最好先畫出函數圖,對最小值的位置有大概的掌握。

  3. 建議讀者嘗試使用其他演算法,譬如 res = opt.minimize_scalar(f, bounds = [0.1, 2], method = ‘bounded’), 其中演算法 methodbounded 搭配的參數 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}

注意事項:

  1. 函數的區域最大值計算被歸類於最小值的計算,換句話說,沒有為最大值準備的指令,只需將函數加上負號即可。即 \max_x f(x) \equiv \min_x -f(x).

  2. 本範例的函數較為複雜,不適合處裡單行函數設定的 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()

習題: 除第 4 題外,皆需畫出函數圖並標示出關鍵值的位置(最小值、最大值或實數根)。

  1. \displaystyle \min_x \sqrt{\frac{x^2+1}{x+1}}

  2. \displaystyle \min_{-4 \leq x \leq 3} (x+1)^5 sin(x-3)

  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.

  4. 計算 \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。畫一張折線圖,呈現每個樣本數的 \lambda 估計值。為得到在每一個樣本數下較穩定的估計值,對每個樣本數皆執行 10,000 次,最後取其平均數作為估計值,如下圖左。也因為執行了 10,000 次,因此可以得出估計值的標準差,同樣也為每個樣本數的標準差畫一張折線圖,如下圖右。
    請注意:本題雖然可以直接以紙筆推演出最後的封閉解(樣本平均值),不過為配合本章的主題,仍採計算對數概似函數最大值的方式進行。另外,也可以嘗試不取對數的做法,即 \displaystyle \max_{\lambda}  \Pi_{i=1}^N f(x_i; \lambda)

商學院  7F16
ccw@gm.ntpu.edu.tw
(02)8674-1111 
ext 66777

部落格統計

  • 107,553 點擊次數
%d bloggers like this: