2.scipy Optimize

1import pylab as pl
2import numpy as np

1import matplotlib as mpl
2mpl.rcParams['font.sans-serif'] = ['SimHei']

 1from math import sin, cos
 2from scipy import optimize
 3
 4def f(x): #❶
 5    x0, x1, x2 = x.tolist() #❷
 6    return [
 7        5*x1+3,
 8        4*x0*x0 - 2*sin(x1*x2),
 9        x1*x2 - 1.5
10    ]
11
12# f计算方程组的误差，[1,1,1]是未知数的初始值
13result = optimize.fsolve(f, [1,1,1]) #❸
14print (result)
15print (f(result))

 1def j(x):  #❶
 2    x0, x1, x2 = x.tolist()
 3    return [[0, 5, 0],
 4            [8 * x0, -2 * x2 * cos(x1 * x2), -2 * x1 * cos(x1 * x2)],
 5            [0, x2, x1]]
 6
 7
 8result = optimize.fsolve(f, [1, 1, 1], fprime=j)  #❷
 9print(result)
10print(f(result))

 1import numpy as np
 2from scipy import optimize
 3
 4X = np.array([ 8.19,  2.72,  6.39,  8.71,  4.7 ,  2.66,  3.78])
 5Y = np.array([ 7.01,  2.78,  6.47,  6.71,  4.1 ,  4.23,  4.05])
 6
 7def residuals(p): #❶
 8    "计算以p为参数的直线和原始数据之间的误差"
 9    k, b = p
10    return Y - (k*X + b)
11
12# leastsq使得residuals()的输出数组的平方和最小，参数的初始值为[1,0]
13r = optimize.leastsq(residuals, [1, 0]) #❷
14k, b = r[0]
15print ("k =",k, "b =",b)

 1#%figonly=最小化正方形面积之和（左），误差曲面（右）
 2scale_k = 1.0
 3scale_b = 10.0
 4scale_error = 1000.0
 5
 6def S(k, b):
 7    "计算直线y=k*x+b和原始数据X、Y的误差的平方和"
 8    error = np.zeros(k.shape)
 9    for x, y in zip(X, Y):
10        error += (y - (k * x + b)) ** 2
11    return error
12
13ks, bs = np.mgrid[k - scale_k:k + scale_k:40j, b - scale_b:b + scale_b:40j]
14error = S(ks, bs) / scale_error
15
16from mpl_toolkits.mplot3d import Axes3D
17from matplotlib.patches import Rectangle
18
19fig = pl.figure(figsize=(12, 5))
20
21ax1 = pl.subplot(121)
22
23ax1.plot(X, Y, "o")
24X0 = np.linspace(2, 10, 3)
25Y0 = k*X0 + b
26ax1.plot(X0, Y0)
27
28for x, y in zip(X, Y):
29    y2 = k*x+b
30    rect = Rectangle((x,y), abs(y-y2), y2-y, facecolor="red", alpha=0.2)
31    ax1.add_patch(rect)
32
33ax1.set_aspect("equal")
34
35
36ax2 = fig.add_subplot(122, projection='3d')
37
38ax2.plot_surface(
39    ks, bs / scale_b, error, rstride=3, cstride=3, cmap="jet", alpha=0.5)
40ax2.scatter([k], [b / scale_b], [S(k, b) / scale_error], c="r", s=20)
41ax2.set_xlabel("$k$")
42ax2.set_ylabel("$b$")
43ax2.set_zlabel("$error$");

 1#%fig=带噪声的正弦波拟合
 2def func(x, p):  #❶
 3    """
 4    数据拟合所用的函数: A*sin(2*pi*k*x + theta)
 5    """
 6    A, k, theta = p
 7    return A * np.sin(2 * np.pi * k * x + theta)
 8
 9
10def residuals(p, y, x):  #❷
11    """
12    实验数据x, y和拟合函数之间的差，p为拟合需要找到的系数
13    """
14    return y - func(x, p)
15
16
17x = np.linspace(0, 2 * np.pi, 100)
18A, k, theta = 10, 0.34, np.pi / 6  # 真实数据的函数参数
19y0 = func(x, [A, k, theta])  # 真实数据
20# 加入噪声之后的实验数据
21np.random.seed(0)
22y1 = y0 + 2 * np.random.randn(len(x))  #❸
23
24p0 = [7, 0.40, 0]  # 第一次猜测的函数拟合参数
25
26# 调用leastsq进行数据拟合
27# residuals为计算误差的函数
28# p0为拟合参数的初始值
29# args为需要拟合的实验数据
30plsq = optimize.leastsq(residuals, p0, args=(y1, x))  #❹
31
32print(u"真实参数:", [A, k, theta])
33print(u"拟合参数", plsq[0])  # 实验数据拟合后的参数
34
35pl.plot(x, y1, "o", label=u"带噪声的实验数据")
36pl.plot(x, y0, label=u"真实数据")
37pl.plot(x, func(x, plsq[0]), label=u"拟合数据")
38pl.legend(loc="best")

1def func2(x, A, k, theta):
2    return A*np.sin(2*np.pi*k*x+theta)   
3
4popt, _ = optimize.curve_fit(func2, x, y1, p0=p0)
5print (popt)

1popt, _ = optimize.curve_fit(func2, x, y1, p0=[10, 1, 0])
2
3print(u"真实参数:", [A, k, theta])
4
5print(u"拟合参数", popt)

 1def target_function(x, y):
 2    return (1 - x)**2 + 100 * (y - x**2)**2
 3
 4
 5class TargetFunction(object):
 6    def __init__(self):
 7        self.f_points = []
 8        self.fprime_points = []
 9        self.fhess_points = []
10
11    def f(self, p):
12        x, y = p.tolist()
13        z = target_function(x, y)
14        self.f_points.append((x, y))
15        return z
16
17    def fprime(self, p):
18        x, y = p.tolist()
19        self.fprime_points.append((x, y))
20        dx = -2 + 2 * x - 400 * x * (y - x**2)
21        dy = 200 * y - 200 * x**2
22        return np.array([dx, dy])
23
24    def fhess(self, p):
25        x, y = p.tolist()
26        self.fhess_points.append((x, y))
27        return np.array([[2 * (600 * x**2 - 200 * y + 1), -400 * x],
28                         [-400 * x, 200]])
29
30
31def fmin_demo(method):
32    target = TargetFunction()
33    init_point = (-1, -1)
34    res = optimize.minimize(
35        target.f,
36        init_point,
37        method=method,
38        jac=target.fprime,
39        hess=target.fhess)
40    return res, [
41        np.array(points) for points in (target.f_points, target.fprime_points,
42                                        target.fhess_points)
43    ]
44
45
46methods = ("Nelder-Mead", "Powell", "CG", "BFGS", "Newton-CG", "L-BFGS-B")
47for method in methods:
48    res, (f_points, fprime_points, fhess_points) = fmin_demo(method)
49    print(
50        "{:12s}: min={:12g}, f count={:3d}, fprime count={:3d}, fhess count={:3d}"
51        .format(method, float(res["fun"]), len(f_points), len(fprime_points),
52                len(fhess_points)))

 1#%figonly=各种优化算法的搜索路径
 2def draw_fmin_demo(f_points, fprime_points, ax):
 3    xmin, xmax = -3, 3
 4    ymin, ymax = -3, 3
 5    Y, X = np.ogrid[ymin:ymax:500j,xmin:xmax:500j]
 6    Z = np.log10(target_function(X, Y))
 7    zmin, zmax = np.min(Z), np.max(Z)
 8    ax.imshow(Z, extent=(xmin,xmax,ymin,ymax), origin="bottom", aspect="auto", cmap="gray")
 9    ax.plot(f_points[:,0], f_points[:,1], lw=1)
10    ax.scatter(f_points[:,0], f_points[:,1], c=range(len(f_points)), s=50, linewidths=0)
11    if len(fprime_points):
12        ax.scatter(fprime_points[:, 0], fprime_points[:, 1], marker="x", color="w", alpha=0.5)
13    ax.set_xlim(xmin, xmax)
14    ax.set_ylim(ymin, ymax)
15
16fig, axes = pl.subplots(2, 3, figsize=(9, 6))
17methods = ("Nelder-Mead", "Powell", "CG", "BFGS", "Newton-CG", "L-BFGS-B")
18for ax, method in zip(axes.ravel(), methods):
19    res, (f_points, fprime_points, fhess_points) = fmin_demo(method)
20    draw_fmin_demo(f_points, fprime_points, ax)
21    ax.set_aspect("equal")
22    ax.set_title(method)

 1def func(x, p):
 2    A, k, theta = p
 3    return A*np.sin(2*np.pi*k*x+theta)
 4
 5def func_error(p, y, x):
 6    return np.sum((y - func(x, p))**2)   
 7
 8x = np.linspace(0, 2*np.pi, 100)
 9A, k, theta = 10, 0.34, np.pi/6 
10y0 = func(x, [A, k, theta])
11np.random.seed(0)
12y1 = y0 + 2 * np.random.randn(len(x))

1result = optimize.basinhopping(func_error, (1, 1, 1),
2                      niter = 10,
3                      minimizer_kwargs={"method":"L-BFGS-B",
4                                        "args":(y1, x)})
5print (result.x)

1#%figonly=用`basinhopping()`拟合正弦波
2pl.plot(x, y1, "o", label=u"带噪声的实验数据")
3pl.plot(x, y0, label=u"真实数据")
4pl.plot(x, func(x, result.x), label=u"拟合数据")
5pl.legend(loc="best");