找到图表的轮廓

Question

我正在尝试使用 matplotlib.pyplot 在 python 中找到图形轮廓的方程（以 μ 和 σ 表示）。我已经用 numpy 和 scipy 概述了该图，但我找不到一种方法来获取两条线的均值和标准差方程。

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull


x_axis = []
y_axis = []

nth_triple = []
average_triple = []

def return_average(lst):
    counter = 0
    for el in lst:
        counter += el
    return counter/len(lst)

def pythagoreanTriplets(limits) : 
    counter = 0
    c, m = 0, 2
    while c < limits :  
        for n in range(1, m) : 
            a = m * m - n * n 
            b = 2 * m * n 
            c = m * m + n * n 
            if c > limits : 
                break

            triple = [a, b, c]

            x_axis.append(a)
            y_axis.append(b)

            nth_triple.append(counter)
            average_triple.append(return_average(triple))
            counter += 1

        m = m + 1

limit = 100000
pythagoreanTriplets(limit)

points = np.column_stack((nth_triple, average_triple))
hull = ConvexHull(points)

mu_x = np.mean(points[hull.vertices, 0])
mu_y = np.mean(points[hull.vertices, 1])
sigma_x = np.std(points[hull.vertices, 0])
sigma_y = np.std(points[hull.vertices, 1])


plt.scatter(nth_triple, average_triple, s=1, marker=',', color='orange', edgecolors='none')
plt.plot(points[hull.vertices, 0], points[hull.vertices, 1], color='black', label='Convex Hull Outline', linestyle='dotted')

plt.show()

虚线是我想要找到的方程：

任何指导将不胜感激。预先感谢。

Answer 1

标准偏差与此事无关；您已经绘制了平均值，因此我们将分析平均值。

你没有提到这一点，但你正在使用欧几里得公式。这是最简单的部分。更棘手的部分是你有一个非单调的 n。您需要做一些代数来计算出 m 的连续函数作为计数器索引的函数；那么我们就有了

import math

import matplotlib
import matplotlib.pyplot as plt
import statistics

import numpy as np

'''
a**2 + b**2 == c**2, all integers;
a = m*m - n*n
b = 2*m*n
c = m*m + n*n

(a + b + c) = m*m - n*n + 2*m*n + m*m + n*n
            = 2*m*m + 2*m*n
            = 2*m*(m + n)
            
(m-1)(m-2)/2 == counter
m2 - 3m + 2 - 2counter = 0
m = (3 + sqrt(9 - 4*(3 - 2counter)))/2
m = 0.5*(3 + math.sqrt(9 - 8*(1 - i))
Minimum when n = 1, (a + b + c)/3 = 2*m*(m + 1)/3
Maximum when n = m - 1, (a + b + c)/3 = 2*m*(m + m - 1)/3
                                       = 2*m*(2*m - 1)/3
'''


def pythagorean_triplets(limits: int) -> tuple[
    list[int],
    list[int],
]:
    nth_triple = []
    average_triple = []

    counter = 0
    c, m = 0, 2

    while c < limits:
        for n in range(1, m):
            # Euclid's formula
            a = m*m - n*n
            b = 2*m*n
            c = m*m + n*n
            if c > limits:
                break

            triple = (a, b, c)
            nth_triple.append(counter)
            mean = statistics.mean(triple)
            average_triple.append(mean)
            counter += 1
        m += 1

    return nth_triple, average_triple


def main() -> None:
    nth_triple, average_triple = pythagorean_triplets(limits=10_000)

    matplotlib.use('TkAgg')
    plt.scatter(
        nth_triple, average_triple,
        s=3, marker=',', color='orange', edgecolors='none',
    )
    m = 0.5*(3 + np.sqrt(9 - 8*(1 - np.array(nth_triple))))
    plt.plot(nth_triple, 2*m*(m + 1)/3)
    plt.plot(nth_triple, 2*m*((m - 1)*2 - 1)/3)
    plt.show()


if __name__ == '__main__':
    main()