Comment savoir si mes données correspondent à la distribution Pareto?

10

J'ai un échantillon qui est un vecteur avec 220 nombres. Voici un lien vers un histogramme de mes données. . Et je souhaite vérifier si mes données correspondent à une distribution de Pareto, mais je ne veux pas voir les tracés QQ avec cette distribution, mais j'ai besoin d'une réponse exacte avec une valeur de p dans R, comme le test de normalité d'Anderson-Darling ( ad.test) . Comment pourrais-je faire ça? Soyez aussi précis que possible.

sévère
la source
1
Le résultat d'un test statistique ne vous dira pas que vos données ont une distribution de Pareto . En fait, vous pouvez être certain que s'il s'agit de données réelles, elles n'ont pas de distribution de Pareto. Tout un test vous montrera si vous disposez de suffisamment de données pour détecter l'écart de quantité de Pareto que vous avez. Autrement dit, s'il rejette tout ce qu'il dit est «oui, la taille de l'échantillon était assez grande pour vous dire ce que vous saviez déjà». Pourquoi voudriez-vous entreprendre un tel exercice, un exercice qui ne peut pas répondre à la vraie question que vous avez?
Glen_b -Reinstate Monica
Votre question n'est-elle rien de plus que «quelles lignes de code I dois-je écrire pour que le programme R exécute la procédure X»? Ensuite, c'est hors sujet ici. Cela pourrait être considéré comme une question de programmation. S'il y a un aspect statistique à votre question (comme «est-ce que cela a du sens?»), Vous devez clarifier et mettre l'accent sur ces aspects
Glen_b -Reinstate Monica
1
Passons maintenant au test d'Anderson-Darling (ou, d'ailleurs, au Kolmogorov-Smirnov que @Zen a suggéré ci-dessus). Ce sont des tests pour des distributions complètement spécifiées . Autrement dit, pour que les tests aient les propriétés souhaitées, vous devez spécifier a priori ( PAS estimer ) tous les paramètres. Vous ne pouvez donc pas utiliser l'un ou l'autre pour cet exercice car vous n'avez pas de paramètres prédéfinis. (Vraisemblablement, vous faites cela à la suggestion de quelqu'un d'autre. Il est très difficile d'expliquer les idées fausses à quelqu'un via un intermédiaire.)
Glen_b -Reinstate Monica
Pourquoi faites-vous ces tests? Par exemple, quelles actions changeront selon que vous refusez ou refusez de rejeter?
Glen_b -Reinstate Monica
Vous devriez toujours regarder un tracé QQ, quel que soit votre motif. Et vous ne devez pas fétichiser une valeur P "exacte". Un test différent vous donnerait une valeur P "exacte" différente.
Nick Cox

Réponses:

13

(PS) Tout d'abord, je pense que Glen_b a raison dans ses commentaires ci-dessus sur l'utilité d'un tel test: les données réelles ne sont certainement pas exactement distribuées par Pareto, et pour la plupart des applications pratiques, la question serait "quelle est la qualité de l'approximation de Pareto?" - et le tracé QQ est un bon moyen de montrer la qualité d'une telle approximation.

De toute façon, vous pouvez faire votre test avec la statistique de Kolmogorov-Smirnov, après avoir estimé les paramètres par maximum de vraisemblance. Cela empêche d'estimation des paramètres à utiliser le -value de , de sorte que vous pouvez faire bootstrap paramétrique pour estimer. Comme Glen_b le dit dans le commentaire, cela peut être connecté au test de Lilliefors .pks.test

Voici quelques lignes de code R.

Définissez d'abord les fonctions de base pour gérer les distributions de Pareto.

# distribution, cdf, quantile and random functions for Pareto distributions
dpareto <- function(x, xm, alpha) ifelse(x > xm , alpha*xm**alpha/(x**(alpha+1)), 0)
ppareto <- function(q, xm, alpha) ifelse(q > xm , 1 - (xm/q)**alpha, 0 )
qpareto <- function(p, xm, alpha) ifelse(p < 0 | p > 1, NaN, xm*(1-p)**(-1/alpha))
rpareto <- function(n, xm, alpha) qpareto(runif(n), xm, alpha)

La fonction suivante calcule le MLE des paramètres (justifications dans Wikipedia ).

pareto.mle <- function(x)
{
  xm <- min(x)
  alpha <- length(x)/(sum(log(x))-length(x)*log(xm))
  return( list(xm = xm, alpha = alpha))
}

Et ces fonctions calculent la statistique KS et utilisent le bootstrap paramétrique pour estimer la valeur .p

pareto.test <- function(x, B = 1e3)
{
  a <- pareto.mle(x)

  # KS statistic
  D <- ks.test(x, function(q) ppareto(q, a$xm, a$alpha))$statistic

  # estimating p value with parametric bootstrap
  B <- 1e5
  n <- length(x)
  emp.D <- numeric(B)
  for(b in 1:B)
  {
    xx <- rpareto(n, a$xm, a$alpha);
    aa <- pareto.mle(xx)
    emp.D[b] <- ks.test(xx, function(q) ppareto(q, aa$xm, aa$alpha))$statistic
  }

  return(list(xm = a$xm, alpha = a$alpha, D = D, p = sum(emp.D > D)/B))
}

Maintenant, par exemple, un échantillon provenant d'une distribution de Pareto:

> # generating 100 values from Pareto distribution
> x <- rpareto(100, 0.5, 2)
> pareto.test(x)
$xm
[1] 0.5007593

$alpha
[1] 2.080203

$D
         D 
0.06020594 

$p
[1] 0.69787

... et à partir d'un :χ2(2)

> # generating 100 values from chi square distribution
> x <- rchisq(100, df=2)
> pareto.test(x)
$xm
[1] 0.01015107

$alpha
[1] 0.2116619

$D
        D 
0.4002694 

$p
[1] 0

Notez que je ne prétends pas que ce test est non biaisé: lorsque l'échantillon est petit, un biais peut exister. Le bootstrap paramétrique ne prend pas bien en compte l'incertitude sur l'estimation des paramètres (pensez à ce qui se passerait lors de l'utilisation de cette stratégie pour tester naïvement si la moyenne d'une variable normale avec une variance inconnue est nulle).

PS Wikipedia en dit quelques mots. Voici deux autres questions pour lesquelles une stratégie similaire a été suggérée: test d'adéquation pour un mélange , test d'adéquation pour une distribution gamma .

Elvis
la source
3
Lorsque vous ajustez la distribution de la statistique de test pour l'estimation des paramètres de cette manière, ce n'est pas un test KS (même s'il est basé sur des statistiques KS) - c'est un type particulier de test Lilliefors . Ce n'est plus non paramétrique, mais on peut en construire un par simulation pour une distribution donnée. Lilliefors l'a fait spécifiquement pour la normale et l'exponentielle ... dans les années 1960.
Glen_b -Reinstate Monica
Merci pour ce commentaire @Glen_b je ne le savais pas.
Elvis
Aucun problème; cela ne change rien au contenu de ce que vous faites (ce qui est bien comme ça), seulement ce qu'il faut appeler.
Glen_b -Reinstate Monica
@Glen_b J'ai apporté des changements substantiels dans ma réponse, merci encore!
Elvis