обратная сторона cdf

Question

Я хотел бы вычислить обратную интегральную функцию плотности (обратный cdf) данного pdf.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);

, что означает, что я выполняю выборку с линейной интерполяциейПереверните и уменьшите гистограммуЭто довольно уродливый код, не очень надежный (если я изменю коэффициент повышения частоты, я могу получить действительно разные результаты), и он бесполезно медленен ... кто-нибудь может предложить лучший подход?

Примечание:Обобщенное обратное Я пытаюсь вычислить (в случае, если PDF не является обратимым):

F^{-1}(t) = \inf{x \in R ; F(x)>t }   

с F функция накопленной плотности

[РЕДАКТИРОВАТЬ: на самом деле, K = 1 (т.е.., без повышения частоты дискретизации), кажется, дает более точные результаты ...]

Спасибо!

Answer 1

Если ваш ввод указан в виде ненормализованной гистограммы, то просто с помощью встроенной функции quantile() автоматически вычисляется точка данных для указанного квантиля, что и делает обратный CDF. Если гистограмма нормализована по количеству точек данных (что делает ее вектором вероятности), то сначала просто умножьте ее на количество точек данных. См. здесь для подробностей quantile(). По сути, вы предполагаете, что с учетом вашей гистограммы / данных первый параметр является фиксированным, что превращает quantiles() в функцию только с указанными значениями вероятности p. Вы можете легко написать функцию-обертку, чтобы сделать ее более удобной в случае необходимости. Это избавляет от необходимости явно вычислять CDF с cumsum().

Добавлена ​​

Если мы предположим, что гистограмма, ячейки и количество точек данных равны 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() для этого, если нежелательно сначала интерполировать.

Answer 2

Хорошо, я думаю, что нашел гораздо более короткую версию, которая работает по крайней мере так же быстро и точно:

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

[РЕДАКТИРОВАТЬ: Нет, на самом деле это все еще в два-три раза медленнее, чем первоначальный код ... не спрашивайте меня, почему! И он не обрабатывает не строго монотонные функции: это приводит к ошибке: «Значения X должны быть различны»]