numpy - 反正切函数的数值计算

Question

我很难解释反正切函数的结果。这种行为对于我遇到的所有实现都是一致的，所以我将在这里限制自己使用 NumPy 和 MATLAB。

这个想法是有一个随机放置的点圈。目标是表示它们在极坐标系中的位置，并且由于它们是均匀分布的，我希望 θ 角（使用atan2函数计算）也可以在区间 -π ... π 上随机分布。

这是MATLAB的代码：

stp = 2*pi/2^8;
siz = 100;
num = 100000000;

x = randi([-siz, siz], [1, num]);
y = randi([-siz, siz], [1, num]);
m = (x.^2+y.^2) < siz^2;

[t, ~] = cart2pol(x(m), y(m));
figure()
histogram(t, -pi:stp:pi);

对于 Python 和 NumPy：

import numpy as np
import matplotlib.pyplot as pl

siz = 100
num = 100000000

rng = np.random.default_rng()
x = rng.integers(low=-siz, high=siz, size=num, endpoint=True)
y = rng.integers(low=-siz, high=siz, size=num, endpoint=True)
m = (x**2+y**2) < siz**2

t = np.arctan2(y[m], x[m]);
pl.hist(t, range=[-np.pi, np.pi], bins=2**8)
pl.show()

在这两种情况下，我得到的结果都是这样的，可以很容易地看到 π/4 的每个倍数的“步数”。

它看起来像是某种精度误差，但奇怪的是我没想到的角度。这种行为也存在于普通atan功能中。

Answer 1

请注意，您使用的是整数

因此，对于每一对 (p,q)，您将拥有floor(sqrt(p**2 + q**2)/gcd(p,q)/r)给出相同角度的对arctan(p,q)。那么对于 (p,q) 的倍数gcd(p,q)是1

另请注意，p**2+q**2对于1的倍数pi/2和2奇数倍数pi/4，我们可以预测，偶数倍数的项目将pi/4多于奇数倍数的项目pi/4。这与我们在你的情节中看到的一致。

例子

让我们用整数坐标绘制位于半径为 10 的圆中的点。

import numpy as np
import matplotlib.pyplot as plt
from collections import Counter
def gcd(a,b):
    if a == 0 or b == 0:
        return max(a,b)
    while b != 0:
        a,b = b, a%b
    return a;
R = 10
x,y = np.indices((R+1, R+1))
m = (x**2 + y**2) <= R**2
x,y = x[m], y[m]
t = np.linspace(0, np.pi / 2)
plt.figure(figsize=(6, 6))
plt.plot(x, y, 'o')
plt.plot(R * np.cos(t), R * np.sin(t))

lines = Counter((xi / gcd(xi,yi), 
                 yi / gcd(xi,yi)) for xi, yi in zip(x,y))
plt.axis('off')
for (x,y),f in lines.items():
    if f != 1:
        r = np.sqrt(x**2 + y**2)
        plt.plot([0, R*x/r], [0, R*y/r], alpha=0.25)
        plt.text(R*1.03*x/r, R*1.03*y/r, f'{int(y)}/{int(x)}: {f}')

在这里，您可以在图中看到一些与其他点具有相同角度的点。45 度有 7 个点，90 的倍数有 10 个。许多点都有一个独特的角度。基本上，您有很多角度，但很少有点，还有一些角度可以达到很多点。

但总的来说，这些点在角度方面几乎是均匀分布的。在这里，我绘制了几乎是一条直线的累积频率（如果分布是均匀的，那将是什么），并且 bin 频率形成了一些三角形分形图案。

R = 20
x,y = np.indices((R+1, R+1))
m = (x**2 + y**2) <= R**2
x,y = x[m], y[m]
plt.figure(figsize=(6,6))
plt.subplot(211)
plt.plot(np.sort(np.arctan2(x,y))*180/np.pi, np.arange(len(x)), '.', markersize=1)

plt.subplot(212)
plt.plot(np.arctan2(x,y)*180/np.pi, np.gcd(x,y), '.', markersize=4)

如果圆圈的大小增加并且您使用足够宽的 bin 制作直方图，您将不会注意到变化，否则您将在直方图中看到这种模式。