cdf的逆

Question

我想计算给定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(即没有上采样)似乎能给出更准确的结果……

谢谢!

Answer 1

如果您的输入是以非标准化直方图的形式指定的，那么只需使用内置的quantile()函数就可以自动计算指定分位数的数据点，这就是逆向CDF所做的事情。如果直方图按数据点的数量进行归一化(使其成为概率向量)，则首先将其乘以数据点的数量。有关quantile()的详细信息，请参阅here。基本上，您将假设给定直方图/数据，第一个参数是固定的，它将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，如果不希望首先进行插值，则在该CDF上使用interp1()。

Answer 2

好的，我想我找到了一个更短的版本，它至少同样快速和准确：

cdf = cumsum(pdf);
M = length(cdf);
xx = linspace(0,1,M);
invcdf = interp1(cdf,xx,xx)

编辑:不，这实际上仍然比初始代码慢两到三倍……别问我为什么！并且它不能处理非严格单调的函数:这会产生错误：“X的值应该是不同的”