我想计算给定pdf的逆累积密度函数(逆cdf)。pdf直接作为直方图给出,即N个等距分量的向量。
我目前的做法是:
cdf = cumsum(pdf);
K = 3; %// some upsampling factor
maxVal = 1; %// just for my own usage - a scaling factor
M = length(cdf);
N = M*K; %// increase resolution for higher accuracy
y = zeros(N, 1);
cursor = 2;
for i=1:N
desiredF = (i-1)/(N-1)*maxVal;
while (cursor<M && cdf(cursor)<desiredF)
cursor = cursor+1;
end;
if (cdf(cursor)==cdf(cursor-1))
y(i) = cursor-1;
else
alpha = min(1, max(0,(desiredF - cdf(cursor-1))/(cdf(cursor)-cdf(cursor-1))));
y(i) = ((cursor-1)*(1-alpha) + alpha*cursor )/maxVal;
end;
end;
y = resample(y, 1, K, 0);
这意味着我用线性插值进行上采样,对直方图进行逆采样和下采样。这是一个相当丑陋的代码,不是很健壮(如果我改变上采样因子,我可以得到非常不同的结果),而且速度非常慢。。。有人能提出更好的方法吗?
注:我试图计算的广义逆(在cdf不可逆的情况下)是:
F^{-1}(t) = inf{x in R ; F(x)>t }
其中F为累积密度函数
[编辑:实际上,K=1(即,没有上采样)似乎可以给出更准确的结果…]
谢谢!
如果您的输入是以非标准化直方图的形式指定的,那么只需使用内置的quantile()
函数即可自动计算指定分位数的数据点,这就是反向CDF的作用。如果直方图由数据点的数量归一化(使其成为概率向量),那么只需先将其乘以数据点的数目。有关quantile()
的详细信息,请参见此处。基本上,您将假设给定直方图/数据,第一个参数是固定的,这将quantiles()
变成仅指定概率值p
的函数。如果需要的话,您可以很容易地编写一个包装器函数,使其更加方便。这消除了使用cumsum()
显式计算CDF的需要。
添加
如果我们假设直方图、仓和数据点的数量分别为h, b, and N
,那么:
h1 = N*h; %// Only if histogram frequencies have been normalized.
data = [];
for kk = 1:length(h1)
data = [data repmat(b(kk), 1, h1(kk))];
end
%// Set p to the probability you want the inv-cdf for...
p = 0.5;
inv_cdf = quantiles(data,p)
添加
对于必须利用现有PDF矢量的解决方案,我们可以执行以下操作。假设x_old
和pdf_old
分别是直方图仓和直方图频率。
p = 0.5; %// the inv-cdf probability that I want
num_points_i_want = 100; %// the number of points I want in my histogram vector
x_new = linspace(min(x_old),max(x_old),num_points_i_want);
pdf_new = interp1(x_old,pdf_old,x_new);
cdf_new = cumsum(pdf_new);
inv_cdf = min(x_new(cdf_new >= p));
或者,如果不希望首先进行插值,我们可以先创建cumsum()
CDF,然后在此基础上使用interp1()
。
好吧,我想我找到了一个更短的版本,它至少能同样快速准确地工作:
cdf = cumsum(pdf);
M = length(cdf);
xx = linspace(0,1,M);
invcdf = interp1(cdf,xx,xx)
[EDIT:不,这实际上仍然比初始代码慢两到三倍…不要问我为什么!而且它不能处理非严格单调的函数:这会产生错误:"X的值应该不同"]