柯西分布采样时如何修正误差

Question

我正在尝试使用所谓的（我认为）inverse transform sampling 从柯西分布中抽样值。我的程序如下：

from numpy import*
from matplotlib.pyplot import*
from scipy.interpolate import interp1d

def f(x):
    return 1 / (pi*(1+x**2)) # Cauchy

x = linspace(-6,+6,100)
dx = x[1]-x[0]
cdf = cumsum(f(x))*dx

plt.plot(x,cdf)
plt.xlabel('X')
plt.ylabel('CDF')
plt.show()

from_cdf_to_x = interp1d(cdf, x)
cdf_new = np.random.uniform(0,1,100)
x_sampled = from_cdf_to_x(cdf_new)

情节看起来不错，但我在第二部分收到以下我不理解的错误：A value in x_new is above the interpolation range。就我而言，这意味着什么，纠正它的正确方法是什么？无论我怎么搜索，我都没有发现我编写的程序有什么不正确的地方...

Answer 1

你的数组 cdf_new 的范围是从 0 到 1.0，但是柯西累积分布渐近接近 1 而没有达到它，因此插值没有覆盖均匀分布的整个范围。生成 Cauchy 累积分布时，您需要减小 cdf_new and/or 的顶部范围

增加 x 的范围

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

def f(x):
    return 1 / (np.pi*(1+x**2)) # Cauchy

x = np.linspace(-100,+100,1000)
dx = x[1]-x[0]
cdf = np.cumsum(f(x))*dx

plt.plot(x,cdf)
plt.xlabel('X')
plt.ylabel('CDF')
plt.show()

from_cdf_to_x = interp1d(cdf, x)
cdf_new = np.random.uniform(0.05,0.95,10000)
x_sampled = from_cdf_to_x(cdf_new)
plt.hist(x_sampled, bins=20)
plt.show()

这是您所期望的正态分布：

柯西分布采样时如何修正误差

How to correct error when sampling Cauchy distribution

python

numpy

transform

sampling

inverse