Численная оценка функций арктангенса

Question

У меня возникают трудности с интерпретацией результатов функций арктангенса. Это поведение одинаково для всех реализаций, с которыми я сталкивался, поэтому здесь я ограничусь NumPy и MATLAB.

Идея состоит в том, чтобы иметь круг случайно расположенных точек. Цель состоит в том, чтобы представить их положения в полярной системе координат, и, поскольку они равномерно распределены, я ожидаю, что угол θ (который рассчитывается с использованием функция) также будет случайным образом распределена по интервалу -π ... π.

Вот код для 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.

Это похоже на какую-то ошибку точности, но странно для углов, где я этого не ожидал. Также такое поведение присутствует для обычных также функционировать.

Answer 1

Обратите внимание, что вы используете целые числа

Итак, для каждой пары (p, q) у вас будет floor(sqrt(p**2 + q**2)/gcd(p,q)/r)пары, дающие один и тот же угол arctan(p,q). Тогда для кратных (p,q) gcd(p,q)является

Заметьте также, что p**2+q**2является 1для кратных pi/2а также 2для нечетных кратных , с этим мы можем предсказать, что будет больше элементов, которые являются четными кратными, чем нечетные кратные 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. Многие из точек имеют уникальный угол. По сути, у вас есть много углов с несколькими точками и несколько углов, которые затрагивают много точек.

Но в целом точки распределены почти равномерно по углу. Здесь я строю кумулятивную частоту, которая представляет собой почти прямую линию (какой она была бы, если бы распределение было равномерным), а частота бина формирует некоторый треугольный фрактальный узор.

      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)

Если размер круга увеличивается, и вы делаете гистограмму с достаточно широкими интервалами, вы не заметите изменений, иначе вы увидите этот шаблон на гистограмме.