逆 CDF 变换采样的分布略有错误
Slighty wrong distribution from inverse CDF transform sampling
背景
我需要从具有已知累积分布函数 (CDF) 的相当复杂的概率密度函数 (PDF) 中随机抽样，我正在尝试使用 inverse transform sampling。这应该很容易做到，因为我有 CDF，只需要在插入统一随机数的同时对其进行数值反转（不可能用代数方式进行）。但是，结果分布的方差低于预期，我在 CDF 中找不到任何错误。
所以我通过从正态分布中抽样来简化和测试我的算法。结果是一样的：位置没问题，但尺度不对。我知道有更好的内置高斯采样方法，但这只是对采样算法的测试。
这个问题最初出现在 Fortran 中，但我在 Python 中复制了这个问题，所以我必须做一些根本性的错误或者有数值问题。
Python
import numpy as np
from scipy.special import erf
from scipy.optimize import brentq
import matplotlib.pyplot as plt
from scipy.stats import norm
def testfunc(x):
## Test case, result should be 6.04880103
# out = 0.5 * (1. + erf((x - 5.) / (2. * np.sqrt(2.)))) - 0.7
r = np.random.uniform()
# hand-built cdf:
# out = 0.5 * (1. + erf((x - 5.) / (2. * np.sqrt(2.)))) - r
# scipy cdf:
out = norm.cdf(x, 5, 2) - r
return out
if __name__ == '__main__':
n = 10000
sol_array = np.zeros(n)
for i in range(0, n):
sol_array[i] = brentq(testfunc, -100.,100.)
print('mean = ' + str(np.mean(sol_array)))
print('std = ' + str(np.std(sol_array)))
plt.hist(sol_array, normed=True, bins='fd')
x = np.linspace(-1, 11, 1000)
plt.plot(x, norm.pdf(x, 5, 2))
plt.show()
正如预期的那样，采样值的平均值约为 5，但对于我手工构建的 CDF 和 scipy 的 CDF，标准差约为 1.28，应该为 2。
这在直方图中也可见：
Fortran
Fortran 中的相同问题，尽管生成的标准差值不同。代码较长，因为包含求解器。该求解器是 Alan Miller of an old FORTRAN 77 netlib-routine (zeroin.f) 翻译的 Fortran 90 版本。
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: n = 1000000
real, dimension(n) :: v
real :: mean, std
integer, dimension(:), allocatable :: seed
integer :: i, seedsize, clock
! seed the PRNG
call random_seed(size=seedsize)
allocate(seed(seedsize))
call system_clock(count=clock)
seed=clock + 37 * (/ (i - 1, i=1, seedsize) /)
call random_seed(put=seed)
deallocate(seed)
do i = 1, n
v(i) = real(zeroin(testfunc, -100._dp, 100._dp, 1e-20_dp, 1e-10_dp))
end do
mean = sum(v) / n
std = sum((v - mean)**2) / n
print*, mean, std
contains
function testfunc(v)
implicit none
real(dp), intent(in) :: v
real(dp) :: testfunc, r
call random_number(r)
! testfunc = 0.5 * (1. + erf((v-5.)/(2.*sqrt(2.)))) - 0.7 ! should be 6.04880
testfunc = 0.5 * (1. + erf((v-5.)/(2.*sqrt(2.)))) - r ! Gaussian test with mu=5 and sigma=2
end function testfunc
function zeroin(f, ax, bx, aerr, rerr) result(fn_val)
! original zeroin.f from netlib.org
! code converted using to_f90 by alan miller
! date: 2003-07-14 time: 12:32:54
!-----------------------------------------------------------------------
! finding a zero of the function f(x) in the interval (ax,bx)
! ------------------------
! INPUT:
! f function subprogram which evaluates f(x) for any x in the
! closed interval (ax,bx). it is assumed that f is continuous,
! and that f(ax) and f(bx) have different signs.
! ax left endpoint of the interval
! bx right endpoint of the interval
! aerr the absolute error tolerance to be satisfied
! rerr the relative error tolerance to be satisfied
! OUTPUT:
! abcissa approximating a zero of f in the interval (ax,bx)
!-----------------------------------------------------------------------
! zeroin is a slightly modified translation of the algol procedure
! zero given by Richard Brent in "Algorithms for Minimization without
! Derivatives", Prentice-Hall, Inc. (1973).
implicit none
real(dp), intent(in) :: ax
real(dp), intent(in) :: bx
real(dp), intent(in) :: aerr
real(dp), intent(in) :: rerr
real(dp) :: fn_val
real(dp) :: a, b, c, d, e, eps, fa, fb, fc, tol, xm, p, q, r, s, atol, rtol
interface
real(selected_real_kind(15, 307)) function f(x)
real(selected_real_kind(15, 307)), intent(in) :: x
end function f
end interface
! compute eps, the relative machine precision
eps = epsilon(0.0_dp)
! initialization
a = ax
b = bx
fa = f(a)
fb = f(b)
if (fb*fa > 0.) then
print*, 'a, b, fa, fb', a, b, fa, fb
stop
end if
atol = 0.5 * aerr
rtol = max(0.5_dp*rerr, 2.0_dp*eps)
! begin step
10 c = a
fc = fa
d = b - a
e = d
20 if (abs(fc) < abs(fb)) then
a = b
b = c
c = a
fa = fb
fb = fc
fc = fa
end if
! convergence test
tol = rtol * max(abs(b),abs(c)) + atol
xm = 0.5 * (c-b)
if (abs(xm) > tol) then
if (fb /= 0.0) then
! is bisection necessary
if (abs(e) >= tol) then
if (abs(fa) > abs(fb)) then
! is quadratic interpolation possible
if (a == c) then
! linear interpolation
s = fb / fc
p = (c-b) * s
q = 1.0 - s
else
! inverse quadratic interpolation
q = fa / fc
r = fb / fc
s = fb / fa
p = s * ((c-b)*q*(q-r)-(b-a)*(r-1.0))
q = (q-1.0) * (r-1.0) * (s-1.0)
end if
! adjust signs
if (p > 0.0) q = -q
p = abs(p)
! is interpolation acceptable
if (2.0*p < (3.0*xm*q-abs(tol*q))) then
if (p < abs(0.5*e*q)) then
e = d
d = p / q
go to 30
end if
end if
end if
end if
! bisection
d = xm
e = d
! complete step
30 a = b
fa = fb
if (abs(d) > tol) b = b + d
if (abs(d) <= tol) b = b + sign(tol,xm)
fb = f(b)
if (fb*(fc/abs(fc)) > 0.0) go to 10
go to 20
end if
end if
! done
fn_val = b
end function zeroin
end
所得样本的平均值约为 5，而标准差约为 1.64。
问题
有谁知道我的算法在什么地方会出现数值问题？ Python 版本和 Fortran 版本都有同样的问题，但程度不同，这让我相信这是浮点数的一些舍入问题，但我想不起来在哪里。即使求解器 returns 是一个四舍五入的值，该差异也不应显示在简单的直方图中。
有没有人发现我的算法有错误？我理解错了吗？
我只检查了Python版本，确实是执行错误
即，您的 testfunc，即 root-finding brentq 例程的目标函数，行为 non-deterministically。在 root-finding 运行 期间（即一次调用 brentq() 直到它 returns），brentq 需要多次调用相同的回调直到达到收敛。但是，每次 brentq 调用此回调时， 目标方程都会发生变化 ，因为 r 会获得一个新的 pseudo-random 值。因此，root-finding 例程无法收敛到您想要的解决方案。
从概念上讲，您需要做的是首先生成均匀随机变量的样本，然后应用相同、确定性变换（即逆分布函数）到他们。当然你不需要做root-solving，因为你可以使用scipy.stats随机变量[=44的ppf（百分位函数，即cdf的反函数）方法=].
作为概念证明，您可以运行以下代码对标准均匀样本数组使用（不必要的昂贵且不太准确）转换方法：
import numpy
import numpy.random
from scipy.optimize import brentq
from scipy.stats import norm
# Setup
n = 10000
numpy.random.seed(0x5eed)
ran_array = numpy.random.uniform(size=n)
sol_array = numpy.empty_like(ran_array)
# Target function for root-finding: F(x) - p = 0 where p is the probability level
# for which the quantile is to be found
def targetfunc(x, p):
return norm.cdf(x, 5, 2) - p
for i in range(n):
sol_array[i] = brentq(targetfunc, -100.0, 100.0, args=(ran_array[i],))
print("mean = %10f" % sol_array.mean())
print("std = %10f" % sol_array.std())
输出：
mean = 5.011041
std = 2.009365
背景
我需要从具有已知累积分布函数 (CDF) 的相当复杂的概率密度函数 (PDF) 中随机抽样，我正在尝试使用 inverse transform sampling。这应该很容易做到，因为我有 CDF，只需要在插入统一随机数的同时对其进行数值反转（不可能用代数方式进行）。但是，结果分布的方差低于预期，我在 CDF 中找不到任何错误。
所以我通过从正态分布中抽样来简化和测试我的算法。结果是一样的：位置没问题，但尺度不对。我知道有更好的内置高斯采样方法，但这只是对采样算法的测试。
这个问题最初出现在 Fortran 中，但我在 Python 中复制了这个问题，所以我必须做一些根本性的错误或者有数值问题。
Python
import numpy as np
from scipy.special import erf
from scipy.optimize import brentq
import matplotlib.pyplot as plt
from scipy.stats import norm
def testfunc(x):
## Test case, result should be 6.04880103
# out = 0.5 * (1. + erf((x - 5.) / (2. * np.sqrt(2.)))) - 0.7
r = np.random.uniform()
# hand-built cdf:
# out = 0.5 * (1. + erf((x - 5.) / (2. * np.sqrt(2.)))) - r
# scipy cdf:
out = norm.cdf(x, 5, 2) - r
return out
if __name__ == '__main__':
n = 10000
sol_array = np.zeros(n)
for i in range(0, n):
sol_array[i] = brentq(testfunc, -100.,100.)
print('mean = ' + str(np.mean(sol_array)))
print('std = ' + str(np.std(sol_array)))
plt.hist(sol_array, normed=True, bins='fd')
x = np.linspace(-1, 11, 1000)
plt.plot(x, norm.pdf(x, 5, 2))
plt.show()
正如预期的那样，采样值的平均值约为 5，但对于我手工构建的 CDF 和 scipy 的 CDF，标准差约为 1.28，应该为 2。
这在直方图中也可见：
Fortran
Fortran 中的相同问题，尽管生成的标准差值不同。代码较长，因为包含求解器。该求解器是 Alan Miller of an old FORTRAN 77 netlib-routine (zeroin.f) 翻译的 Fortran 90 版本。
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: n = 1000000
real, dimension(n) :: v
real :: mean, std
integer, dimension(:), allocatable :: seed
integer :: i, seedsize, clock
! seed the PRNG
call random_seed(size=seedsize)
allocate(seed(seedsize))
call system_clock(count=clock)
seed=clock + 37 * (/ (i - 1, i=1, seedsize) /)
call random_seed(put=seed)
deallocate(seed)
do i = 1, n
v(i) = real(zeroin(testfunc, -100._dp, 100._dp, 1e-20_dp, 1e-10_dp))
end do
mean = sum(v) / n
std = sum((v - mean)**2) / n
print*, mean, std
contains
function testfunc(v)
implicit none
real(dp), intent(in) :: v
real(dp) :: testfunc, r
call random_number(r)
! testfunc = 0.5 * (1. + erf((v-5.)/(2.*sqrt(2.)))) - 0.7 ! should be 6.04880
testfunc = 0.5 * (1. + erf((v-5.)/(2.*sqrt(2.)))) - r ! Gaussian test with mu=5 and sigma=2
end function testfunc
function zeroin(f, ax, bx, aerr, rerr) result(fn_val)
! original zeroin.f from netlib.org
! code converted using to_f90 by alan miller
! date: 2003-07-14 time: 12:32:54
!-----------------------------------------------------------------------
! finding a zero of the function f(x) in the interval (ax,bx)
! ------------------------
! INPUT:
! f function subprogram which evaluates f(x) for any x in the
! closed interval (ax,bx). it is assumed that f is continuous,
! and that f(ax) and f(bx) have different signs.
! ax left endpoint of the interval
! bx right endpoint of the interval
! aerr the absolute error tolerance to be satisfied
! rerr the relative error tolerance to be satisfied
! OUTPUT:
! abcissa approximating a zero of f in the interval (ax,bx)
!-----------------------------------------------------------------------
! zeroin is a slightly modified translation of the algol procedure
! zero given by Richard Brent in "Algorithms for Minimization without
! Derivatives", Prentice-Hall, Inc. (1973).
implicit none
real(dp), intent(in) :: ax
real(dp), intent(in) :: bx
real(dp), intent(in) :: aerr
real(dp), intent(in) :: rerr
real(dp) :: fn_val
real(dp) :: a, b, c, d, e, eps, fa, fb, fc, tol, xm, p, q, r, s, atol, rtol
interface
real(selected_real_kind(15, 307)) function f(x)
real(selected_real_kind(15, 307)), intent(in) :: x
end function f
end interface
! compute eps, the relative machine precision
eps = epsilon(0.0_dp)
! initialization
a = ax
b = bx
fa = f(a)
fb = f(b)
if (fb*fa > 0.) then
print*, 'a, b, fa, fb', a, b, fa, fb
stop
end if
atol = 0.5 * aerr
rtol = max(0.5_dp*rerr, 2.0_dp*eps)
! begin step
10 c = a
fc = fa
d = b - a
e = d
20 if (abs(fc) < abs(fb)) then
a = b
b = c
c = a
fa = fb
fb = fc
fc = fa
end if
! convergence test
tol = rtol * max(abs(b),abs(c)) + atol
xm = 0.5 * (c-b)
if (abs(xm) > tol) then
if (fb /= 0.0) then
! is bisection necessary
if (abs(e) >= tol) then
if (abs(fa) > abs(fb)) then
! is quadratic interpolation possible
if (a == c) then
! linear interpolation
s = fb / fc
p = (c-b) * s
q = 1.0 - s
else
! inverse quadratic interpolation
q = fa / fc
r = fb / fc
s = fb / fa
p = s * ((c-b)*q*(q-r)-(b-a)*(r-1.0))
q = (q-1.0) * (r-1.0) * (s-1.0)
end if
! adjust signs
if (p > 0.0) q = -q
p = abs(p)
! is interpolation acceptable
if (2.0*p < (3.0*xm*q-abs(tol*q))) then
if (p < abs(0.5*e*q)) then
e = d
d = p / q
go to 30
end if
end if
end if
end if
! bisection
d = xm
e = d
! complete step
30 a = b
fa = fb
if (abs(d) > tol) b = b + d
if (abs(d) <= tol) b = b + sign(tol,xm)
fb = f(b)
if (fb*(fc/abs(fc)) > 0.0) go to 10
go to 20
end if
end if
! done
fn_val = b
end function zeroin
end
所得样本的平均值约为 5，而标准差约为 1.64。
问题
有谁知道我的算法在什么地方会出现数值问题？ Python 版本和 Fortran 版本都有同样的问题，但程度不同，这让我相信这是浮点数的一些舍入问题，但我想不起来在哪里。即使求解器 returns 是一个四舍五入的值，该差异也不应显示在简单的直方图中。
有没有人发现我的算法有错误？我理解错了吗？
我只检查了Python版本，确实是执行错误
即，您的 testfunc，即 root-finding brentq 例程的目标函数，行为 non-deterministically。在 root-finding 运行 期间（即一次调用 brentq() 直到它 returns），brentq 需要多次调用相同的回调直到达到收敛。但是，每次 brentq 调用此回调时， 目标方程都会发生变化 ，因为 r 会获得一个新的 pseudo-random 值。因此，root-finding 例程无法收敛到您想要的解决方案。
从概念上讲，您需要做的是首先生成均匀随机变量的样本，然后应用相同、确定性变换（即逆分布函数）到他们。当然你不需要做root-solving，因为你可以使用scipy.stats随机变量[=44的ppf（百分位函数，即cdf的反函数）方法=].
作为概念证明，您可以运行以下代码对标准均匀样本数组使用（不必要的昂贵且不太准确）转换方法：
import numpy
import numpy.random
from scipy.optimize import brentq
from scipy.stats import norm
# Setup
n = 10000
numpy.random.seed(0x5eed)
ran_array = numpy.random.uniform(size=n)
sol_array = numpy.empty_like(ran_array)
# Target function for root-finding: F(x) - p = 0 where p is the probability level
# for which the quantile is to be found
def targetfunc(x, p):
return norm.cdf(x, 5, 2) - p
for i in range(n):
sol_array[i] = brentq(targetfunc, -100.0, 100.0, args=(ran_array[i],))
print("mean = %10f" % sol_array.mean())
print("std = %10f" % sol_array.std())
输出：
mean = 5.011041 std = 2.009365