7

I'm attempting to create an ECDF (and a confidence bound) from data in Python. I can generate the ECDF fairly easily with numpy by sorting and using linspace. However, I'm not entirely certain what the appropriate confidence bounds are, and there don't seem to be any built-in libraries that calculate the bounds (statsmodels seems to just give the ECDF).

If I want a point-wise confidence bound of $1-\alpha$ is it appropriate to use the DKW inequality to calculate my region with

$$C_n(\alpha) = \sqrt{\frac{1}{2n}\log\left(\frac{2}{\alpha}\right)} \,,$$

where $n$ is the number of observations in my sample? Thus if $F(x)$ is my ECDF, my upper and lower bounds would be

$$\mathrm{UB}(x) = \min\left(1, F(x)+C_n(\alpha)\right)$$ $$\mathrm{LB}(x) = \max\left(0, F(x)-C_n(\alpha)\right)$$

MATLAB has a built-in function ECDF, but I didn't have much luck understanding how to apply Greenwood's Formula (referenced at the bottom) to generate the bounds.

Jim
  • 2,152
ALollz
  • 205
  • 1
    Have a look at https://stats.stackexchange.com/questions/298290/plotting-non-parametric-ecdef-confidence-envelopes-for-comparison and my comments there. Look at the code of ecdf.ksCI in CRAN package sfsmisc. That code is simple (I wrote it ...) so should be easy to translate to python ... – kjetil b halvorsen Jun 28 '18 at 08:47
  • 1
    What is not entirely clear is that you ask for point-wise bounds, but suggest to use an inequality, DKW, that is mostly employed to find simultaneous confidence bounds. Please note: Greenwood's approach is point-wise. – Jim Jul 04 '18 at 09:10
  • @Jim So would applying that method above just be incorrect, or would it instead give me the simultaneous confidence bound? – ALollz Jul 04 '18 at 19:24
  • 1
    It would give you the simultaneous one. Which is fine if that's what you're after. But note that a simultaneous CI is wider than a point-wise one. – Jim Jul 04 '18 at 19:37

1 Answers1

3

In Matlab's console type:

edit ecdf

It opens the source code in the editor.

Go to line 194:

if nargout>2 || (nargout==0 && isequal(bounds,'on'))

This is the start of the code block that calculates the lower - and upper (confidence) bounds: [Flo, Fup]. The code block is 30 lines long and pretty straightforward. Posted below for your convenience:

if nargout>2 || (nargout==0 && isequal(bounds,'on'))
     % Get standard error of requested function
     if cdf_sf % 'cdf' or 'survivor'
         se = NaN(size(D));
         if N(end)==D(end)
            t = 1:length(N)-1;
         else
            t = 1:length(N);
         end
         se(t) = S(t) .* sqrt(cumsum(D(t) ./ (N(t) .* (N(t)-D(t))))); % <--- line 203
     else % 'cumhazard'
         se = sqrt(cumsum(D ./ (N .* N)));
     end

     % Get confidence limits
     zalpha = -norminv(alpha/2);
     halfwidth = zalpha*se;
     Flo = max(0, Func - halfwidth);
     Flo(isnan(halfwidth)) = NaN; % max drops NaNs, put them back
     if cdf_sf % 'cdf' or 'survivor'
         Fup = min(1, Func + halfwidth);
         Fup(isnan(halfwidth)) = NaN; % max drops NaNs
     else % 'cumhazard'
         Fup = Func + halfwidth; % no restriction on upper limit
     end
         Flo = [NaN; Flo];
         Fup = [NaN; Fup];
else 
     Flo = [];
     Fup = [];
end

The square root of Greenwood's formula, i.e.

$$ S(t) \sqrt{\sum_{t_i < T} \frac{d_i}{r_i(r_i - d_i)}} \,, $$

is implemented in line 203 as:

se(t) = S(t) .* sqrt(cumsum(D(t) ./ (N(t) .* (N(t)-D(t)))));

Can you take it from here? Let me know.

Jim
  • 2,152