实现贝茨分布

Question

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

（我在区间[-1;1]上工作，我对变量做了一个简单的更改）

这条曲线在n之后会不稳定，这会阻止我前进。 为了考虑变量x是连续的，在10×6样本中对采样间隔进行采样。以下是不同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

主要问题是，具有交替和的公式极易出现数值精度问题

避免右侧问题的一个技巧是假设分布是对称的，只计算一半

一个简单的精度优化是通过调用scipy.special.comb来替换combinaison公式中的阶乘。这避免了需要划分非常大的数字

较小精度的优化是同时计算奇偶数的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）
• 对F使用append将非常缓慢。最好创建一个大小正确的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's kdeplot直接绘制。您还可以单独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()