27、Python 中的数学游戏与 NumPy 模块应用

Python 中的数学游戏与 NumPy 模块应用

在 Python 的世界里，有许多有趣的数学游戏和强大的工具可以帮助我们解决各种问题。本文将介绍一些使用 Python 进行数学模拟和数据处理的示例，同时深入探讨 NumPy 模块的使用。

随机模拟与概率问题

首先，我们来看一些随机模拟的例子。在磁盘头寻道问题中，我们可以使用随机数来模拟磁盘头的移动。

from random import random
N = 1000  # number of files the head seeks
tot_dist = 0
for i in range(N):
    tot_dist += random()*2
print(tot_dist/N)

在这个例子中，我们模拟了磁盘头的寻道过程，通过多次随机移动并计算总距离，最后得到平均移动距离。需要注意的是， N 的值越大，结果越准确。

另一种情况是磁盘头不回到位置 0，而是停留在当前位置等待下一次寻道。我们可以用以下代码来模拟：

from random import random
N = 1000  # number of files the head seeks
tot_dist, cur_loc = 0, 0
for i in range(N):
    new_loc = random()
    tot_dist += abs(cur_loc - new_loc)
    cur_loc = new_loc
print(tot_dist/N)

这个结果接近 1/3。

接下来是一个有趣的概率问题：朋友见面问题。假设两个朋友决定在晚上 8 点到 9 点之间见面，其中一个朋友到达后会等待 10 分钟。我们可以使用随机模块来模拟这个过程，并通过可视化来计算他们见面的概率。

from random import random
from pylab import *

N = 40000  # number of events

# generate N events of friends times
friend1, friend2 = [], []
for i in range(N):
    friend1.append(random())
    friend2.append(random())

# find all occurrences of friends meeting
met = array([(x, y) for (x, y) in zip(friend1, friend2) if abs(y - x) < 1.0/6])
not_met = array([(x, y) for (x, y) in zip(friend1, friend2) if abs(y - x) >= 1.0/6])

# plot the result, this might shed some light on the problem!
plot(met[:, 0], met[:, 1], 'y+', mec='y')
plot(not_met[:, 0], not_met[:, 1], 'bo', mec='b')
title("Probability of meeting: %1.3f" % (float(len(met))/N))
xlabel('Time of arrival of Friend 1 [minutes]')
ylabel('Time of arrival of Friend 2 [minutes]')
xticks([n/6 for n in range(7)], [n*10 for n in range(7)])
yticks([n/6 for n in range(7)], [n*10 for n in range(7)])
axis('scaled')
show()

通过可视化，我们可以看到中间的区域代表两个朋友见面的事件。通过计算这个区域的面积，我们可以得到见面的概率，理论值为 1 - (5/6)^2 = 11/36 ≈ 0.3055，与数值估计结果接近。

随机序列操作

随机模块还提供了一些操作序列的函数，如下表所示：

函数	描述
choice(s)	从序列 s 中返回一个随机元素
shuffle(s)	对序列 s 进行洗牌操作
sample(s, n)	从序列 s 中返回一个大小为 n 的子序列

我们可以用一个扑克牌的例子来演示这些函数的使用：

from random import choice, shuffle, sample
cards = list(zip(list(range(1, 14))*4, 'S'*13+'H'*13+'D'*13+'C'*13))
print(cards[:5])
print(choice(cards))
shuffle(cards)
print(cards[:5])
print(sample(cards, 5))

NumPy 模块介绍

NumPy 是 Python 中一个强大的数值计算库，其 ndarray 对象是许多数据处理和可视化脚本的基础。

数组创建

创建 NumPy 数组有多种方法，如下表所示：

函数	描述	示例
array(s)	基于序列 s 创建一个数组	array(((1, 2),(3, 4))) 返回 array([[1, 2], [3, 4]])
ones(t)	基于元组 t 创建一个用 1 初始化的 N 维数组	ones(2) 返回 array([ 1., 1.]) ； ones((3,3)) 返回 array([[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]])
arange([start,] stop[, step])	创建一个从 start 开始，到 stop 结束（不包括 stop ），步长为 step 的数组	arange(1, 2, 0.5) 返回 array([ 1., 1.5])
linspace(start, stop, num=50)	创建一个从 start 到 stop 线性间隔的大小为 num 的向量	linspace(1, 10, 3) 返回 array([ 1., 5.5, 10.])
logspace(start, stop, num=50)	类似于 linspace() ，但值在对数尺度上从 10^ start 到 10^ stop 均匀间隔	logspace(0, 1, 3) 返回 array([ 1., 3.16227766, 10.])

以下是一些创建数组的示例代码：

from numpy import *
v = array([1, 2])
print(v)
print(type(v))
m = array([[1, 0],[0, 4]])
print(m)
print(type(m))

数组切片、索引和重塑

数组可以使用 reshape() 和 resize() 函数进行重塑，使用 Python 的切片和索引运算符 [] 和 [:] 进行索引和切片。

a = arange(12).reshape(4, 3)
print(a)
print(a[1])
print(a[-1])
print(a[1, 1])
print(a[:, 1])
print(a[1, :2])

NumPy 数组是 N 维数组，我们可以用以下代码创建一个 3 维数组：

print(ones((2, 3, 4)))

在处理多维数据时， ... 运算符非常有用，它表示“所有剩余维度”。

a = ones((2, 3, 4))
print(a[0, ...])
print(a[0, 1, ...])

比较不同的数据存储方式

在存储数据时，我们可以选择使用 N 维数组或扁平的数据结构。以比较不同抵押贷款为例，固定抵押贷款支付取决于贷款金额、利率和还款次数三个参数。我们可以使用 arange() 函数来生成参数范围，并使用不同的数据结构来存储结果。

from numpy import *
loans = arange(100000, 160000, 20000)
num_payments = arange(5, 30, 5)*12
interests = arange(3, 5.5, 0.5)/100.0/12.0

# method 1, storing results in a list
res1 = []

# method 2, storing results in an array
res2 = zeros([len(loans), len(num_payments), len(interests)])

for i, loan in enumerate(loans):
    for j, num_pay in enumerate(num_payments):
        for k, interest in enumerate(interests):
            res1.append([loan, num_pay, interest, -pmt(interest, num_pay, loan)])
            res2[i][j][k] = -pmt(interest, num_pay, loan)

使用 N 维数组的好处是索引更快，而扁平列表的结果更易读。选择哪种数据结构取决于具体问题。

数学函数与傅里叶展开可视化

NumPy 数组支持简单的算术运算和大多数数学函数。我们可以用一个傅里叶展开的例子来展示数组作为向量进行操作的强大之处。

from pylab import *

# prepare the plot
figure()
hold(True)

# number of points to display the wave
N = 2**8
t = linspace(0, 1, N)
y = zeros(N)

for n in range(1, 8, 2):
    # the sine waves, added
    y += 4/(pi*n)*sin(2*pi*n*t*2)

    # plot the graph
    plot(t, y)

# annotate the graph
axis([0, 1, -1.4, 1.4])
grid()
xlabel('Time [seconds]')
ylabel('Value []')
title('Fourier expansion of a rectangular wave')
show()

在这个例子中，我们使用正弦波的傅里叶展开来生成矩形波。通过对数组进行操作，我们可以轻松地实现这个过程，并将结果可视化。

总之，Python 中的随机模拟和 NumPy 模块为我们解决各种数学问题和数据处理任务提供了强大的工具。通过合理选择数据结构和使用数学函数，我们可以更高效地完成任务。

Python 中的数学游戏与 NumPy 模块应用

数据存储方式的进一步探讨

在前面我们已经了解了使用 N 维数组和扁平列表来存储抵押贷款相关数据。接下来，我们深入探讨一下它们在不同场景下的特点以及如何根据具体问题进行选择。

索引效率对比

使用 N 维数组时，索引操作非常高效。例如，当我们固定贷款金额为某一值时，在 res2 数组中可以直接通过 res2[0, ...] 来获取相关结果。这个操作是直接根据数组的维度和索引进行定位的，时间复杂度较低。

exec(open('ndflat.py').read())
for row in res1:
    if(row[0] == 120000 and row[1] == 120):
        print(row)
res2[1, 1, ...]

而对于扁平列表 res1 ，要获取相同条件下的结果，需要遍历整个列表并进行条件判断。这意味着时间复杂度会随着列表长度的增加而线性增长，尤其是当数据量较大时，这种性能差异会更加明显。

可读性分析

扁平列表 res1 的结果更加直观和易于理解。它以一种人类可读的形式展示了所有参数的组合，每一行都清晰地列出了贷款金额、还款次数、利率和每月还款额等信息。

values = res2[1, 1, :]
for i, v in enumerate(values):
    row = [loans[1], num_payments[1], interests[i], v]
    print(row)

相比之下，N 维数组 res2 虽然在存储和索引上有优势，但要将其结果以可读的形式呈现出来，需要额外的处理，如使用 for 循环进行遍历和格式化。

数据存储与文件处理

扁平列表非常适合存储为 CSV 文件，因为 CSV 本身就是一种扁平的数据结构。我们可以直接将列表中的每一行写入 CSV 文件，操作简单方便。

而对于 N 维数组，虽然也可以将其扁平化后存储为 CSV 文件，但需要进行额外的处理。例如，使用列表推导式将 3 维数组 res2 扁平化：

exec(open('ndflat.py').read())
res2.shape
[n1, n2, n3] = [range(i) for i in res2.shape]
flat_res2 = [res2[a, b, c] for a in n1 for b in n2 for c in n3]

数学函数在数组上的应用

NumPy 数组支持丰富的数学函数，这些函数可以直接应用于数组的每个元素，极大地提高了计算效率。

简单算术运算

数组可以进行加法、减法、乘法、除法和指数运算等简单算术操作。例如：

from numpy import *
a = array([1, 2, 3])
b = array([4, 5, 6])
print(a + b)
print(a * b)

数学模块函数的应用

大多数 math 和 cmath 模块中的数学函数在 NumPy 中都有实现。例如，我们可以使用 sin() 、 cos() 等三角函数对数组进行操作：

t = linspace(0, 2*pi, 100)
y = sin(t)
plot(t, y)
show()

傅里叶展开的深入理解

在前面的例子中，我们展示了如何使用傅里叶展开来生成矩形波。下面我们进一步分析这个过程。

傅里叶展开原理

傅里叶展开是将一个周期函数表示为一系列正弦和余弦函数的和。在我们的例子中，使用正弦波的叠加来逼近矩形波。公式为 $y = \sum_{n=1,3,5,…}^{\infty} \frac{4}{\pi n} \sin(2\pi n t N)$，其中 $N$ 表示周期数。

代码实现分析

from pylab import *
figure()
hold(True)
N = 2**8
t = linspace(0, 1, N)
y = zeros(N)
for n in range(1, 8, 2):
    y += 4/(pi*n)*sin(2*pi*n*t*2)
    plot(t, y)
axis([0, 1, -1.4, 1.4])
grid()
xlabel('Time [seconds]')
ylabel('Value []')
title('Fourier expansion of a rectangular wave')
show()

在代码中，我们首先创建了时间向量 t ，然后初始化结果向量 y 。在 for 循环中，每次将一个正弦波添加到 y 中，并绘制当前的结果。随着循环的进行，正弦波的叠加越来越接近矩形波。

总结与建议

通过以上的讨论，我们可以总结出以下几点：

• 随机模拟 ：随机模块可以帮助我们解决各种概率问题，通过模拟大量事件并进行统计分析，得到问题的近似解。同时，可视化可以帮助我们更好地理解问题和验证结果。
• 随机序列操作 ：随机模块提供的序列操作函数可以方便地处理序列数据，如扑克牌的随机抽取和洗牌等。
• NumPy 数组 ：NumPy 的 ndarray 对象是强大的数据处理工具，支持多种数组创建、切片、索引和重塑操作。在处理多维数据时，N 维数组可以提高存储和索引效率，但在可读性方面可能需要额外处理。
• 数学函数 ：NumPy 支持丰富的数学函数，可以直接应用于数组，提高计算效率。
• 数据存储 ：选择 N 维数组还是扁平数据结构取决于具体问题。如果需要高效的索引和计算，N 维数组是更好的选择；如果更注重数据的可读性和与 CSV 文件的兼容性，扁平列表可能更合适。

在实际应用中，我们应该根据问题的特点和需求，合理选择数据结构和算法，以达到最佳的性能和效果。同时，不断学习和掌握 Python 中的各种工具和技术，将有助于我们更好地解决复杂的数学和数据处理问题。