2
$\begingroup$

I have a model for data in my experiment that states that the data has a Gumbel distribution with known location and scale.

I am then looking at observations with very high scores that I suspect to be outliers.

Is there a Bayesian approach (or otherwise) to get the probability that this observation is from the underlying Gumbel distribution?

As an end result I would like to decide whether or not I accept this observation as coming from the underlying Gumbel distribution.

Edit

I am in fact working with biological data and have demonstrated both empirically and through a set of assumptions that the distribution is in fact of Gumbel type with parameters that can be estimated from the data itself and/or from a theoretical perspective given some assumptions.

Now, while most data behaves in-line with this distribution, I would like to find the probability that observations with very large values might not be from the same Gumbel distribution but due to some other phenomenon that is not described in this distribution that I want to detect.

In this case, updating the data using the new observation would not be useful (especially considering that this observation might come from another random variable with a different distribution).

I am able to calculate the one-side probability of how extreme the observation is by using.

$$1-F(x)$$

From my (very limited) understanding of statistics and probability, I fear that quoting this value might be erroneous. My main concern is that this is the "probability of observing more extreme observations" and NOT the "probability that this observation is not from the distribution".

Any thoughts on this?

$\endgroup$
  • $\begingroup$ How could you know, that your data comes from Gumbel distribution with known parameters? Could you clarify this issue? If this is a prior distribution, that you (or some other expert) have validated with other data, then the Bayesian approach could be applied to update those parameters with new data. For example if you know a priori that previous studies have confirmed that distribution is Gumbel with parameters $\mu =1, \beta =2$, and you have a sample on new data $d$, then you can find posterior distribution for Gumbel parameters, i.e. update them with new information from new sample. $\endgroup$ – Tomas Apr 2 '12 at 20:27
  • $\begingroup$ the probability $1-F(x)$ is a sensible quantification of how extreme the observation is. Exactly what do you mean by the "probability that this observation is not from the distribution"? Is it the (Bayesian) conditional probability of the Gumbel distribution given the observation, in which case you need a prior on a set of distributions to begin with? If you have a single alternative to the Gumbel in mind a you can compute the likelihood ratio. You might also consider the possibility of a mixture model with one Gumbel component. $\endgroup$ – NRH Apr 3 '12 at 10:57
  • $\begingroup$ Yes, I am in fact looking for the (Bayesian) conditional probability of the Gumbel distribution given the observation. In this case, I might set $P(A)$ to, say 0.999, the prior probability of the observation coming from a Gumbel distribution, while $P(B) = 1 - P(A)$ is the evidence. $P(B|A) = 1 - F(x)$ is the probability of the observation given the Gumbel distribution. Bayes' theorem can now allow us to find $P(A|B)$. Does this make sense? $\endgroup$ – Andrew Apr 3 '12 at 12:45
  • $\begingroup$ In light of the clarifying edit, I have retagged this question with the outliers tag and edited the title. I also changed the term "random variable" in the original--which has a well understood conventional meaning--to "observation," which is in line with the apparent intent of the question. $\endgroup$ – whuber Apr 3 '12 at 15:53
  • $\begingroup$ @Andrew, you cannot compute the probability sought from the prior and the Gumbel alone. You need to specify the alternative distribution as well for this computation. You also need the density and not distribution function. I have edited my answer to reflect this question. $\endgroup$ – NRH Apr 4 '12 at 8:03
2
$\begingroup$

If you have a known distribution, Gumbel or not, with distribution function $F$ and an observation $x$ a simple one-sided probability of how extreme the observation is simply the probability $$1-F(x)$$ of an observation larger than the actual observation $x$.

If $F$ is given in terms of a scale-location transformation of a distribution with distribution function $G$ with scale parameter $\sigma$ and location parameter $\mu$ then the probability is $$1 - F(x) = 1 - G\left( \frac{x - \mu}{\sigma}\right).$$

For the Gumbel distribution with $G(x) = \exp(-\exp(-x))$ this amounts to computing $$1 - \exp\left(-\exp\left(-\frac{x-\mu}{\sigma}\right)\right).$$ Note that in this parametrization of a location-scale transformation of $G$ the location parameter $\mu$ and the scale parameter $\sigma$ are not the mean and standard deviation of the distribution.

With a single, known distribution the probability could have a Bayesian or a frequency interpretation. If the parameters $\mu$ and $\sigma$ are estimated based on a data set a Bayesian approach could be to integrate the formula above over the posterior on $(\mu, \sigma)$ whereas a typical plug-in frequency approach is to simply plug in the estimates of $\mu$ and $\sigma$. In the latter case, a confidence interval on the estimated tail probability is preferred to be able to judge the uncertainty in the estimated probability.

In some applications, e.g. in biological sequence analysis, the probability is used to judge a very large number of observations, and though the computation is the same, we do face a multiple testing problem, which should be accounted for in the evaluation of the computed probabilities.

Edit

The following was added to answer a question in a comment by the OP. If there are two known distributions, or models, with densities $f_1$ and $f_2$, respectively, and we assign prior probabilities $\pi_1$ and $\pi_2$ to these the conditional probability of model 1 given $x$ is $$\frac{\pi_1 f_1(x)}{\pi_1 f_1(x) + \pi_2 f_2(x)}.$$ This easily generalizes to a finite number of models, but we have to specify these other models in addition to the prior probabilities to carry out the computation of the posterior probability of model 1.

$\endgroup$
5
$\begingroup$

Edit

The following material responds to an earlier version of the question, which I interpreted as asking how to test whether a collection of observations is consistent with an assumed Gumbel distribution with known parameters. Edits to the question now suggest the interest lies in identifying individual (high) outlying observations, which is a slightly different question. The entire collection of data could depart from a Gumbel shape, for instance, without containing any outliers. However, the presence of an outlier or outlier also constitutes a departure from the assumed distribution, so the methods described here remain relevant to the problem, albeit indirectly.


Usually this kind of question is answered inversely: you test the hypothesis that the data came from the specified distribution.

The Kolmogorov-Smirnov test tends to work well: so well, in fact, that with largish datasets (numbering in the hundreds or greater) it can detect such small deviations that it can be too powerful in some applications. But here it sounds like such power is desired.

The ks.test routine in R can be persuaded to conduct this test even though R does not appear to support Gumbel distributions. Use the fact that the exponential of a Gumbel variate is a location and scale version of the Exponential distribution (a Gamma distribution with shape 1 and scale 1).

For example, suppose x is an array of numbers and the supposed Gumbel parameters are mu for location and beta for scale. Then

x <- mu - beta * log(-log(runif(n)))
ks.test(exp(-(x - mu)/beta), "pgamma", 1, 1)$p.value

yields a p-value for this test. As a check, we can generate n iid values from a Gumbel distribution with parameters mu.actual and beta.actual and test them against supposed values mu and beta:

trial <- function(n, mu, beta, mu.actual, beta.actual) {
    x <- mu.actual - beta.actual * log(-log(runif(n)))
    ks.test(exp(-(x - mu)/beta), "pgamma", 1, 1)$p.value
}

When the parameter values agree, we should obtain a uniform distribution of p-values:

data <- replicate(10000, trial(8, 0, 1, 0, 1))
hist(data)

Null distribution of p-values

It looks good. To assess the power of this test, we draw random variables from a slightly different distribution and look again:

data <- replicate(10000, trial(8, 0, 1, 1, 1))
hist(data)

Alternative distribution of p-values

In this case, using (say) a test level of 0.05 and drawing only 8 values from a Gumbel(1,1) distribution and supposing it is a Gumbel(0,1) distribution--a shift of +1--this plot shows we would have about a 60% chance of detecting the difference, because about 60% of the results are to the left of 0.05. That chance increases to 90+% for 16 draws (not shown here). Detecting a change in shape is harder--but you can run whatever simulations are appropriate for your hypotheses to learn the details.

$\endgroup$
  • $\begingroup$ as always, a detailed and accurate answer. However, the OP could actually mean exactly what he asks, that you have a known Gumbel distribution and want to judge if a single observation is extreme in this distribution. This is a simpler question than testing the Gumbel assumption, which is encountered in, for instance, biological sequence analysis when using local alignment tools like BLAST. $\endgroup$ – NRH Apr 3 '12 at 7:57
  • $\begingroup$ Thank you, NRH. I believe this reply did answer "exactly what he asks," but in light of recent edits and comments, it is clear I did not pick up on the ambiguous language and that you read the OP's mind better than I did :-). $\endgroup$ – whuber Apr 3 '12 at 16:01
  • $\begingroup$ I do apologise for my lack of clarity when trying to speak "statistics". Your answer is however much appreciated and it would be interesting (and possibly useful for other readers) if we could create another question that requires the answer above. $\endgroup$ – Andrew Apr 3 '12 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.