实现贝茨分布

Question

我一直在尝试绘制贝茨分布曲线，贝茨分布是n独立标准均匀变量（从 0 到 1）的平均值的分布。

（我在时间间隔上工作[-1;1]，我对变量进行了简单的更改）。

在 n 次之后，曲线不稳定，这阻止了我继续前进。为了考虑变量x是连续的，我在10**6个样本中采样了interval。以下是不同的一些示例n：

但是n大于29，曲线发散，并且越大n，发散引起的变形越接近曲线的（平均）中心：

Bates 概率分布定义如下：

我的代码：

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仅用于告诉 x 轴上的除法。一个小得多的数字就足够了。（在下面的代码中，我将变量重命名为xaxis_steps）。
• 使用appendforF会非常慢。最好创建一个正确大小的 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图。此类曲线的平滑度取决于样本的数量。kde 可以通过seaborn 的 kdeplot直接绘制。您还可以单独计算 kde 函数，然后将其应用于给定的 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()