scipy题目（高级编程技术 week-13）

╰+哭是因爲堅強的太久メ 2022-05-23 05:34 102阅读 0赞

# 第一题 Least squares #

![Acrobat_2018-06-02_23-51-12.png][]

关于这一道题，使用`scipy`库中的leastsq函数即可求解。

## 第一步：先生成可用的矩阵和向量 ##

在本题中，我需要先生成矩阵A  A ，向量b b，考虑到后期该问题需要有解，向量b  b 不随机生成，而是先随机生成x，然后算出b。代码如下

import scipy as sp
    import numpy as np
    import matplotlib.pyplot as plt
    # 维数
    m = 20
    n = 10
    # 生成矩阵A
    A = np.random.normal(size=(m,n), scale=15, loc=10)
    # 生成假想的该问题应有的解
    xx = np.random.normal(size=n,scale=15, loc=10)
    # 生成该解对应的向量b
    b = np.dot(A, xx)

在生成了可用的数据后，便使用`leastsq`函数解决问题，代码如下

from scipy.optimize import leastsq
    # 问题求解初始向量
    x = np.ones(10)
    # 误差函数
    def error(xx):
        return np.dot(A, xx) - b
    
    x_result,cov_x = leastsq(error, x)
    
    print(x_result, cov_x)
    # 计算误差向量的范数
    norm_residual = np.linalg.norm(np.dot(A, x_result) - b)
    print(norm_residual)

运行以上代码后，可以得到下面的解：

![Snipaste_2018-06-02_23-57-51.png][]

# 第二题 Optimization #

求函数最大值，使用scipy自带的`fmin`函数即可，这里稍微转换一下，由于只有fmin这个函数可用，而我们找的又是最大值，我们就只需要给函数加个负数，就变成了求最小值，这样子就把问题转换过来了。下面是代码：

def func(x):
        return -(np.sin(x-2) * np.sin(x-2) * np.exp(-x**2))
    minimum = sp.optimize.fmin(func, 1)
    print(minimum)

求得的结果是

![chrome_2018-06-03_00-04-36.png][]

经过绘图检验，可见的确是在x≈0.2162 x ≈ 0.2162左右这个点处取得最大值。

![chrome_2018-06-03_00-02-34.png][]

# 第三题 Pairwise distances #

![Acrobat_2018-06-03_00-07-43.png][]

在scipy文档中，找到了这样的函数可以实现计算行之间的距离的功能：`pdist`

![zeal_2018-06-03_00-29-28.png][]

因此实现的代码如下：

from scipy.spatial.distance import pdist
    m = 30
    n = 20
    X = np.random.normal(size=(n, m), scale=10, loc=10)
    Y = pdist(X, 'sqeuclidean')
    print(Y.shape)
    print(Y)

注意，这里需要说明的数，这里使用的是欧拉距离公式，因此在函数参数中有一个`sqeuclidean`的项。

文档内容可见

![zeal_2018-06-03_00-32-28.png][]

结果可见

![chrome_2018-06-03_00-30-27.png][]

经过检验，的确应该是190个数（毕竟 20 \* 19 /2 )， 其余的结果也是正确的。

[Acrobat_2018-06-02_23-51-12.png]: https://lh3.googleusercontent.com/-dDrRXVqhk3Q/WxK883E-l5I/AAAAAAAAIGw/aqsnsx414eUyVTsb4Er8ROkleqtKVnIWQCHMYCw/s0/Acrobat_2018-06-02_23-51-12.png
[Snipaste_2018-06-02_23-57-51.png]: https://lh3.googleusercontent.com/-s2hBkNNeRjc/WxK-gTRv7DI/AAAAAAAAIHA/JRLBcM8fMsQueKjikFqli8j4ofSbpyGGACHMYCw/s0/Snipaste_2018-06-02_23-57-51.png
[chrome_2018-06-03_00-04-36.png]: https://lh3.googleusercontent.com/-IvETNo_juTo/WxLAFdrZDwI/AAAAAAAAIHc/7_Py8ig2HPAkhIxYjZLLkd_o7WOzy339gCHMYCw/s0/chrome_2018-06-03_00-04-36.png
[chrome_2018-06-03_00-02-34.png]: https://lh3.googleusercontent.com/-SkeVUHyos0w/WxK_mgpVPXI/AAAAAAAAIHQ/jse6SOy5lcoIsnnyhtVsCvLig_ziCXnAwCHMYCw/s0/chrome_2018-06-03_00-02-34.png
[Acrobat_2018-06-03_00-07-43.png]: https://lh3.googleusercontent.com/-5aEeTC6bIP8/WxLA0AuD_aI/AAAAAAAAIHk/MJxRkLpbFcsJY3WMCFoq2ipq6rS4H6PiQCHMYCw/s0/Acrobat_2018-06-03_00-07-43.png
[zeal_2018-06-03_00-29-28.png]: https://lh3.googleusercontent.com/-Tsoa0YGUm1E/WxLF6UXmGWI/AAAAAAAAIHw/jhuLv28CyiYp9fFCZxAIb0kHXuhNTCVRwCHMYCw/s0/zeal_2018-06-03_00-29-28.png
[zeal_2018-06-03_00-32-28.png]: https://lh3.googleusercontent.com/-NRkKbXszA4A/WxLGnXHik_I/AAAAAAAAIIE/mxnQxzD8N7YJp2HtGXbzLBUe4IG3SGLrQCHMYCw/s0/zeal_2018-06-03_00-32-28.png
[chrome_2018-06-03_00-30-27.png]: https://lh3.googleusercontent.com/-tX2q0DHAWfw/WxLGJLypZYI/AAAAAAAAIH0/VEq0zx9SVJoUgWjUmFzC9Kz6VvEeZJAKACHMYCw/s0/chrome_2018-06-03_00-30-27.png