# My data is not normally distributed: what can I do to estimate a tail probability?

Continuing on from my earlier question, I'm attempting to analyse the data qualitatively.

In the following plot, I make $10000$ samples where I count "the number of clashes". I plot $n$ vs. the number of times $n$ clashes occurred.

(The number of clashes is a measure of "how wrong" an attempted attack was [on the secret sharing scheme I'm looking at]).

(Drawn using tikzDevice for R, then edited manually.)

In R, it fails the shapiro.test, so it's not normally distributed:

> shapiro.test(z[1:5000])

Shapiro-Wilk normality test

data:  z[1:5000]
W = 0.9947, p-value = 1.597e-12


So:

Q: How can I estimate the probability $p$ of $0$ clashes from the above distribution?

It should be very small, around $10^{-14}$:

• I have a theoretical lower bound of $1.046 \times 10^{-14}$, and I expect it to be close to the actual value.
• I have made $10^{11}$ samples, and all had at least one clash.

I attempted to fit an exponential curve to the left hand side (drawn above): the curve is $$3.29 \times 10^{-12} \exp(0.56n)$$ which, when $n=0$ gives the estimate $\hat{p}=3.29 \times 10^{-16}$. But I know that this estimate is off by around a factor of $100$, which makes me think this is not the best approach. (Or maybe I should fit some other curve, or use more samples. Or maybe this level of confidence is to be expected.)

• I'm trying to show that $\mathrm{Pr}[0 \text{ clashes}]$ is small (say less than $10^{-8}$ or $10^{-9}$). So the estimate doesn't need to be precise, but I need to have confidence in the estimate.

• The theoretical maximum number of clashes is $220$ (this number can be achieved).

• "Do you know the statistical power of the Shapiro-Wilk test for such a large sample size?" In short, no, I don't. But we can compare the results to random data from a normal distribution:

> shapiro.test(rnorm(5000, mean = mean(z), sd = sd(z)))

Shapiro-Wilk normality test

data:  rnorm(5000, mean = mean(z), sd = sd(z))
W = 0.9996, p-value = 0.4053


While the results fluctuate between runs, they don't seem comparable to my data.

I also tried with fewer samples included and it didn't seem to "help".

> shapiro.test(z[1:100])

Shapiro-Wilk normality test

data:  z[1:100]
W = 0.9757, p-value = 0.06116


compared to

> shapiro.test(rnorm(100, mean = mean(z), sd = sd(z)))

Shapiro-Wilk normality test

data:  rnorm(100, mean = mean(z), sd = sd(z))
W = 0.9845, p-value = 0.2932


(Here, it fluctuates quite a lot.)

• I'm capable of making around $10^{10}$ samples, if it would help.

• This question might be a better fit for stats.stackexchange.com . – Rahul Jul 12 '14 at 2:16
• I'm currently having technical difficulties logging into stats.SE. Let's give it a chance here; if it doesn't get answered I'll vote to migrate. (Oh, and before anyone comments, I know "data" is plural, but I couldn't bare to write "My data are not normally distributed".) – Rebecca J. Stones Jul 12 '14 at 2:24
• If you have fat tails probably you don't want to use normall or exponential. Pateto? Modified Beta? i dont know but actuarial people probably would have good suggestions. Other good reason to try the stats site. – Sergio Parreiras Jul 12 '14 at 2:37
• There are a number of approaches better suited to tail probabilities than fitting an exponential curve. However, before I can offer suggestions, could you provide some more info? In particular, what level of precision do you need? What is the theoretical max number of clashes? Do you know the statistical power of the Shapiro-Wilk test for such a large sample size? (5,000 data points is a pretty big sample, so even very small deviations from normality will get ridiculously low p-values). – user76844 Jul 15 '14 at 3:53
• I updated my question. Hopefully this is useful. – Rebecca J. Stones Jul 15 '14 at 6:21

Let $C=$ Number of Clashes and $p_0 = P(C=0)$. In general, your distribution is a discrete probability mass function. However, since you only care about a single value $(C=0)$ the rest of the distribution is somewhat irrelevant.

Therefore, one simple approach is to develop a conservative confidence interval for $p_0$ by inverting the hypothesis test for a binomial proportion:

$H_0: p_0=p\;\;\text{vs.}\;\; H_a: p_0<p$

There are two common (approximate) approaches to constructing binomial proportion CIs : The normal approximation and the Wilson Score Interval. Unfortunately, since we expect $p_0 \ll .001$, these approximate intervals will likely not have the desired true confidence level.

Instead, I would recommend using a Clopper-Pearson confidence interval. Not only is it based on the exact binomial distribution, but it is guaranteed to achieve its nominal confidence level (as opposed to only approximately for the other two). In fact, it is usually quite conservative, in that a 95% CI may actually achieve a 99% coverage probability, but it will not be less than 95%. For "typical" values of p, (e.g., between .01 and .99) the approximate intervals are generally preferred because they are usually close to their nominal coverage, while the Clopper-Pearson interval is too conservative. However, this is not your situation, so the Clopper-Pearson is a safer interval for bounding purposes.

Now, with $10^{11}$ data points where each has at least one clash, your CI will be quite easy to construct, and rather close to 0, as you expected. Below is the $1-\alpha$ one-sided Clopper-Pearson CI for $p_0$: (Note that the LB=$1.046 \times 10^{-14}$, as per your theoretical calculation), assuming you have observed an uninterrupted string of N non-zero results.

Upper Bound = $\{p:(1-p)^{N}= \alpha\}=\{p:N\ln(1-p)= \ln(\alpha)\} \rightarrow p=1-\exp(\frac{\ln(\alpha)}{N})$.

Now, with $N=10^{11}$ you will have numerical issues for the usual values of $\alpha =\{.05,.01\}$However, since you indicated that you wanted to know if it is less than, say, $10^{-8}$, we can reverse this process and find the level of confidence with which we can say that $p_0\leq 10^{-8}$:

$10^{-8} = 1-\exp(\frac{\ln(\alpha)}{10^{11}}) \rightarrow \frac{\ln(\alpha)}{10^{11}} = \ln(1-10^{-8}) \rightarrow \alpha = (1-10^{-8})^{10^{11}}\ll 10^{-10}$. Therefore, you can be virtually certain that $p_0<10^{-8}$

If you need actual numbers, I went onto WolframAlpha and got the following values for quantities that you may be interested in:

Confidence level for $p_0<10^{-8}$ = $1-5 \times 10^{-435} \approx 1$

99% one sided CI $\approx 4.6 \times 10^{-11}$

95% one sided CI $\approx 3\times 10^{-11}$

Confidence level for $p_0 < 10^{-13}$ $\approx 0.01$

The last result may be somewhat disappointing, in that we can only bound $p_0$ within 3 orders of magnitude of the theoretical minimum with any reasonable degree of confidence, even with your enormous sample size. Such is the difficulty of extremely small probabilities.

Hope this helps.

• Thanks for taking the time to write this answer. It looks like just what I'm after; I'll look into it. – Rebecca J. Stones Jul 15 '14 at 23:45