如何实现积分 Chi^2 函数的反函数?
How to implement the inverse of the integrated Chi^2 function?
我在一篇名为 A Neural Bayesian Estimator for Conditional Probability Densities 的论文中实现了一些变量预处理。
它指出:
1.) 给定非线性,应用单调变量变换 F:t->s 使得 s 均匀分布。这可以通过以下方式实现,如论文中所述:
>>> sorting the target vector in ascending order
>>> fitting the spline to data, e.g. using interpolate from scipy
2.) 之后 s 被缩放到 -1 和 1 之间。这可以通过 interp
:
来实现
>>> from numpy import interp
>>> interp(256,[1,512],[5,10])
3.) 最后,平坦分布需要转换为高斯分布,以零为中心,std 1。
虽然前两个步骤很清楚如何实施,但我正在为第三个步骤而苦苦挣扎。
关于3.),作者进一步说明可以使用积分X^2(X...chi)函数的反函数。有没有图书馆,最好是 Python,适合这个工作?
更新 1:
又看了一遍论文,好像X^2和chi没有直接关系,而是这样计算的:
X^2 = P*(1-o)^2+(1-P)*((-1)-o)^2
P
作为纯度(给定一些变量可以很容易地计算出来)和 o
变量本身。
对于缩放到 -1 和 1 之间的给定 s,我可能只计算下限 = -1 和上限 = s 的积分,然后得到它的倒数。
问题:如何用数字表示?
如果您的意思是 X2 与描述的 PDF 分发 here, then what you're looking is X2 CDF. It is expressed via incomplete Gamma function, see same reference, and you could use SciPy to compute it, this or that 应该符合要求。不要忘记分母中的完整 Gamma 函数。
要找到不完全 Gamma 的反函数,您可以查看 SciPy 中的反函数:this or that。
因此,我认为您不需要所有这些插值内容
更新
可以分析地计算该表达式,例如,使用在线积分器
喜欢 that。只需计算上限的结果和下限的结果之间的差异,你就设置好了
更新二
你必须自己设置间隔
以下是(绝对未经测试!)您可以尝试使用的代码。注意,我使用通用
求根例程,虽然因为积分是多项式,但更优化的方法可能
是使用 here 的多项式根,甚至
自己编码 - 它只是一个 cubic equation
def intgrl(x):
return x*(x*(3.0 + x - 6.0*p) + 3.0)/3.0
def CDF(x, norm):
return (intgrl(x) - intgrl(-1.0))/norm
def f(x, norm, rn):
return CDF(x, norm) - rn
norm = intgrl(1.0) - intgrl(-1.0)
rn = 0.12345
res = scipy.optimize.brentq(f, -1.0, 1.0, args=(norm, rn))
更新三
变量 rn
被定义为从 0 到 1 的一些(随机 U(0,1))数。
from scipy.optimize import brentq
import numpy as np
import matplotlib.pyplot as plt
def denormPDF(x, p):
return p*(1.0-x)**2 + (1.0-p)*((-1.0)-x)**2
def intgrl(x, p):
return x*(x*(3.0 + x - 6.0*p) + 3.0)/3.0
def CDF(x, p, norm):
return (intgrl(x, p) - intgrl(-1.0, p))/norm
def PDF(x, p, norm):
return denormPDF(x, p)/norm
def f(x, p, norm, rn):
return CDF(x, p, norm) - rn
p = 0.25
norm = intgrl(1.0, p) - intgrl(-1.0, p)
x = np.linspace(-1.0, 1.0, 100)
y = [PDF(x, p, norm) for x in x]
z = [CDF(x, p, norm) for x in x]
# plot PDF
plt.plot(x, y)
plt.show()
# plot CDF
plt.plot(x, z)
plt.show()
rn = np.linspace(0.000001, 1.0-0.000001, 50)
iCDF = [brentq(f, -1.0, 1.0, args=(p, norm, rn)) for rn in rn]
# plot inverse CDF
plt.plot(rn, iCDF)
plt.show()
我在一篇名为 A Neural Bayesian Estimator for Conditional Probability Densities 的论文中实现了一些变量预处理。
它指出: 1.) 给定非线性,应用单调变量变换 F:t->s 使得 s 均匀分布。这可以通过以下方式实现,如论文中所述:
>>> sorting the target vector in ascending order
>>> fitting the spline to data, e.g. using interpolate from scipy
2.) 之后 s 被缩放到 -1 和 1 之间。这可以通过 interp
:
>>> from numpy import interp
>>> interp(256,[1,512],[5,10])
3.) 最后,平坦分布需要转换为高斯分布,以零为中心,std 1。
虽然前两个步骤很清楚如何实施,但我正在为第三个步骤而苦苦挣扎。
关于3.),作者进一步说明可以使用积分X^2(X...chi)函数的反函数。有没有图书馆,最好是 Python,适合这个工作?
更新 1:
又看了一遍论文,好像X^2和chi没有直接关系,而是这样计算的:
X^2 = P*(1-o)^2+(1-P)*((-1)-o)^2
P
作为纯度(给定一些变量可以很容易地计算出来)和 o
变量本身。
对于缩放到 -1 和 1 之间的给定 s,我可能只计算下限 = -1 和上限 = s 的积分,然后得到它的倒数。
问题:如何用数字表示?
如果您的意思是 X2 与描述的 PDF 分发 here, then what you're looking is X2 CDF. It is expressed via incomplete Gamma function, see same reference, and you could use SciPy to compute it, this or that 应该符合要求。不要忘记分母中的完整 Gamma 函数。
要找到不完全 Gamma 的反函数,您可以查看 SciPy 中的反函数:this or that。
因此,我认为您不需要所有这些插值内容
更新
可以分析地计算该表达式,例如,使用在线积分器 喜欢 that。只需计算上限的结果和下限的结果之间的差异,你就设置好了
更新二
你必须自己设置间隔
以下是(绝对未经测试!)您可以尝试使用的代码。注意,我使用通用 求根例程,虽然因为积分是多项式,但更优化的方法可能 是使用 here 的多项式根,甚至 自己编码 - 它只是一个 cubic equation
def intgrl(x):
return x*(x*(3.0 + x - 6.0*p) + 3.0)/3.0
def CDF(x, norm):
return (intgrl(x) - intgrl(-1.0))/norm
def f(x, norm, rn):
return CDF(x, norm) - rn
norm = intgrl(1.0) - intgrl(-1.0)
rn = 0.12345
res = scipy.optimize.brentq(f, -1.0, 1.0, args=(norm, rn))
更新三
变量 rn
被定义为从 0 到 1 的一些(随机 U(0,1))数。
from scipy.optimize import brentq
import numpy as np
import matplotlib.pyplot as plt
def denormPDF(x, p):
return p*(1.0-x)**2 + (1.0-p)*((-1.0)-x)**2
def intgrl(x, p):
return x*(x*(3.0 + x - 6.0*p) + 3.0)/3.0
def CDF(x, p, norm):
return (intgrl(x, p) - intgrl(-1.0, p))/norm
def PDF(x, p, norm):
return denormPDF(x, p)/norm
def f(x, p, norm, rn):
return CDF(x, p, norm) - rn
p = 0.25
norm = intgrl(1.0, p) - intgrl(-1.0, p)
x = np.linspace(-1.0, 1.0, 100)
y = [PDF(x, p, norm) for x in x]
z = [CDF(x, p, norm) for x in x]
# plot PDF
plt.plot(x, y)
plt.show()
# plot CDF
plt.plot(x, z)
plt.show()
rn = np.linspace(0.000001, 1.0-0.000001, 50)
iCDF = [brentq(f, -1.0, 1.0, args=(p, norm, rn)) for rn in rn]
# plot inverse CDF
plt.plot(rn, iCDF)
plt.show()