NumPy-初学者指南中文第三版-二-

NumPy 初学者指南中文第三版（二）

原文：NumPy: Beginner’s Guide - Third Edition

协议：CC BY-NC-SA 4.0

四、为您带来便利的便利函数

如我们所见，NumPy 具有大量函数。 这些函数中的许多函数只是为了方便起见，知道这些函数将大大提高您的生产率。 这包括选择数组某些部分（例如，基于布尔条件）或处理多项式的函数。 本章提供了一个计算相关性示例，使您可以使用 NumPy 进行数据分析。

在本章中，我们将涵盖以下主题：

• 数据选择与提取
• 简单的数据分析
• 收益相关的示例
• 多项式
• 线性代数函数

在第 3 章，“熟悉常用函数”中，我们有一个数据文件可以使用。 在本章中，情况有所改善-我们现在有两个数据文件。 让我们使用 NumPy 探索数据。

相关

您是否注意到某些公司的股价会紧随其后，通常是同一行业的竞争对手？ 理论上的解释是，由于这两家公司属于同一类型的业务，因此它们面临着相同的挑战，需要相同的材料和资源，并争夺相同类型的客户。

您可能想到了许多可能的对，但是您需要检查一下真实的关系。 一种方法是查看两种股票的股票收益的相关性和因果关系。 高相关性意味着某种关系。 但是，这并不是因果关系的证明，尤其是如果您没有使用足够的数据。

实战时间 – 交易相关货币对

在本节中，我们将使用两个样本数据集，其中包含日末价格数据。 第一家公司是必和必拓（BHP），该公司活跃于石油，金属和钻石的开采。 第二个是淡水河谷（VALE），这也是一家金属和采矿公司。 因此，活动有一些重叠，尽管不是 100%。 要评估相关偶对，请按照下列步骤操作：

1. 首先，从本章示例代码目录中的 CSV 文件加载数据，特别是两种证券的收盘价，并计算收益。 如果您不记得该怎么做，请参阅第 3 章，“熟悉常用函数”中的示例。

2. 协方差告诉我们两个变量如何一起变化；无非就是相关性：

   

   使用cov()函数从返回值计算协方差矩阵（并非严格如此，但这可以让我们演示一些矩阵运算）：

   covariance = np.cov(bhp_returns, vale_returns)
   print("Covariance", covariance)
   

   协方差矩阵如下：

   Covariance [[ 0.00028179  0.00019766]
    [ 0.00019766  0.00030123]]
   
   

3. 使用diagonal()方法查看对角线上的值：

   print("Covariance diagonal", covariance.diagonal())
   

   协方差矩阵的对角线值如下：

   Covariance diagonal [ 0.00028179  0.00030123]
   
   

   请注意，对角线上的值彼此不相等。 这与相关矩阵不同。

4. 使用trace()方法计算轨迹，即对角线值的总和：

   print("Covariance trace", covariance.trace())
   

   协方差矩阵的跟踪值如下：

   Covariance trace 0.00058302354992
   
   

5. 将两个向量的相关性定义为协方差，除以向量各自标准偏差的乘积。 向量a和b的等式如下：

   print(covariance/ (bhp_returns.std() * vale_returns.std()))
   

   相关矩阵如下：

   [[ 1.00173366  0.70264666]
   [ 0.70264666  1.0708476 ]]
   
   

6. 我们将用相关系数来衡量我们偶对的相关性。 相关系数取介于 -1 和 1 之间的值。 根据定义，一组值与自身的相关性为 1。 这将是理想值； 但是，我们也会对较低的值感到满意。 使用corrcoef()函数计算相关系数（或更准确地说，相关矩阵）：

   print("Correlation coefficient", np.corrcoef(bhp_returns, vale_returns))
   

   系数如下：

   [[ 1\.          0.67841747]
   [ 0.67841747  1\.        ]]
   
   

   对角线上的值仅是BHP和VALE与它们自身的相关性，因此等于 1。在任何可能性下，都不会进行任何实际计算。 由于相关性是对称的，因此其他两个值彼此相等，这意味着BHP与VALE的相关性等于VALE与BHP的相关性。 似乎这里的相关性不是那么强。

7. 另一个要点是，正在考虑的两只股票是否同步。 如果两只股票的差额是与均值之差的两个标准差，则认为它们不同步。

   如果它们不同步，我们可以发起交易，希望它们最终能够再次恢复同步。 计算两种证券的收盘价之间的差异，以检查同步：

   difference = bhp - vale
   

   检查最后的价格差异是否不同步； 请参阅以下代码：

   avg = np.mean(difference)
   dev = np.std(difference)
   print("Out of sync", np.abs(difference[-1] – avg) > 2 * dev)
   

   不幸的是，我们还不能交易：

   Out of sync False
   
   

8. 绘图需要matplotlib；这将在第 9 章"matplotlib 绘图”中讨论。 可以按以下方式进行绘制：

   t = np.arange(len(bhp_returns))
   plt.plot(t, bhp_returns, lw=1, label='BHP returns')
   plt.plot(t, vale_returns, '--', lw=2, label='VALE returns')
   plt.title('Correlating arrays')
   
   plt.xlabel('Days')
   plt.ylabel('Returns')
   plt.grid()
   plt.legend(loc='best')
   plt.show()
   

   结果图如下所示：

   

刚刚发生了什么？

我们分析了BHP和VALE收盘价的关系。 确切地说，我们计算了他们的股票收益的相关性。 我们通过 corrcoef()函数实现了这一目标。 此外，我们看到了如何计算可以从中得出相关性的协方差矩阵。 另外，我们演示了diagonal()和trace()方法，它们分别为我们提供对角线值和矩阵迹线。 有关源代码，请参见本书代码包中的correlation.py文件：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

bhp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)

bhp_returns = np.diff(bhp) / bhp[ : -1]

vale = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)

vale_returns = np.diff(vale) / vale[ : -1]

covariance = np.cov(bhp_returns, vale_returns)
print("Covariance", covariance)

print("Covariance diagonal", covariance.diagonal())
print("Covariance trace", covariance.trace())

print(covariance/ (bhp_returns.std() * vale_returns.std()))

print("Correlation coefficient", np.corrcoef(bhp_returns, vale_returns))

difference = bhp - vale
avg = np.mean(difference)
dev = np.std(difference)

print("Out of sync", np.abs(difference[-1] - avg) > 2 * dev)

t = np.arange(len(bhp_returns))
plt.plot(t, bhp_returns, lw=1, label='BHP returns')
plt.plot(t, vale_returns, '--', lw=2, label='VALE returns')
plt.title('Correlating arrays')
plt.xlabel('Days')
plt.ylabel('Returns')
plt.grid()
plt.legend(loc='best')
plt.show()

小测验 - 计算协方差

Q1. 哪个函数返回两个数组的协方差？

1. covariance
2. covar
3. cov
4. cvar

多项式

您喜欢微积分吗？ 好吧，我喜欢它！ 微积分学中的一种思想是泰勒展开，即代表无穷级数的可微函数（请参见这里和这里）。

注意

泰勒级数的定义如下所示：

在此定义中，是在点a上计算的函数f的n阶导数。

实际上，这意味着我们可以使用高阶多项式来估计任何可微的，因此是连续的函数。 然后，我们假设较高学位的条款可以忽略不计。

实战时间 – 拟合多项式

NumPy polyfit()函数拟合多项式的一组数据点，即使基础函数不是连续的：

1. 继续使用BHP和VALE的价格数据，查看其收盘价之间的差异，并将其拟合为三次方的多项式：

   bhp=np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
   vale=np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
   t = np.arange(len(bhp))
   poly = np.polyfit(t, bhp - vale, 3)
   print("Polynomial fit", poly)
   

   多项式拟合（在此示例中，选择了三次多项式）如下：

   Polynomial fit [  1.11655581e-03  -5.28581762e-02   5.80684638e-01   5.79791202e+01]
   
   

2. 您看到的数字是多项式的系数。 用polyval()函数和我们从拟合中得到的多项式对象外插到下一个值：

   print("Next value", np.polyval(poly, t[-1] + 1))
   

   我们预测的下一个值将是：

   Next value 57.9743076081
   
   

3. 理想情况下，BHP和VALE的收盘价之间的差异应尽可能小。 在极端情况下，它有时可能为零。 使用roots()函数找出多项式拟合何时达到零：

   print( "Roots", np.roots(poly))
   

   多项式的根如下：

   Roots [ 35.48624287+30.62717062j  35.48624287-30.62717062j -23.63210575 +0.j        ]
   
   

4. 在微积分课上，您可能学到的另一件事是找到“极值”，这些极值可能是最大值或最小值。 从微积分中记住，这些是我们函数的导数为零的点。 用polyder()函数来微分多项式拟合：

   der = np.polyder(poly)
   print("Derivative", der)
   

   导数多项式的系数如下：

   Derivative [ 0.00334967 -0.10571635  0.58068464]
   
   

5. 获取导数的根：

   print("Extremas", np.roots(der))
   

   我们得到的极值如下：

   Extremas [ 24.47820054   7.08205278]
   
   

   让我们使用polyval()函数仔细检查并计算拟合值：

   vals = np.polyval(poly, t)
   

6. 现在，使用argmax()和argmin()函数找到最大值和最小值：

   vals = np.polyval(poly, t)
   print(np.argmax(vals))
   print(np.argmin(vals))
   

   这为我们提供了以下屏幕快照中所示的预期结果。 好的，结果并不完全相同，但是，如果我们退回到步骤 1，我们可以看到t是通过 arange()函数定义的：

   7
   24
   
   

   绘制数据并对其拟合 ，以得到以下曲线：

   

   显然，平滑线是拟合的，而锯齿线是基础的数据。 但是由于不太适合，您可能需要尝试更高阶的多项式。

刚刚发生了什么？

我们使用 polyfit()函数将数据拟合为多项式。 我们了解了用于计算多项式值的polyval()函数，用于返回多项式的根的roots()函数以及用于返回多项式导数的polyder()函数（ 参见polynomials.py）：

from __future__ import print_function
import numpy as np
import sys
import matplotlib.pyplot as plt

bhp=np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
vale=np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)

t = np.arange(len(bhp))
poly = np.polyfit(t, bhp - vale, 3)
print("Polynomial fit", poly)

print("Next value", np.polyval(poly, t[-1] + 1))

print("Roots", np.roots(poly))

der = np.polyder(poly)
print("Derivative", der)

print("Extremas", np.roots(der))
vals = np.polyval(poly, t)
print(np.argmax(vals))
print(np.argmin(vals))

plt.plot(t, bhp - vale, label='BHP - VALE')
plt.plot(t, vals, '-—', label='Fit')
plt.title('Polynomial fit')
plt.xlabel('Days')
plt.ylabel('Difference ($)')
plt.grid()
plt.legend()
plt.show()

勇往直前 – 改进拟合

您可以做很多事情来改进拟合。 例如，尝试使用其他幂，因为在本节中选择了三次多项式。 考虑在拟合之前对数据进行平滑处理。 平滑数据的一种方法是移动平均值。 您可以在第 3 章，“熟悉常用函数”中找到简单和 EMA 计算的示例。

余量

交易量是在投资中非常重要的变量； 它表明价格走势有多大。 平衡交易量指标是最简单的股票价格指标之一。 它基于当日和前几日的收盘价以及当日的交易量。 对于每一天，如果今天的收盘价高于昨天的收盘价，那么余额表上的交易量的值等于今天的交易量。 另一方面，如果今天的收盘价低于昨天的收盘价，那么资产负债表上交易量指标的值就是资产负债表上的交易量与今天的交易量之差。 但是，如果收盘价没有变化，那么余额表上的交易量的值为零。

实战时间 – 平衡交易量

换句话说，我们需要乘以收盘价和交易量的符号。 在本节中，我们研究解决此问题的两种方法：一种使用 NumPy sign()函数，另一种使用 NumPy piecewise()函数。

1. 将BHP数据加载到close和volume数组中：

   c, v=np.loadtxt('BHP.csv', delimiter=',', usecols=(6, 7), unpack=True)
   

   计算绝对值的变化。 用diff()函数计算收盘价的变化。 diff()函数计算两个连续数组元素之间的差，并返回包含这些差的数组：

   change = np.diff(c)
   print("Change", change)
   

   收盘价变化如下：

   Change [ 1.92 -1.08 -1.26  0.63 -1.54 -0.28  0.25 -0.6   2.15  0.69 -1.33  1.16
    1.59 -0.26 -1.29 -0.13 -2.12 -3.91  1.28 -0.57 -2.07 -2.07  2.5   1.18
   -0.88  1.31  1.24 -0.59]
   
   

2. NumPy 的sign()函数返回数组中每个元素的符号。 -1 表示负数，1 表示正数，否则返回 0。 将sign()函数应用于change数组：

   signs = np.sign(change)
   print("Signs", signs)
   

   更改数组的符号如下：

   Signs [ 1\. -1\. -1\.  1\. -1\. -1\.  1\. -1\.  1\.  1\. -1\.  1\.  1\. -1\. -1\. -1\. -1\. -1.
   -1\. -1\. -1\.  1\.  1\.  1\. -1\.  1\.  1\. -1.]
   
   

   另外，我们可以用piecewise()函数来计算 。 顾名思义，piecewise()函数逐段求值函数 。 使用适当的返回值和条件调用该函数：

   pieces = np.piecewise(change, [change < 0, change > 0], [-1, 1])
   print("Pieces", pieces)
   

   这些标志再次显示如下：

   Pieces [ 1\. -1\. -1\.  1\. -1\. -1\.  1\. -1\.  1\.  1\. -1\.  1\.  1\. -1\. -1\. -1\. -1\. -1.
   -1\. -1\. -1\.  1\.  1\.  1\. -1\.  1\.  1\. -1.]
   
   

   检查结果是否相同：

   print("Arrays equal?", np.array_equal(signs, pieces))
   

   结果如下：

   Arrays equal? True
   
   

3. 余额的大小取决于前一个收盘价的变化，因此我们无法在样本的第一天进行计算：

   print("On balance volume", v[1:] * signs)
   

   余额余额如下：

   [ 2620800\. -2461300\. -3270900\.  2650200\. -4667300\. -5359800\.  7768400.
    -4799100\.  3448300\.  4719800\. -3898900\.  3727700\.  3379400\. -2463900.
    -3590900\. -3805000\. -3271700\. -5507800\.  2996800\. -3434800\. -5008300.
    -7809799\.  3947100\.  3809700\.  3098200\. -3500200\.  4285600\.  3918800.
    -3632200.]
   
   

刚刚发生了什么？

我们计算了取决于收盘价变化的余额数量。 使用 NumPy sign()和piecewise()函数，我们遍历了两种不同的方法来确定更改的符号（请参见obv.py），如下所示：

from __future__ import print_function
import numpy as np

c, v=np.loadtxt('BHP.csv', delimiter=',', usecols=(6, 7), unpack=True)

change = np.diff(c)
print("Change", change)

signs = np.sign(change)
print("Signs", signs)

pieces = np.piecewise(change, [change < 0, change > 0], [-1, 1])
print("Pieces", pieces)

print("Arrays equal?", np.array_equal(signs, pieces))

print("On balance volume", v[1:] * signs)

模拟

通常，您会想先尝试一下。 玩耍，试验，但最好不要炸东西或变脏！ NumPy 非常适合进行实验。 我们将使用 NumPy 模拟交易日，而不会实际亏损。 许多人喜欢逢低买入，换句话说，等待股票价格下跌之后再购买。 一个变种是等待价格下跌一小部分，例如比当天的开盘价低 0.1%。

实战时间 – 使用vectorize()避免循环

vectorize()函数是 ，这是另一个可以减少程序循环次数的技巧 。 请按照以下步骤计算一个交易日的利润：

1. 首先，加载数据：

   o, h, l, c = np.loadtxt('BHP.csv', delimiter=',', usecols=(3, 4, 5, 6), unpack=True)
   

2. vectorize()函数与 Python map()函数的 NumPy 等效。 调用vectorize()函数，并将其作为参数作为calc_profit()函数：

   func = np.vectorize(calc_profit)
   

3. 现在，我们可以应用func()，就好像它是一个函数一样。 将我们获得的的func()函数结果应用于价格数组：

   profits = func(o, h, l, c)
   

4. calc_profit()函数非常简单。 首先，我们尝试以较低的开盘价购买。 如果这超出每日范围，那么很明显，我们的尝试失败，没有获利，或者我们蒙受了损失，因此将返回 0。否则，我们以收盘价卖出利润仅仅是买入价和收盘价之间的差。 实际上，查看相对利润实际上更有趣：

   def calc_profit(open, high, low, close):
      #buy just below the open
      buy = open * 0.999
      # daily range
      if low <  buy < high:
         return (close - buy)/buy
      else:
         return 0
   
   print("Profits", profits)
   

5. 假设有两天的利润为零，既没有净收益也没有亏损。 选择交易日并计算平均值：

   real_trades = profits[profits != 0]
   print("Number of trades", len(real_trades), round(100.0 * len(real_trades)/len(c), 2), "%")
   print("Average profit/loss %", round(np.mean(real_trades) * 100, 2))
   

   交易摘要如下所示：

   Number of trades 28 93.33 %
   Average profit/loss % 0.02
   
   

6. 作为乐观主义者，我们对赢得大于零的交易感兴趣。 选择获胜交易的天数并计算平均值：

   winning_trades = profits[profits > 0]
   print("Number of winning trades", len(winning_trades), round(100.0 * len(winning_trades)/len(c), 2), "%")
   print("Average profit %", round(np.mean(winning_trades) * 100, 2))
   

   获胜行业统计如下：

   Number of winning trades 16 53.33 %
   Average profit % 0.72
   
   

7. 另外，作为悲观主义者，我们对小于零的亏损交易感兴趣。 选择亏损交易的天数并计算平均值：

   losing_trades = profits[profits < 0]
   print("Number of losing trades", len(losing_trades), round(100.0 * len(losing_trades)/len(c), 2), "%")
   print("Average loss %", round(np.mean(losing_trades) * 100, 2))
   

   亏损交易统计如下：

   Number of losing trades 12 40.0 %
   Average loss % -0.92
   
   

刚刚发生了什么？

我们对函数进行向量化，这是避免使用循环的另一种方法。 我们用一个函数模拟了一个交易日，该函数返回了每天交易的相对利润。 我们打印了亏损交易和获胜交易的统计摘要（请参见simulation.py）：

from __future__ import print_function
import numpy as np

o, h, l, c = np.loadtxt('BHP.csv', delimiter=',', usecols=(3, 4, 5, 6), unpack=True)

def calc_profit(open, high, low, close):
   #buy just below the open
   buy = open * 0.999

   # daily range
   if low <  buy < high:
      return (close - buy)/buy
   else:
      return 0

func = np.vectorize(calc_profit)
profits = func(o, h, l, c)
print("Profits", profits)

real_trades = profits[profits != 0]
print("Number of trades", len(real_trades), round(100.0 * len(real_trades)/len(c), 2), "%")
print("Average profit/loss %", round(np.mean(real_trades) * 100, 2))

winning_trades = profits[profits > 0]
print("Number of winning trades", len(winning_trades), round(100.0 * len(winning_trades)/len(c), 2), "%")
print("Average profit %", round(np.mean(winning_trades) * 100, 2))

losing_trades = profits[profits < 0]
print("Number of losing trades", len(losing_trades), round(100.0 * len(losing_trades)/len(c), 2), "%")
print("Average loss %", round(np.mean(losing_trades) * 100, 2))

勇往直前 – 分析连续的获胜和失败

尽管平均利润为正，但了解我们是否必须承受连续亏损也很重要。 如果是这种情况，我们可能只剩下很少甚至没有资本，那么平均利润就无关紧要。

找出是否有这样的损失。 如果需要，您还可以找出是否有长时间的连胜纪录。

平滑

嘈杂的数据很难处理，因此我们经常需要进行一些平滑处理。 除了计算移动平均值外，我们还可以使用 NumPy 函数之一来平滑数据。

hanning()函数是由加权余弦形成的窗口函数 ）：

在上式中， N对应于窗口的大小。 在后面的章节中，我们将介绍其他窗口函数。

实战时间 – 使用hanning()函数进行平滑处理

我们将使用 hanning()函数来平滑股票收益数组，如以下步骤所示：

1. 调用hanning()函数来计算特定长度窗口的权重（在本示例中为 8），如下所示：

   N = 8
   weights = np.hanning(N)
   print("Weights", weights)
   

   权重如下：

   Weights [ 0\.          0.1882551   0.61126047  0.95048443  0.95048443  0.61126047
    0.1882551   0\.        ]
   
   

2. 使用具有标准化权重的convolve()计算BHP和VALE报价的股票收益：

   bhp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
   bhp_returns = np.diff(bhp) / bhp[ : -1]
   smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns)[N-1:-N+1]
   
   vale = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
   vale_returns = np.diff(vale) / vale[ : -1]
   smooth_vale = np.convolve(weights/weights.sum(), vale_returns)[N-1:-N+1]
   

3. 使用以下代码使用matplotlib进行绘图：

   t = np.arange(N - 1, len(bhp_returns))
   plt.plot(t, bhp_returns[N-1:], lw=1.0)
   plt.plot(t, smooth_bhp, lw=2.0)
   plt.plot(t, vale_returns[N-1:], lw=1.0)
   plt.plot(t, smooth_vale, lw=2.0)
   plt.show()
   

   该图表如下所示：

   

   上图的细线是股票收益，粗线是平滑的结果。 如您所见，这些线交叉了几次。 这些点可能很重要，因为趋势可能在那里更改了 。 或者至少，BHP与VALE的关系可能已更改。 这些拐点可能经常发生，因此我们可能希望展望未来。

4. 将平滑步骤的结果拟合为多项式，如下所示：

   K = 8
   t = np.arange(N - 1, len(bhp_returns))
   poly_bhp = np.polyfit(t, smooth_bhp, K)
   poly_vale = np.polyfit(t, smooth_vale, K)
   

5. 接下来，我们需要评估在上一步中找到的多项式彼此相等的情况。 归结为减去多项式并找到所得多项式的根。 使用polysub()减去多项式：

   poly_sub = np.polysub(poly_bhp, poly_vale)
   xpoints = np.roots(poly_sub)
   print("Intersection points", xpoints)
   

   这些点如下所示：

   Intersection points [ 27.73321597+0.j          27.51284094+0.j          24.32064343+0.j
    18.86423973+0.j          12.43797190+1.73218179j  12.43797190-1.73218179j
    6.34613053+0.62519463j   6.34613053-0.62519463j]
   
   

6. 我们得到的数字很复杂，这对我们不利（除非存在假想时间）。 使用isreal()函数检查哪些数字是实数：

   reals = np.isreal(xpoints)
   print("Real number?", reals)
   

   结果如下：

   Real number? [ True  True  True  True False False False False]
   
   

   一些数字是实数，因此请使用select()函数选择它们。 select()函数根据条件列表从选项列表中获取元素来形成数组：

   xpoints = np.select([reals], [xpoints])
   xpoints = xpoints.real
   print("Real intersection points", xpoints)
   

   实际交点如下：

   Real intersection points [ 27.73321597  27.51284094  24.32064343  18.86423973   0\.           0\.   0\.  0.]
   
   

7. 我们设法得到一些零。 trim_zeros()函数从一维数组中去除前导零和尾随零。 使用trim_zeros（）函数消除零：

   print("Sans 0s", np.trim_zeros(xpoints))
   

   零消失了，输出如下所示：

   Sans 0s [ 27.73321597  27.51284094  24.32064343  18.86423973]
   
   

刚刚发生了什么？

我们将hanning()函数应用于包含股票收益的数组。 我们用 polysub()函数减去了两个多项式。 然后，我们使用isreal()函数检查实数 ，然后使用select()函数选择实数。 最后，我们使用trim_zeros()函数从数组中剥离了零（请参见smoothing.py）：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

N = 8

weights = np.hanning(N)
print("Weights", weights)

bhp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp) / bhp[ : -1]
smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns)[N-1:-N+1]

vale = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
vale_returns = np.diff(vale) / vale[ : -1]
smooth_vale = np.convolve(weights/weights.sum(), vale_returns)[N-1:-N+1]

K = 8
t = np.arange(N - 1, len(bhp_returns))
poly_bhp = np.polyfit(t, smooth_bhp, K)
poly_vale = np.polyfit(t, smooth_vale, K)

poly_sub = np.polysub(poly_bhp, poly_vale)
xpoints = np.roots(poly_sub)
print("Intersection points", xpoints)

reals = np.isreal(xpoints)
print("Real number?", reals)

xpoints = np.select([reals], [xpoints])
xpoints = xpoints.real
print("Real intersection points", xpoints)

print("Sans 0s", np.trim_zeros(xpoints))

plt.plot(t, bhp_returns[N-1:], lw=1.0, label='BHP returns')
plt.plot(t, smooth_bhp, lw=2.0, label='BHP smoothed')

plt.plot(t, vale_returns[N-1:], '--', lw=1.0, label='VALE returns')
plt.plot(t, smooth_vale, '-.', lw=2.0, label='VALE smoothed')
plt.title('Smoothing')
plt.xlabel('Days')
plt.ylabel('Returns')
plt.grid()
plt.legend(loc='best')
plt.show()

勇往直前 – 平滑变化

试用其他平滑函数 -- hamming()，blackman()，bartlett()和kaiser()。 它们的工作方式几乎与hanning()函数相同。

初始化

到目前为止，在本书中，我们遇到了一些用于初始化数组的便捷函数。 full()和full_like()函数最近被添加到 NumPy，以使初始化更加容易。

以下简短的 Python 会话显示了这两个函数的（缩写）文档：

$ python
>>> import numpy as np
>>> help(np.full)
Return a new array of given shape and type, filled with `fill_value`.
>>> help(np.full_like)

返回形状和类型与给定数组相同的完整数组。

实战时间 – 使用full()和full_like()函数创建值初始化的数组

让我们演示full()和full_like()函数的工作方式。 如果您还不在 Python shell 中，请输入以下 ：

$ python
>>> import numpy as np

1. 使用充满数字 42 的full()函数创建一个二分之一的数组，如下所示：

   >>> np.full((1, 2), 42)
   array([[ 42.,  42.]])
   
   

   从输出可以推断出，数组元素是浮点数，这是 NumPy 数组的默认数据类型。 指定整数数据类型，如下所示：

   >>> np.full((1, 2), 42, dtype=np.int)
   array([[42, 42]])
   
   

2. full_like()函数查看输入数组的元数据，并使用该信息创建一个填充有指定值的新数组。 例如，在使用linspace()函数创建数组之后，将其用作full_like()函数的模板：

   >>> a = np.linspace(0, 1, 5)
   >>> a
   array([ 0\.  ,  0.25,  0.5 ,  0.75,  1\.  ])
   >>> np.full_like(a, 42)
   array([ 42.,  42.,  42.,  42.,  42.])
   
   

   同样，我们有一个充满42的数组。 要将数据类型更改为整数，请键入以下内容：

   >>> np.full_like(a, 42, dtype=np.int)
   array([42, 42, 42, 42, 42])
   
   

刚刚发生了什么？

我们使用full()和full_like()函数创建了数组。 full()函数用数字42填充数组。 full_like()函数使用输入数组的元数据来创建新数组。 这两个函数都可以指定数据类型。

总结

我们使用corrcoef()函数计算了两只股票的股票收益率的相关性。 另外，我们演示了diagonal() 和trace()函数，它们可以为我们提供矩阵的对角线和迹线。

我们使用polyfit()函数将数据拟合为多项式。 我们了解了用于计算多项式值的polyval()函数，用于返回多项式根的roots()函数以及用于返回多项式导数的polyder()函数。

我们看到full()函数用数字填充数组，full_like()函数使用输入数组的元数据创建一个新数组。 这两个函数都可以指定数据类型。

希望您提高了工作效率，因此我们可以在下一章继续使用矩阵和通用函数（ufuncs）。

五、使用矩阵和ufunc

本章介绍矩阵和通用函数（ufuncs）。 矩阵在数学上是众所周知的，在 NumPy 中也具有表示。 通用函数适用于数组，逐元素或标量。 ufuncs 期望一组标量作为输入，并产生一组标量作为输出。 通用函数通常可以映射到它们的数学对应物，例如加，减，除，乘等。 我们还将介绍三角函数，按位和比较通用函数。

在本章中，我们将介绍以下主题：

• 矩阵创建
• 矩阵运算
• 基本函数
• 三角函数
• 按位函数
• 比较函数

矩阵

NumPy 中的矩阵是ndarray的子类。 我们可以使用特殊的字符串格式创建矩阵。 就像在数学中一样，它们是二维的。 正如您期望的那样，矩阵乘法不同于正常的 NumPy 乘法。 幂运算符也是如此。 我们可以使用mat()，matrix()和bmat()函数创建矩阵。

实战时间 – 创建矩阵

如果输入已经是矩阵或ndarray，则mat()函数不会复制。 调用此函数等效于调用matrix(data, copy=False)。 我们还将演示转置和求逆矩阵。

1. 行用分号分隔，值用空格分隔。 使用以下字符串调用mat()函数以创建矩阵：

   A = np.mat('1 2 3; 4 5 6; 7 8 9')
   print("Creation from string", A)
   

   矩阵输出应为以下矩阵：

   Creation from string [[1 2 3]
    [4 5 6]
    [7 8 9]]
   
   

2. 如下将矩阵转换为具有T属性的矩阵：

   print("transpose A", A.T)
   

   以下是转置矩阵：

   transpose A [[1 4 7]
    [2 5 8]
    [3 6 9]]
   
   

3. 可以使用I属性将矩阵反转，如下所示：

   print("Inverse A", A.I)
   

   逆矩阵打印如下（请注意，这是O(n<sup class="calibre54">3</sup>)操作，这意味着需要平均三次时间）：

   Inverse A [[ -4.50359963e+15   9.00719925e+15  -4.50359963e+15]
    [  9.00719925e+15  -1.80143985e+16   9.00719925e+15]
    [ -4.50359963e+15   9.00719925e+15  -4.50359963e+15]]
   
   

4. 不使用字符串创建矩阵，而是使用数组：

   print("Creation from array", np.mat(np.arange(9).reshape(3, 3)))
   

   新创建的数组如下所示：

   Creation from array [[0 1 2]
    [3 4 5]
    [6 7 8]]
   
   

刚刚发生了什么？

我们使用mat()函数创建了矩阵。 我们将使用T属性将矩阵转置，并使用I属性将其反转（请参见matrixcreation.py）：

from __future__ import print_function
import numpy as np

A = np.mat('1 2 3; 4 5 6; 7 8 9')
print("Creation from string", A)
print("transpose A", A.T)
print("Inverse A", A.I)
print("Check Inverse", A * A.I)

print("Creation from array", np.mat(np.arange(9).reshape(3, 3)))

从其他矩阵创建矩阵

有时，我们想由其他较小的矩阵创建矩阵。 我们可以通过bmat()函数来实现。 b在这里代表块矩阵。

实战时间 – 从其他矩阵创建矩阵

我们将从两个较小的矩阵创建一个矩阵，如下所示：

1. 首先，创建一个2×2单位矩阵：

   A = np.eye(2)
   print("A", A)
   

   单位矩阵如下所示：

   A [[ 1\.  0.]
    [ 0\.  1.]]
   
   

2. 创建另一个类似于A的矩阵，并将其乘以 2：

   B = 2 * A
   print("B", B)
   

   第二个矩阵如下：

   B [[ 2\.  0.]
    [ 0\.  2.]]
   
   

3. 从字符串创建复合矩阵。 字符串使用与mat()函数相同的格式-使用矩阵而不是数字：

   print("Compound matrix\n", np.bmat("A B; A B"))
   

   复合矩阵如下所示：

   Compound matrix
   [[ 1\.  0\.  2\.  0.]
    [ 0\.  1\.  0\.  2.]
    [ 1\.  0\.  2\.  0.]
    [ 0\.  1\.  0\.  2.]]
   
   

刚刚发生了什么？

我们使用bmat()函数从两个较小的矩阵创建了一个块矩阵。 我们给该函数一个字符串，其中包含矩阵名称，而不是数字（请参见bmatcreation.py）：

from __future__ import print_function
import numpy as np

A = np.eye(2)
print("A", A)
B = 2 * A
print("B", B)
print("Compound matrix\n", np.bmat("A B; A B"))

小测验 – 使用字符串定义矩阵

Q1. mat()和bmat()函数接受的字符串中的行分隔符是什么？

1. 分号
2. 句号
3. 逗号
4. 空格

通用函数

通用函数（ufuncs）期望一组标量作为输入并产生一组标量作为输出。 它们实际上是 Python 对象，它们封装了函数的行为。 通常，我们可以将ufunc映射到它们的数学对应项，例如加，减，除，乘等。 通常，通用函数由于其特殊的优化以及在本机级别上运行而更快。

实战时间 – 创建通用函数

我们可以使用 NumPy 和frompyfunc()函数从 Python 函数创建ufunc，如下所示：

1. 定义一个 Python 函数，该函数回答关于宇宙，存在和其他问题的最终问题（来自《银河系漫游指南》 ，道格拉斯·亚当，Pan Books，如果您还没有阅读，可以安全地忽略这一点！）：

   def ultimate_answer(a):
   

   到目前为止，没有什么特别的。 我们给函数命名为ultimate_answer()并定义了一个参数a。

2. 使用zeros_like()函数，创建一个形状与a相同的所有零组成的结果：

   result = np.zeros_like(a)
   

3. 现在，将初始化数组的元素设置为答案42，然后返回结果。 完整的函数应显示在下面的代码片段中。 flat属性使我们可以访问平面迭代器，该迭代器允许我们设置数组的值。

   def ultimate_answer(a):
      result = np.zeros_like(a)
      result.flat = 42
      return result
   

4. 用frompyfunc（）创建一个ufunc；指定1作为输入参数的数量，然后指定1作为输出参数的数量：

   ufunc = np.frompyfunc(ultimate_answer, 1, 1)
   print("The answer", ufunc(np.arange(4)))
   

   一维数组的结果如下所示：

   The answer [42 42 42 42]
   
   

   使用以下代码对二维数组执行相同的操作：

   print("The answer", ufunc(np.arange(4).reshape(2, 2)))
   

   二维数组的输出如下所示：

   The answer [[42 42]
    [42 42]]
   
   

刚刚发生了什么？

我们定义了一个 Python 函数。 在此函数中，我们使用zeros_like()函数，根据输入参数的形状将数组的元素初始化为零。 然后，使用ndarray的flat属性，将数组元素设置为最终答案42（请参见answer42.py）：

from __future__ import print_function
import numpy as np

def ultimate_answer(a):
   result = np.zeros_like(a)
   result.flat = 42
   return result

ufunc = np.frompyfunc(ultimate_answer, 1, 1)
print("The answer", ufunc(np.arange(4)))

print("The answer", ufunc(np.arange(4).reshape(2, 2)))

通用函数的方法

函数如何具有方法？ 如前所述，通用函数不是函数，而是表示函数的 Python 对象。 通用函数具有五种重要方法，如下：

1. ufunc.reduce(a[, axis, dtype, out, keepdims])
2. ufunc.accumulate(array[, axis, dtype, out])
3. ufunc.reduceat(a, indices[, axis, dtype, out])
4. ufunc.outer(A, B)
5. ufunc.at(a, indices[, b])])])

实战时间 – 将ufunc方法应用于add函数

让我们在add()函数上调用前四个方法：

1. 通用函数沿指定元素的连续轴递归地减少输入数组。 对于add()函数，约简的结果类似于计算数组的总和。 调用reduce()方法：

   a = np.arange(9)
   print("Reduce", np.add.reduce(a))
   

   精简数组应如下所示：

   Reduce 36
   
   

2. accumulate()方法也递归地遍历输入数组。 但是，与reduce()方法相反，它将中间结果存储在一个数组中并返回它。 如果使用add()函数，则结果等同于调用cumsum()函数。 在add()函数上调用accumulate()方法：

   print("Accumulate", np.add.accumulate(a))
   

   累积的数组如下：

   Accumulate [ 0  1  3  6 10 15 21 28 36]
   
   

3. reduceat()方法的解释有点复杂，因此让我们对其进行调用并逐步了解其算法。 reduceat()方法需要输入数组和索引列表作为参数：

   print("Reduceat", np.add.reduceat(a, [0, 5, 2, 7]))
   

   结果如下所示：

   Reduceat [10  5 20 15]
   
   

   第一步与索引0和5有关。 此步骤可减少索引0和5之间的数组元素的运算：

   print("Reduceat step I", np.add.reduce(a[0:5]))
   

   步骤 1 的输出如下：

   Reduceat step I 10
   
   

   第二步涉及索引5和2。 由于2小于5，因此返回索引为5的数组元素：

   print("Reduceat step II", a[5])
   

   第二步产生以下输出：

   Reduceat step II 5
   
   

   第三步涉及索引2和7。 此步骤可减少索引2和7之间的数组元素的运算：

   print("Reduceat step III", np.add.reduce(a[2:7]))
   

   第三步的结果如下所示：

   Reduceat step III 20
   
   

   第四步涉及索引7。 此步骤导致从索引7到数组末尾的数组元素减少操作：

   print("Reduceat step IV", np.add.reduce(a[7:]))
   

   第四步结果如下所示：

   Reduceat step IV 15
   
   

4. outer()方法返回一个具有等级的数组，该等级是其两个输入数组的等级的总和。 该方法适用于所有可能的输入数组元素对。 在add()函数上调用outer()方法：

   print("Outer", np.add.outer(np.arange(3), a))
   

   外部总和的输出结果如下：

   Outer [[ 0  1  2  3  4  5  6  7  8]
    [ 1  2  3  4  5  6  7  8  9]
    [ 2  3  4  5  6  7  8  9 10]]
   
   

刚刚发生了什么？

我们将通用函数的前四种方法reduce()，accumulate()，reduceat()和outer()应用于add()函数（请参见ufuncmethods.py）：

from __future__ import print_function
import numpy as np

a = np.arange(9)

print("Reduce", np.add.reduce(a))
print("Accumulate", np.add.accumulate(a))
print("Reduceat", np.add.reduceat(a, [0, 5, 2, 7]))
print("Reduceat step I", np.add.reduce(a[0:5]))
print("Reduceat step II", a[5])
print("Reduceat step III", np.add.reduce(a[2:7]))
print("Reduceat step IV", np.add.reduce(a[7:]))
print("Outer", np.add.outer(np.arange(3), a))

算术函数

通用算术运算符+，-和*分别隐式链接到通用函数的加，减和乘。 这意味着在 NumPy 数组上使用这些运算符之一时，将调用相应的通用函数。 除法涉及一个稍微复杂的过程。 与数组分割有关的三个通用函数是divide()，true_divide()和floor_division()。 两个运算符对应于除法：/和//。

实战时间 – 分割数组

让我们来看一下数组分割的作用：

1. divide()函数将截断整数除法和普通浮点除法：

   a = np.array([2, 6, 5])
   b = np.array([1, 2, 3])
   print("Divide", np.divide(a, b), np.divide(b, a))
   

   divide()函数的结果如下所示：

   Divide [2 3 1] [0 0 0]
   
   

   如您所见，截断发生了。

2. 函数true_divide()更接近于除法的数学定义。 整数除法返回浮点结果，并且不会发生截断：

   print("True Divide", np.true_divide(a, b), np.true_divide(b, a))
   

   true_divide()函数的结果如下：

   True Divide [ 2\.          3\.          1.66666667] [ 0.5         0.33333333  0.6       ]
   
   

3. floor_divide()函数总是返回整数结果。 等效于调用divide()函数之后调用floor()函数。 floor()函数丢弃浮点数的小数部分并返回一个整数：

   print("Floor Divide", np.floor_divide(a, b), np.floor_divide(b, a))
   c = 3.14 * b
   print("Floor Divide 2", np.floor_divide(c, b), np.floor_divide(b, c))
   

   floor_divide()函数调用导致：

   Floor Divide [2 3 1] [0 0 0]
   Floor Divide 2 [ 3\.  3\.  3.] [ 0\.  0\.  0.]
   
   

4. 默认情况下，/运算符等效于调用divide()函数：

   from __future__ import division
   

   但是，如果在 Python 程序的开头找到此行，则将调用true_divide()函数。 因此，此代码将如下所示：

   print("/ operator", a/b, b/a)
   

   The result is shown as follows:

   / operator [ 2\.          3\.          1.66666667] [ 0.5         0.33333333  0.6       ]
   
   

5. //运算符等效于调用floor_divide()函数。 例如，查看以下代码片段：

   print("// operator", a//b, b//a)
   print("// operator 2", c//b, b//c)
   

   //运算符结果如下所示：

   // operator [2 3 1] [0 0 0]
   // operator 2 [ 3\.  3\.  3.] [ 0\.  0\.  0.]
   
   

刚刚发生了什么？

divide()函数会执行截断整数除法和常规浮点除法。 true_divide()函数始终返回浮点结果而没有任何截断。 floor_divide()函数始终返回整数结果； 结果与通过连续调用divide()和floor()函数（请参见dividing.py）获得的相同：

from __future__ import print_function
from __future__ import division
import numpy as np

a = np.array([2, 6, 5])
b = np.array([1, 2, 3])

print("Divide", np.divide(a, b), np.divide(b, a))
print("True Divide", np.true_divide(a, b), np.true_divide(b, a))
print("Floor Divide", np.floor_divide(a, b), np.floor_divide(b, a))
c = 3.14 * b
print("Floor Divide 2", np.floor_divide(c, b), np.floor_divide(b, c))
print("/ operator", a/b, b/a)
print("// operator", a//b, b//a)
print("// operator 2", c//b, b//c)

勇往直前 – 尝试__future__.division

通过实验确认导入__future__.division的影响。

模运算

我们可以使用 NumPy mod()，remainder()和fmod()函数来计算模或余数。 同样，我们可以使用%运算符。 这些函数之间的主要区别在于它们如何处理负数。 该组中的奇数是fmod()函数。

实战时间 – 计算模数

让我们调用前面提到的函数：

1. remainder()函数以元素为单位返回两个数组的其余部分。 如果第二个数字为 0，则返回 0：

   a = np.arange(-4, 4)
   print("Remainder", np.remainder(a, 2))
   

   remainder()函数的结果如下所示：

   Remainder [0 1 0 1 0 1 0 1]
   
   

2. mod()函数的功能与remainder()函数的功能完全相同：

   print("Mod", np.mod(a, 2))
   

   mod()函数的结果如下所示：

   Mod [0 1 0 1 0 1 0 1]
   
   

3. %运算符只是remainder()函数的简写：

   print("% operator", a % 2)
   

   %运算符的结果如下所示：

   % operator [0 1 0 1 0 1 0 1]
   
   

4. fmod()函数处理负数的方式与mod()和%的处理方式不同。 余数的符号是被除数的符号，除数的符号对结果没有影响：

   print("Fmod", np.fmod(a, 2))
   

   fmod()结果打印如下：

   Fmod [ 0 -1  0 -1  0  1  0  1]
   
   

刚刚发生了什么？

我们向 NumPy 演示了mod()，remainder()和fmod()函数，它们计算模或余数（请参见modulo.py）：

from __future__ import print_function
import numpy as np

a = np.arange(-4, 4)

print("Remainder", np.remainder(a, 2))
print("Mod", np.mod(a, 2))
print("% operator", a % 2)
print("Fmod", np.fmod(a, 2))

斐波那契数

斐波那契数 基于递归关系：

用 NumPy 代码直接表达这种关系是困难的。 但是，我们可以用矩阵形式表示这种关系，也可以按照黄金比例公式：

与

这将介绍matrix()和rint()函数。 matrix()函数创建矩阵， 数字四舍五入到最接近的整数，但结果不是整数。

实战时间 – 计算斐波纳契数

矩阵可以表示斐波那契递归关系。 我们可以将斐波纳契数的计算表示为重复的矩阵乘法：

1. 如下创建斐波那契矩阵：

   F = np.matrix([[1, 1], [1, 0]])
   print("F", F)
   

   斐波那契矩阵如下所示：

   F [[1 1]
    [1 0]]
   
   

2. 通过从 8 减去 1 并取矩阵的幂，计算出第 8 个斐波那契数（忽略 0）。 斐波那契数然后出现在对角线上：

   print("8th Fibonacci", (F ** 7)[0, 0])
   

   斐波那契数如下：

   8th Fibonacci 21
   
   

3. “黄金比例”公式（又称为“Binet”公式）使我们能够计算斐波纳契数，并在最后进行四舍五入。 计算前八个斐波那契数：

   n = np.arange(1, 9)
   sqrt5 = np.sqrt(5)
   phi = (1 + sqrt5)/2
   fibonacci = np.rint((phi**n - (-1/phi)**n)/sqrt5)
   print("Fibonacci", fibonacci)
   

   前八个斐波那契数如下：

   Fibonacci [  1\.   1\.   2\.   3\.   5\.   8\.  13\.  21.]
   
   

刚刚发生了什么？

我们用两种方法计算了斐波那契数。 在此过程中，我们了解了用于创建矩阵的matrix()函数。 我们还了解了rint()函数，该函数将数字四舍五入到最接近的整数，但不将类型更改为整数（请参见fibonacci.py）：

from __future__ import print_function
import numpy as np

F = np.matrix([[1, 1], [1, 0]])
print("F", F)
print("8th Fibonacci", (F ** 7)[0, 0])
n = np.arange(1, 9)

sqrt5 = np.sqrt(5)
phi = (1 + sqrt5)/2
fibonacci = np.rint((phi**n - (-1/phi)**n)/sqrt5)
print("Fibonacci", fibonacci)

勇往直前 – 时间计算

您可能想知道哪种方法更快，因此请继续并确定时间。 用frompyfunc()创建通用的斐波那契函数，并对其计时。

利萨如曲线

所有标准的三角函数，例如sin，cos，tan等，都由 NumPy 中的通用函数表示）。利萨如曲线是使用三角函数的一种有趣方式。 我记得在物理实验室的示波器上制作了李沙育的数字。 两个参数方程式描述了这些图形：

x = A sin(at + π/2)
y = B sin(bt)

实战时间 – 绘制利萨如曲线

利萨如的数字是由A，B，a和b四个参数确定的 。 为了简单起见，我们将A和B设置为1：

1. 将具有linspace()函数的t从-pi初始化为具有201点的pi：

   a = 9
   b = 8
   t = np.linspace(-np.pi, np.pi, 201)
   

2. 使用sin()函数和np.pi计算x：

   x = np.sin(a * t + np.pi/2)
   

3. 使用sin()函数计算y：

   y = np.sin(b * t)
   

4. 如下图所示：

   plt.plot(x, y)
   plt.title('Lissajous curves')
   plt.grid()
   plt.show()
   

   a = 9和b = 8的结果如下：

   

刚刚发生了什么？

我们用上述参数方程式绘制了利萨如曲线，其中A=B=1，a=9和b=8。 我们使用了sin()和linspace()函数，以及 NumPy pi常量（请参见lissajous.py）：

import numpy as np
import matplotlib.pyplot as plt

a = 9
b = 8
t = np.linspace(-np.pi, np.pi, 201)
x = np.sin(a * t + np.pi/2)
y = np.sin(b * t)
plt.plot(x, y)
plt.title('Lissajous curves')
plt.grid()
plt.show()

方波

方波也是您可以在示波器上查看的那些整洁的东西之一。 正弦波可以很好地将其近似为 。 毕竟，方波是可以用无限傅立叶级数表示的信号。

注意

傅立叶级数是以著名数学家让·巴蒂斯特·傅立叶（Jean-Baptiste Fourier）命名的，一系列正弦和余弦项之和。

代表方波的该特定系列的公式如下：

实战时间 – 绘制方波

就像上一节一样，我们将初始化t。 我们需要总结一些术语。 术语数量越多，结果越准确； k = 99应该足够。 为了绘制方波，请按照下列步骤操作：

1. 我们将从初始化t和k开始。 将该函数的初始值设置为0：

   t = np.linspace(-np.pi, np.pi, 201)
   k = np.arange(1, 99)
   k = 2 * k - 1
   f = np.zeros_like(t)
   

2. 使用sin()和sum()函数计算函数值：

   for i, ti in enumerate(t):
      f[i] = np.sum(np.sin(k * ti)/k)
   
   f = (4 / np.pi) * f
   

3. 要绘制的代码与上一节中的代码几乎相同：

   plt.plot(t, f)
   plt.title('Square wave')
   plt.grid()
   plt.show()
   

   用k = 99生成的所得方波如下：

   

刚刚发生了什么？

我们使用 sin()函数生成了方波，或者至少是它的近似值。 输入值通过 linspace()函数进行组装，而k值通过arange()函数进行组装（请参见squarewave.py）：

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(-np.pi, np.pi, 201)
k = np.arange(1, 99)
k = 2 * k - 1
f = np.zeros_like(t)

for i, ti in enumerate(t):
   f[i] = np.sum(np.sin(k * ti)/k)

f = (4 / np.pi) * f

plt.plot(t, f)
plt.title('Square wave')
plt.grid()
plt.show()

勇往直前 – 摆脱循环

您可能已经注意到代码中存在一个循环。 使用 NumPy 函数摆脱它，并确保性能也得到改善。

锯齿波和三角波

锯齿波和三角波也是在示波器上容易看到的现象。 就像方波一样，我们可以定义无限傅立叶级数。 三角波可以通过获取锯齿波的绝对值来找到。 表示一系列锯齿波的公式如下：

实战时间 – 绘制锯齿波和三角波

就像上一节一样，我们初始化t。 同样，k = 99应该足够。 为了绘制锯齿波和三角波，请按照下列步骤操作：

1. 将该函数的初始值设置为zero：

   t = np.linspace(-np.pi, np.pi, 201)
   k = np.arange(1, 99)
   f = np.zeros_like(t)
   

2. 使用sin()和sum()函数计算函数值：

   for i, ti in enumerate(t):
      f[i] = np.sum(np.sin(2 * np.pi * k * ti)/k)
   
   f = (-2 / np.pi) * f
   

3. 绘制锯齿波和三角波很容易，因为三角波的值应等于锯齿波的绝对值。 绘制波形，如下所示：

   plt.plot(t, f, lw=1.0, label='Sawtooth')
   plt.plot(t, np.abs(f), '--', lw=2.0, label='Triangle')
   plt.title('Triangle and sawtooth waves')
   plt.grid()
   plt.legend(loc='best')
   plt.show()
   

   在下图中，三角形波是带有虚线的波：

   

刚刚发生了什么？

我们使用sin()函数绘制了锯齿波 。 我们将输入值与linspace()函数组合在一起，并将k值与arange()函数组合在一起。 通过取绝对值从锯齿波中产生一个三角波（请参见sawtooth.py）：

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(-np.pi, np.pi, 201)
k = np.arange(1, 99)
f = np.zeros_like(t)

for i, ti in enumerate(t):
   f[i] = np.sum(np.sin(2 * np.pi * k * ti)/k)

f = (-2 / np.pi) * f
plt.plot(t, f, lw=1.0, label='Sawtooth')
plt.plot(t, np.abs(f), '--', lw=2.0, label='Triangle')
plt.title('Triangle and sawtooth waves')
plt.grid()
plt.legend(loc='best')
plt.show()

勇往直前 – 摆脱循环

如果您选择接受 ，那么您面临的挑战是摆脱程序中的循环。 它应该可以与 NumPy 函数一起使用，并且性能应该得到改善。

按位和比较函数

按位函数是整数或整数数组的位，因为它们是通用函数。 运算符^，&，|，<<，>>等具有其 NumPy 对应物。 比较运算符，例如<，>和==等也是如此。 这些运算符使您可以做一些巧妙的技巧，从而提高性能； 但是，它们会使您的代码难以理解，因此请谨慎使用。

实战时间 – 翻转位

现在，我们将介绍三个技巧：检查整数的符号是​​否不同，检查数字是否为2的幂，以及计算作为2的幂的数字的模数。 我们将展示一个仅用于运算符的符号，以及一个使用相应的 NumPy 函数的符号：

1. 第一个技巧取决于XOR或^运算符。 XOR运算符也称为不等式运算符； 因此，如果两个操作数的符号位不同，则XOR运算将导致负数。

   下面的真值表说明了XOR运算符：

   输入 1	输入 2	异或
   True	True	False
   False	True	True
   True	False	True
   False	False	False

   ^运算符对应于bitwise_xor()函数，<运算符对应于less()函数：

   x = np.arange(-9, 9)
   y = -x
   print("Sign different?", (x ^ y) < 0)
   print("Sign different?", np.less(np.bitwise_xor(x, y), 0))
   

   结果为 ，如下所示：

   Sign different? [ True  True  True  True  True  True  True  True  True False  True  True
    True  True  True  True  True  True]
   Sign different? [ True  True  True  True  True  True  True  True  True False  True  True
    True  True  True  True  True  True]
   
   

   正如预期的那样，除零以外，所有符号均不同。

2. 2 的幂由 1 表示，后跟一系列二进制表示的尾随零。 例如，10，100或1000。 比 2 的幂小 1 的数字将由一排二进制 1 表示。 例如，11，111或1111（或十进制中的3，7和15）。 现在，如果我们对 2 的幂，和比它小 1 的整数进行“与”运算，则应该得到 0。

   AND运算符的真值表如下所示：

   输入 1	输入 2	AND
   True	True	True
   False	True	False
   True	False	False
   False	False	False

   &的 NumPy 对应项是bitwise_and()，==的对应项是equal()通用函数：

   print("Power of 2?\n", x, "\n", (x & (x - 1)) == 0)
   print("Power of 2?\n", x, "\n", np.equal(np.bitwise_and(x,  (x - 1)), 0))
   

   The result is shown as follows:

   Power of 2?
   **[-9 -8 -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8]** 
   [False False False False False False False False False  True  True  True
    False  True False False False  True]
   Power of 2?
   **[-9 -8 -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8]** 
   [False False False False False False False False False  True  True  True
    False  True False False False  True]
   
   

3. 当计算整数的 2 的幂的模数时，例如 4、8、16 等，计算 4 的模数的技巧实际上起作用。 左移导致值加倍。我们在上一步中看到，从 2 的幂中减去 1 会导致二进制表示形式的数字带有一行诸如 11、111 或 1111 之类的数字。 这基本上给了我们一个掩码。 用这样的数字按位与，得到的余数为 2。 NumPy 的<<的等价物是left_shift()通用函数：

   print("Modulus 4\n", x, "\n", x & ((1 << 2) - 1))
   print("Modulus 4\n", x, "\n", np.bitwise_and(x, np.left_shift(1, 2) - 1))
   

   The result is shown as follows:

   Modulus 4
   **[-9 -8 -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8]** 
   [3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0]
   Modulus 4
   **[-9 -8 -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8]** 
   [3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0]
   
   

刚刚发生了什么？

我们介绍了三点技巧：检查整数的符号是​​否不同，检查数字是否为2的幂，并计算数字的模数为2的幂。 我们看到了运算符^，&，<<和<的 NumPy 对应项（请参见bittwidling.py）：

from __future__ import print_function
import numpy as np

x = np.arange(-9, 9)
y = -x
print("Sign different?", (x ^ y) < 0)
print("Sign different?", np.less(np.bitwise_xor(x, y), 0))
print("Power of 2?\n", x, "\n", (x & (x - 1)) == 0)
print("Power of 2?\n", x, "\n", np.equal(np.bitwise_and(x,  (x - 1)), 0))
print("Modulus 4\n", x, "\n", x & ((1 << 2) - 1))
print("Modulus 4\n", x, "\n", np.bitwise_and(x, np.left_shift(1, 2) - 1))

花式索引

at()方法是在 NumPy 1.8 中添加的 。 此方法允许原地建立花式索引。 花式索引是不涉及整数或切片的索引 ，这是正常的索引。 原地意味着将对我们操作的数组进行修改。

at()方法的签名为ufunc.at(a, indices[, b])。  indices数组指定要操作的元素。 我们仅需要b数组用于具有两个操作数的通用函数。 以下“实战时间”部分给出了at()方法的示例 。

实战时间 – 使用at()方法为 ufuncs 原地建立索引

要演示方法的工作方式，请启动 Python 或 IPython shell 并导入 NumPy。 您现在应该知道如何执行此操作。

1. 创建一个由七个随机整数组成的数组，该整数从-3到3，种子为42：

   >>> a = np.random.random_integers(-3, 3, 7)
   >>> a
   array([ 1,  0, -1,  2,  1, -2,  0])
   
   

   当我们在编程中谈论随机数字时，我们通常会谈论伪随机数。 这些数字看起来是随机的，但实际上是使用种子来计算的。

2. 将sign()通用函数的at()方法应用于第四和第六个数组元素：

   >>> np.sign.at(a, [3, 5])
   >>> a
   array([ 1, 0, -1,  1,  1, -1,  0])
   
   

刚刚发生了什么？

我们使用at()方法来选择数组元素，并执行原地操作-确定符号。 我们还学习了如何创建随机整数。

总结

在本章中，您学习了关于矩阵和通用函数的知识。 我们介绍了如何创建矩阵，并研究了通用函数如何工作。 您简要介绍了算术，三角函数，按位和比较通用函数。

在下一章中，您将介绍 NumPy 模块。

六、深入探索 NumPy 模块

NumPy 具有许多从其前身 Numeric 继承的模块。 其中一些包具有 SciPy 对应版本，可能具有更完整的功能。 我们将在下一章中讨论 SciPy。

在本章中，我们将介绍以下主题：

• linalg包
• fft包
• 随机数
• 连续和离散分布

线性代数

线性代数是数学的重要分支。 numpy.linalg包包含线性代数函数。 使用此模块，您可以求矩阵求逆，计算特征值，求解线性方程式和确定行列式等。

实战时间 – 转换矩阵

线性代数中矩阵A的逆是矩阵A^(-1)，当与原始矩阵相乘时，它等于单位矩阵I。 可以这样写：

A A^(-1) = I

numpy.linalg包中的inv()函数可以通过以下步骤反转示例矩阵：

1. 使用前面章节中使用的mat()函数创建示例矩阵：

   A = np.mat("0 1 2;1 0 3;4 -3 8")
   print("A\n", A)
   

   A矩阵如下所示：

   A
   [[ 0  1  2]
    [ 1  0  3]
    [ 4 -3  8]]
   
   

2. 用inv()函数将矩阵求逆：

   inverse = np.linalg.inv(A)
   print("inverse of A\n", inverse)
   

   逆矩阵如下所示：

   inverse of A
   [[-4.5  7\.  -1.5]
    [-2\.   4\.  -1\. ]
    [ 1.5 -2\.   0.5]]
   
   

   提示

   如果矩阵是奇异的，或者不是正方形，则引发LinAlgError。 如果需要，可以用笔和纸手动检查结果。 这留给读者练习。

3. 通过将原始矩阵乘以inv()函数的结果来检查结果：

   print("Check\n", A * inverse)
   

   结果是单位矩阵，如预期的那样：

   Check
   [[ 1\.  0\.  0.]
    [ 0\.  1\.  0.]
    [ 0\.  0\.  1.]]
   
   

刚刚发生了什么？

我们用numpy.linalg包的inv()函数计算了矩阵的逆。 我们使用矩阵乘法检查了这是否确实是逆矩阵（请参见inversion.py）：

from __future__ import print_function
import numpy as np

A = np.mat("0 1 2;1 0 3;4 -3 8")
print("A\n", A)

inverse = np.linalg.inv(A)
print("inverse of A\n", inverse)

print("Check\n", A * inverse)

小测验 - 创建矩阵

Q1. 哪个函数可以创建矩阵？

1. array
2. create_matrix
3. mat
4. vector

勇往直前 – 反转自己的矩阵

创建自己的矩阵并将其求逆。 逆仅针对方阵定义。 矩阵必须是正方形且可逆； 否则，将引发LinAlgError异常。

求解线性系统

矩阵以线性方式将向量转换为另一个向量。 该变换在数学上对应于线性方程组。 numpy.linalg函数solve()求解形式为Ax = b的线性方程组，其中A是矩阵，b可以是一维或二维数组，而x是未知数变量。 我们将看到dot()函数的使用。 此函数返回两个浮点数组的点积。

dot()函数计算点积。 对于矩阵A和向量b，点积等于以下总和：

实战时间 – 解决线性系统

通过以下步骤解决线性系统的示例：

1. 创建A和b：

   A = np.mat("1 -2 1;0 2 -8;-4 5 9")
   print("A\n", A)
   b = np.array([0, 8, -9])
   print("b\n", b)
   

   A和b出现如下：

   

2. 用solve()函数求解线性系统：

   x = np.linalg.solve(A, b)
   print("Solution", x)
   

   线性系统的解如下：

   Solution [ 29\.  16\.   3.]
   
   

3. 使用dot()函数检查解决方案是否正确：

   print("Check\n", np.dot(A , x))
   

   结果是预期的：

   Check
   [[ 0\.  8\. -9.]]
   
   

刚刚发生了什么？

我们使用 NumPy linalg模块的solve()函数求解了线性系统，并使用dot()函数检查了解。 请参考本书代码捆绑中的solution.py文件：

from __future__ import print_function
import numpy as np

A = np.mat("1 -2 1;0 2 -8;-4 5 9")
print("A\n", A)

b = np.array([0, 8, -9])
print("b\n", b)

x = np.linalg.solve(A, b)
print("Solution", x)

print("Check\n", np.dot(A , x))

查找特征值和特征向量

特征值是方程Ax = ax的标量解，其中A是二维矩阵，x是一维向量。特征向量是对应于特征值的向量。 numpy.linalg包中的eigvals()函数计算特征值。 eig()函数返回包含特征值和特征向量的元组。

实战时间 – 确定特征值和特征向量

让我们计算矩阵的特征值：

1. 创建一个矩阵，如下所示：

   A = np.mat("3 -2;1 0")
   print("A\n", A)
   

   我们创建的矩阵如下所示：

   A
   [[ 3 -2]
    [ 1  0]]
   
   

2. 调用eigvals()函数：

   print("Eigenvalues", np.linalg.eigvals(A))
   

   矩阵的特征值如下：

   Eigenvalues [ 2\.  1.]
   
   

3. 使用eig（）函数确定特征值和特征向量。 此函数返回一个元组，其中第一个元素包含特征值，第二个元素包含相应的特征向量，按列排列：

   eigenvalues, eigenvectors = np.linalg.eig(A)
   print("First tuple of eig", eigenvalues)
   print("Second tuple of eig\n", eigenvectors)
   

   特征值和特征向量如下所示：

   First tuple of eig [ 2\.  1.]
   Second tuple of eig
   [[ 0.89442719  0.70710678]
    [ 0.4472136   0.70710678]]
   
   

4. 通过计算特征值方程Ax = ax的右侧和左侧，使用dot()函数检查结果：

   for i, eigenvalue in enumerate(eigenvalues):
         print("Left", np.dot(A, eigenvectors[:,i]))
         print("Right", eigenvalue * eigenvectors[:,i])
         print()
   

   输出如下：

   Left [[ 1.78885438]
    [ 0.89442719]]
   Right [[ 1.78885438]
    [ 0.89442719]]
   
   

刚刚发生了什么？

我们发现了具有numpy.linalg模块的eigvals()和eig()函数的矩阵的特征值和特征向量。 我们使用dot()函数检查了结果（请参见eigenvalues.py）：

from __future__ import print_function
import numpy as np

A = np.mat("3 -2;1 0")
print("A\n", A)

print("Eigenvalues", np.linalg.eigvals(A) )

eigenvalues, eigenvectors = np.linalg.eig(A)
print("First tuple of eig", eigenvalues)
print("Second tuple of eig\n", eigenvectors)

for i, eigenvalue in enumerate(eigenvalues):
      print("Left", np.dot(A, eigenvectors[:,i]))
      print("Right", eigenvalue * eigenvectors[:,i])
      print()

奇异值分解

奇异值分解（SVD）是一种分解因子，可以将矩阵分解为三个矩阵的乘积。 SVD 是先前讨论的特征值分解的概括。 SVD 对于像这样的伪逆算法非常有用，我们将在下一部分中进行讨论。 numpy.linalg包中的 svd()函数可以执行此分解。 此函数返回三个矩阵U，∑和V，使得U和V为一元且∑包含输入矩阵的奇异值：

星号表示 Hermitian 共轭或共轭转置。复数的共轭改变复数虚部的符号，因此与实数无关。

注意

如果A*A = AA* = I（单位矩阵），则复方矩阵 A 是单位的。 我们可以将 SVD 解释为三个操作的序列-旋转，缩放和另一个旋转。

我们已经在本书中转置了矩阵。 转置翻转矩阵，将行变成列，然后将列变成行。

实战时间 – 分解矩阵

现在该使用以下步骤用 SVD 分解矩阵：

1. 首先，创建如下所示的矩阵：

   A = np.mat("4 11 14;8 7 -2")
   print("A\n", A)
   

   我们创建的矩阵如下所示：

   A
   [[ 4 11 14]
    [ 8  7 -2]]
   
   

2. 用svd()函数分解矩阵：

   U, Sigma, V = np.linalg.svd(A, full_matrices=False)
   print("U")
   print(U)
   print("Sigma")
   print(Sigma)
   print("V")
   print(V)
   

   由于full_matrices=False规范，NumPy 执行了简化的 SVD 分解，计算速度更快。 结果是一个元组，在左侧和右侧分别包含两个单位矩阵U和V，以及中间矩阵的奇异值：

   U
   [[-0.9486833  -0.31622777]
    [-0.31622777  0.9486833 ]]
   Sigma
   [ 18.97366596   9.48683298]
   V
   [[-0.33333333 -0.66666667 -0.66666667]
    [ 0.66666667  0.33333333 -0.66666667]]
   
   

3. 实际上，我们没有中间矩阵，只有对角线值。 其他值均为 0。 用diag()函数形成中间矩阵。 将三个矩阵相乘如下：

   print("Product\n", U * np.diag(Sigma) * V)
   

   这三个矩阵的乘积等于我们在第一步中创建的矩阵：

   Product
   [[  4\.  11\.  14.]
    [  8\.   7\.  -2.]]
   
   

刚刚发生了什么？

我们分解矩阵，并通过矩阵乘法检查结果。 我们使用了 NumPy linalg模块中的svd()函数（请参见decomposition.py）：

from __future__ import print_function
import numpy as np

A = np.mat("4 11 14;8 7 -2")
print("A\n", A)

U, Sigma, V = np.linalg.svd(A, full_matrices=False)

print("U")
print(U)

print("Sigma")
print(Sigma)

print("V")
print(V)

print("Product\n", U * np.diag(Sigma) * V)

伪逆

矩阵的 Moore-Penrose 伪逆的计算公式为numpy.linalg模块的pinv()函数。 使用 SVD 计算伪逆（请参见前面的示例）。 inv()函数仅接受方阵； pinv()函数确实没有的限制，因此被认为是反函数的推广。

实战时间 – 计算矩阵的伪逆

让我们计算矩阵的伪逆：

1. 首先，创建一个矩阵：

   A = np.mat("4 11 14;8 7 -2")
   print("A\n", A)
   

   我们创建的矩阵如下所示：

   A
   [[ 4 11 14]
    [ 8  7 -2]]
   
   

2. 用pinv()函数计算伪逆矩阵：

   pseudoinv = np.linalg.pinv(A)
   print("Pseudo inverse\n", pseudoinv)
   

   伪逆结果如下：

   Pseudo inverse
   [[-0.00555556  0.07222222]
    [ 0.02222222  0.04444444]
    [ 0.05555556 -0.05555556]]
   
   

3. 将原始和伪逆矩阵相乘：

   print("Check", A * pseudoinv)
   

   我们得到的不是一个恒等矩阵，但是很接近它：

   Check [[  1.00000000e+00   0.00000000e+00]
    [  8.32667268e-17   1.00000000e+00]]
   
   

刚刚发生了什么？

我们使用numpy.linalg模块的pinv()函数计算了矩阵的伪逆。 通过矩阵乘法检查得出的矩阵大约是单位矩阵（请参见pseudoinversion.py）：

from __future__ import print_function
import numpy as np

A = np.mat("4 11 14;8 7 -2")
print("A\n", A)

pseudoinv = np.linalg.pinv(A)
print("Pseudo inverse\n", pseudoinv)

print("Check", A * pseudoinv)

行列式

行列式是与方阵相关的值。 在整个数学中都使用它； 有关更多详细信息，请参见这里。 对于n x n实值矩阵，行列式对应于矩阵变换后 n 维体积所经历的缩放。 行列式的正号表示体积保留它的方向（顺时针或逆时针），而负号表示方向相反。 numpy.linalg模块具有det()函数，该函数返回矩阵的行列式。

实战时间 – 计算矩阵的行列式

要计算矩阵的行列式 ，请按照下列步骤操作：

1. 创建矩阵：

   A = np.mat("3 4;5 6")
   print("A\n", A)
   

   我们创建的矩阵如下所示：

   A
   [[ 3\.  4.]
    [ 5\.  6.]]
   
   

2. 用det()函数计算行列式：

   print("Determinant", np.linalg.det(A))
   

   行列式如下所示：

   Determinant -2.0
   
   

刚刚发生了什么？

我们从numpy.linalg模块（请参见determinant.py）使用det()函数计算了矩阵的行列式：

from __future__ import print_function
import numpy as np

A = np.mat("3 4;5 6")
print("A\n", A)

print("Determinant", np.linalg.det(A))

快速傅立叶变换

快速傅里叶变换（FFT）是一种用于计算离散傅立叶变换（DFT）的有效算法。

注意

傅立叶变换与傅立叶级数相关，在上一章中提到了第 5 章，“处理矩阵和函数”。 傅里叶级数将信号表示为正弦和余弦项之和。

FFT 在更多朴素算法上进行了改进，其阶数为O(N log N)。 DFT 在信号处理，图像处理，求解偏微分方程等方面具有应用。 NumPy 有一个名为fft的模块，该模块提供 FFT 函数。 该模块中的许多函数已配对。 对于那些函数，另一个函数执行逆运算。 例如，fft()和ifft()函数形成这样的一对。

实战时间 – 计算傅立叶变换

首先，我们将创建一个要转换的信号。 通过以下步骤计算傅立叶变换：

1. 创建具有30点的余弦波，如下所示：

   x =  np.linspace(0, 2 * np.pi, 30)
   wave = np.cos(x)
   

2. 用fft()函数变换余弦波：

   transformed = np.fft.fft(wave)
   

3. 使用ifft()函数应用逆变换。 它应该大致返回原始信号。 检查以下行：

   print(np.all(np.abs(np.fft.ifft(transformed) - wave) < 10 ** -9))
   

   结果显示如下：

   True
   
   

4. 使用 matplotlib 绘制转换后的信号：

   plt.plot(transformed)
   plt.title('Transformed cosine')
   plt.xlabel('Frequency')
   plt.ylabel('Amplitude')
   plt.grid()
   plt.show()
   

   下图显示了 FFT 结果：

   

刚刚发生了什么？

我们将fft()函数应用于余弦波。 应用ifft()函数后，我们得到了信号（请参阅fourier.py）：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

x =  np.linspace(0, 2 * np.pi, 30)
wave = np.cos(x)
transformed = np.fft.fft(wave)
print(np.all(np.abs(np.fft.ifft(transformed) - wave) < 10 ** -9))

plt.plot(transformed)
plt.title('Transformed cosine')
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.grid()
plt.show()

移位

numpy.linalg模块的fftshift()函数将零频率分量移到频谱中心。 零频率分量对应于信号的平均值 。 ifftshift()函数可逆转此操作。

实战时间 – 变换频率

我们将创建一个信号，对其进行转换，然后将其移位。 按以下步骤移动频率：

1. 创建具有30点的余弦波：

   x =  np.linspace(0, 2 * np.pi, 30)
   wave = np.cos(x)
   

2. 使用fft()函数变换余弦波：

   transformed = np.fft.fft(wave)
   

3. 使用fftshift()函数移动信号：

   shifted = np.fft.fftshift(transformed)
   

4. 用ifftshift()函数反转移位。 这应该撤消这种转变。 检查以下代码段：

   print(np.all((np.fft.ifftshift(shifted) - transformed) < 10 ** -9))
   

   The result appears as follows:

   True
   
   

5. 绘制信号并使用 matplotlib 对其进行转换：

   plt.plot(transformed, lw=2, label="Transformed")
   plt.plot(shifted, '--', lw=3, label="Shifted")
   plt.title('Shifted and transformed cosine wave')
   plt.xlabel('Frequency')
   plt.ylabel('Amplitude')
   plt.grid()
   plt.legend(loc='best')
   plt.show()
   

   下图显示了移位和 FFT 的效果：

   

刚刚发生了什么？

我们将fftshift()函数应用于余弦波。 应用ifftshift()函数后，我们返回我们的信号（请参阅fouriershift.py）：

import numpy as np
import matplotlib.pyplot as plt

x =  np.linspace(0, 2 * np.pi, 30)
wave = np.cos(x)
transformed = np.fft.fft(wave)
shifted = np.fft.fftshift(transformed)
print(np.all(np.abs(np.fft.ifftshift(shifted) - transformed) < 10 ** -9))

plt.plot(transformed, lw=2, label="Transformed")
plt.plot(shifted, '--', lw=3, label="Shifted")
plt.title('Shifted and transformed cosine wave')
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.grid()
plt.legend(loc='best')
plt.show()

随机数

蒙特卡罗方法，随机演算等中使用了随机数。 真正的随机数很难生成，因此在实践中，我们使用伪随机数字，除了某些非常特殊的情况外，对于大多数意图和目的来说都是足够随机的。 这些数字似乎是随机的，但是如果您更仔细地分析它们，则将意识到它们遵循一定的模式。 与随机数相关的函数位于 NumPy 随机模块中。 核心随机数字生成器基于 Mersenne Twister 算法 -- 一种标准且众所周知的算法。 我们可以从离散或连续分布中生成随机数。 分布函数具有一个可选的size参数，该参数告诉 NumPy 生成多少个数字。 您可以指定整数或元组作为大小。 这将导致数组中填充适当形状的随机数。 离散分布包括几何分布，超几何分布和二项分布。

实战时间 – 使用二项来赌博

二项分布模型是整数个独立试验中的成功的次数，其中每个实验中成功的概率是固定的数字。

想象一下一个 17 世纪的赌场，您可以在上面掷 8 个筹码。 九枚硬币被翻转。 如果少于五个，那么您将损失八分之一，否则将获胜。 让我们模拟一下，从拥有的 1,000 个硬币开始。 为此，可以使用随机模块中的binomial()函数。

要了解 binomial()函数，请查看以下部分：

1. 将代表现金余额的数组初始化为零。 调用大小为 10000 的binomial()函数。这表示在我们的赌场中有 10,000 次硬币翻转：

   cash = np.zeros(10000)
   cash[0] = 1000
   outcome = np.random.binomial(9, 0.5, size=len(cash))
   

2. 查看硬币翻转的结果并更新现金数组。 打印结果的最小值和最大值，只是为了确保我们没有任何奇怪的异常值：

   for i in range(1, len(cash)):
      if outcome[i] < 5:
         cash[i] = cash[i - 1] - 1
      elif outcome[i] < 10:
         cash[i] = cash[i - 1] + 1
      else:
         raise AssertionError("Unexpected outcome " + outcome)
   
   print(outcome.min(), outcome.max())
   

   不出所料，该值在0和9之间。 在下图中，您可以看到现金余额执行随机游走：

   

刚刚发生了什么？

我们使用 NumPy 随机模块中的binomial()函数进行了随机游走实验（请参见headortail.py）：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

cash = np.zeros(10000)
cash[0] = 1000
np.random.seed(73)
outcome = np.random.binomial(9, 0.5, size=len(cash))

for i in range(1, len(cash)):
   if outcome[i] < 5:
      cash[i] = cash[i - 1] - 1
   elif outcome[i] < 10:
      cash[i] = cash[i - 1] + 1
   else:
      raise AssertionError("Unexpected outcome " + outcome)

print(outcome.min(), outcome.max())

plt.plot(np.arange(len(cash)), cash)
plt.title('Binomial simulation')
plt.xlabel('# Bets')
plt.ylabel('Cash')
plt.grid()
plt.show()

超几何分布

超几何分布对其中装有两种对象的罐进行建模。 该模型告诉我们，如果我们从罐子中取出指定数量的物品而不更换它们，可以得到一种类型的对象。 NumPy 随机模块具有模拟这种情况的hypergeometric()函数。

实战时间 – 模拟游戏节目

想象一下，游戏会显示出参赛者每次正确回答问题时，都会从罐子中拉出三个球，然后放回去。 现在，有一个陷阱，罐子里的一个球不好。 每次拔出时，参赛者将失去 6 分。 但是，如果他们设法摆脱 25 个普通球中的 3 个，则得到 1 分。 那么，如果我们总共有 100 个问题，将会发生什么？ 查看以下部分以了解解决方案：

1. 使用hypergeometric()函数初始化游戏结果。 此函数的第一个参数是做出正确选择的方法数量，第二个参数是做出错误选择的方法数量，第三个参数是采样的项目数量：

   points = np.zeros(100)
   outcomes = np.random.hypergeometric(25, 1, 3, size=len(points))
   

2. 根据上一步的结果设置评分：

   for i in range(len(points)):
      if outcomes[i] == 3:
         points[i] = points[i - 1] + 1
      elif outcomes[i] == 2:
         points[i] = points[i - 1] - 6
      else:
         print(outcomes[i])
   

   下图显示了评分如何演变：

   

刚刚发生了什么？

我们使用 NumPy random模块中的 hypergeometric()函数模拟了游戏节目。 游戏得分取决于每次比赛参与者从罐子中抽出多少好球和坏球（请参阅urn.py）：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

points = np.zeros(100)
np.random.seed(16)
outcomes = np.random.hypergeometric(25, 1, 3, size=len(points))

for i in range(len(points)):
   if outcomes[i] == 3:
      points[i] = points[i - 1] + 1
   elif outcomes[i] == 2:
      points[i] = points[i - 1] - 6
   else:
      print(outcomes[i])

plt.plot(np.arange(len(points)), points)
plt.title('Game show simulation')
plt.xlabel('# Rounds')
plt.ylabel('Score')
plt.grid()
plt.show()

连续分布

我们通常使用概率密度函数（PDF）对连续分布进行建模。 值处于特定间隔的可能性由 PDF 的积分确定）。 NumPy random模块具有表示连续分布的函数-beta()，chisquare()，exponential()，f()，gamma()，gumbel()，laplace()，lognormal() ，logistic()，multivariate_normal()，noncentral_chisquare()，noncentral_f()，normal()等。

实战时间 – 绘制正态分布

我们可以从正态分布中生成随机数，并通过直方图可视化其分布）。 通过以下步骤绘制正态分布：

1. 使用random NumPy 模块中的normal()函数，为给定的样本量生成随机的数字：

   N=10000
   normal_values = np.random.normal(size=N)
   

2. 绘制直方图和理论 PDF，其中心值为 0，标准偏差为 1。 为此，请使用 matplotlib：

   _, bins, _ = plt.hist(normal_values, np.sqrt(N), normed=True, lw=1)
   sigma = 1
   mu = 0
   plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins - mu)**2 / (2 * sigma**2) ),lw=2)
   plt.show()
   

   在下面的图表中，我们看到了熟悉的钟形曲线：

   

刚刚发生了什么？

我们使用来自随机 NumPy 模块的normal()函数可视化正态分布。 为此，我们绘制了钟形曲线和随机生成的值的直方图（请参见normaldist.py）：

import numpy as np
import matplotlib.pyplot as plt

N=10000

np.random.seed(27)
normal_values = np.random.normal(size=N)
_, bins, _ = plt.hist(normal_values, np.sqrt(N), normed=True, lw=1, label="Histogram")
sigma = 1
mu = 0
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins - mu)**2 / (2 * sigma**2) ), '--', lw=3, label="PDF")
plt.title('Normal distribution')
plt.xlabel('Value')
plt.ylabel('Normalized Frequency')
plt.grid()
plt.legend(loc='best')
plt.show()

对数正态分布

对数正态分布是自然对数呈正态分布的随机变量的分布。 随机 NumPy 模块的 lognormal()函数可对该分布进行建模。

实战时间 – 绘制对数正态分布

让我们用直方图可视化对数正态分布及其 PDF：

1. 使用random NumPy 模块中的normal()函数生成随机数：

   N=10000
   lognormal_values = np.random.lognormal(size=N)
   

2. 绘制直方图和理论 PDF，其中心值为 0，标准偏差为 1：

   _, bins, _ = plt.hist(lognormal_values, np.sqrt(N), normed=True, lw=1)
   sigma = 1
   mu = 0
   x = np.linspace(min(bins), max(bins), len(bins))
   pdf = np.exp(-(numpy.log(x) - mu)**2 / (2 * sigma**2))/ (x * sigma * np.sqrt(2 * np.pi))
   plt.plot(x, pdf,lw=3)
   plt.show()
   

   直方图和理论 PDF 的拟合非常好，如下图所示：

   

刚刚发生了什么？

我们使用random NumPy 模块中的lognormal()函数可视化了对数正态分布。 我们通过绘制理论 PDF 曲线和随机生成的值的直方图（请参见lognormaldist.py）来做到这一点：

import numpy as np
import matplotlib.pyplot as plt

N=10000
np.random.seed(34)
lognormal_values = np.random.lognormal(size=N)
_, bins, _ = plt.hist(lognormal_values, np.sqrt(N), normed=True, lw=1, label="Histogram")
sigma = 1
mu = 0
x = np.linspace(min(bins), max(bins), len(bins))
pdf = np.exp(-(np.log(x) - mu)**2 / (2 * sigma**2))/ (x * sigma * np.sqrt(2 * np.pi))
plt.xlim([0, 15])
plt.plot(x, pdf,'--', lw=3, label="PDF")
plt.title('Lognormal distribution')
plt.xlabel('Value')
plt.ylabel('Normalized frequency')
plt.grid()
plt.legend(loc='best')
plt.show()

统计量自举

自举是一种用于估计方差，准确性和其他样本估计量度的方法，例如算术平均值。 最简单的自举过程包括以下步骤：

1. 从具有相同大小N的原始数据样本中生成大量样本。 您可以将原始数据视为包含数字的罐子。 我们通过N次从瓶子中随机选择一个数字来创建新样本。 每次我们将数字返回到罐子中时，一个生成的样本中可能会多次出现一个数字。
2. 对于新样本，我们为每个样本计算要调查的统计估计值（例如，算术平均值）。 这为我们提供了估计器可能值的样本。

实战时间 – 使用numpy.random.choice()进行采样

我们将使用numpy.random.choice()函数对执行自举。

1. 启动 IPython 或 Python Shell 并导入 NumPy：

   $ ipython
   In [1]: import numpy as np
   
   

2. 按照正态分布生成数据样本：

   In [2]: N = 500
   
   In [3]: np.random.seed(52)
   
   In [4]: data = np.random.normal(size=N)
   
   

3. 计算数据的平均值：

   In [5]: data.mean()
   Out[5]: 0.07253250605445645
   
   

   从原始数据生成100样本并计算其平均值（当然，更多样本可能会导致更准确的结果）：

   In [6]: bootstrapped = np.random.choice(data, size=(N, 100))
   
   In [7]: means = bootstrapped.mean(axis=0)
   
   In [8]: means.shape
   Out[8]: (100,)
   
   

4. 计算得到的算术平均值的均值，方差和标准偏差：

   In [9]: means.mean()
   Out[9]: 0.067866373318115278
   
   In [10]: means.var()
   Out[10]: 0.001762807104774598
   
   In [11]: means.std()
   Out[11]: 0.041985796464692651
   
   

   如果我们假设均值的正态分布，则可能需要了解 z 得分，其定义如下：

   

   In [12]: (data.mean() - means.mean())/means.std()
   Out[12]: 0.11113598238549766
   
   

   从 z 得分值，我们可以了解实际均值的可能性。

刚刚发生了什么？

我们通过生成样本并计算每个样本的平均值来自举数据样本。 然后，我们计算了均值，标准差，方差和均值的 z 得分。 我们使用numpy.random.choice()函数进行自举。

总结

您在本章中学到了很多有关 NumPy 模块的知识。 我们介绍了线性代数，快速傅立叶变换，连续和离散分布以及随机数。

在下一章中，我们将介绍专门的例程。 这些函数可能不经常使用，但是在需要时非常有用。