python - 实现贝茨分布

Question

我一直在尝试绘制贝茨分布曲线，贝茨分布是 n 的均值分布独立的标准统一变量(从 0 到 1)。
(我在区间 [-1;1] 上工作，我对变量进行了简单的更改)。
在 n 次之后，曲线不稳定，这阻止了我继续前进。
为了考虑变量x是连续的，我在10**6个样本中采样了interval。以下是不同 n 的一些示例:

但是对于 n大于 29，曲线发散，更大的 n ，发散引起的变形越接近曲线的(平均)中心:

贝茨概率分布定义如下:

我的代码:

samples=10**6

def combinaison(n,k):   # combination of K out of N
  cnk=fac(n)/(fac(k)*fac(abs(n-k))) # fac is factoriel 
  return cnk


def dens_probas(a,b,n):
  x=np.linspace(a, b, num=samples)
  y=(x-a)/(b-a)
  F=list()
  for i in range(0,len(y)):
    g=0
    for k in range(0,int(n*y[i]+1)):
      g=g+pow(-1,k)*combinaison(n,k)*pow(y[i]-k/n,n-1)
    d=(n**n/fac(n-1))*g
    F.append(d)         
  return F 

任何想法来纠正更大的分歧 n ?

Answer 1

主要问题是具有交替总和的公式极易出现数值精度问题。
避免右侧问题的一个技巧是假设分布是对称的并且只计算它的一半。
一个直接的精度优化是替换公式中的阶乘 combinaison调用 scipy.special.comb .这避免了需要划分非常大的数字。
一个较小的精度优化是计算 g偶数和奇数放在一起。但是乍一看公式不能减少太多，所以替换:

        for k in range(0, int(floor(n * y[i] + 1))):
            g += pow(-1, k) * combinaison(n, k) * pow(y[i] - k / n, n - 1)

经过:

        last_k = int(floor(n * y[i]))
        for k in range(0, last_k + 1, 2): # note that k increments in steps of 2
            if k == last_k:
                g += combinaison(n, k) * (pow(y[i] - k / n, n - 1))
            else:
                g += combinaison(n, k) * (pow(y[i] - k / n, n - 1) - pow(y[i] - (k + 1)/ n, n - 1) * (n - k) / (k + 1))

其他一些注意事项:

变量 samples仅用于告诉 xaxis 中的除法。一个小得多的数字就足够了。 (在下面的代码中，我将变量重命名为 xaxis_steps )。

使用 append为 F将非常缓慢。最好创建一个大小正确的 numpy 数组，然后将其填充。(这也使复制一半更容易。)

from matplotlib import pyplot as plt
import numpy as np
from scipy.special import comb
from math import factorial as fac
from math import floor

xaxis_steps = 500

def combinaison(n, k):  # combination of K out of N
    return comb(n, k)

def dens_probas(a, b, n):
    x = np.linspace(a, b, num=xaxis_steps)
    y = (x - a) / (b - a)
    F = np.zeros_like(y)
    for i in range(0, (len(y)+1) // 2):
        g = 0
        for k in range(0, int(floor(n * y[i] + 1))):
            g += pow(-1, k) * combinaison(n, k) * pow(y[i] - k / n, n - 1)
        F[i] = (n ** n / fac(n - 1)) * g
        F[-i-1] = F[i]  # symmetric graph
    plt.plot(x, F, label=f'n={n}')
    return F

for n in (5, 30, 50, 80, 90):
    dens_probas(-1, 1, n)
plt.legend()
plt.show()

所有这些优化共同解决了 n=30 的准确性问题。到附近 n=80 :

一种完全不同的方法是生成大量统一的样本并采取手段。从这些样本 kde可以生成图。此类曲线的平滑度取决于样本数量。可以通过 seaborn's kdeplot 直接绘制 kde .您也可以单独 calculate the kde function ，然后将其应用于给定的 x 范围并通过标准 matplotlib 绘制它。

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kde

num_samples = 10 ** 5

def dens_probas(a, b, n):
    samples = np.random.uniform(a, b, size=(num_samples, n)).mean(axis=1)
    samples = np.hstack([samples, a + b - samples])  # force symmetry; this is not strictly necessary
    return gaussian_kde(samples)

for n in (5, 30, 50, 80, 90, 200):
    kde = dens_probas(-1, 1, n)
    xs = np.linspace(-1, 1, 1000)
    F = kde(xs)
    plt.plot(xs, F, label=f'n={n}')
plt.legend()
plt.show()