逆CDF变换采样的分布略有错误
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了逆CDF变换采样的分布略有错误相关的知识，希望对你有一定的参考价值。
Background
我需要从具有已知累积分布函数（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，但标准偏差约为1.28，对于我手工制作的CDF和scipy的CDF，它应该是2。这在直方图中也可见：
Fortran
Fortran中的问题相同，但结果标准差的值不同。代码更长，因为包含了解算器。该解算器是由旧的FORTRAN 77 netlib例程（Alan Miller）的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。
Question
有没有人知道我的算法在哪些方面可能会成为数字问题？事实上，Python版本和Fortran版本都有相同的问题，但是在不同程度上让我相信它是浮点数问题的一些舍入，但我无法想象在哪里。即使求解器返回舍入值，该差异也不应显示在简单的直方图中。
有没有人在我的算法中看到错误？我理解错了吗？
答案
我只检查了Python版本，实际上确实存在错误。
也就是说，你的testfunc，即找根brentq例程的目标函数，表现非确定性。在根查找运行期间（即一次调用brentq()直到它返回），brentq需要多次调用相同的回调，直到达到收敛。但是，每当brentq调用此回调时，目标方程随着r获得新的伪随机值而变化。因此，根查找例程无法收敛到您想要的解决方案。
在概念上，您需要做的是首先生成均匀随机变量的样本，并对它们应用相同的确定性变换（即逆分布函数）。当然你不需要做根求解，因为你可以使用ppf随机变量类的cdf（百分位函数，即scipy.stats的逆）方法。
作为概念验证，您可以在标准统一样本数组上使用（不必要的昂贵且不太精确）转换方法运行以下代码：
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变换采样的分布略有错误的主要内容，如果未能解决你的问题，请参考以下文章
谁知道利用逆变换法(反变换法），用matlab编程正态分布随机变量随机数？
压缩感知——沃尔什-哈达玛(WHT)变换与逆变换的Matlab代码实现