Intervalles de tolérance non paramétriques pour les variables discrètes

8

Supposons que plusieurs personnes évaluent à quel point elles ont aimé un film sur une échelle discrète de 1 à 10, et que vous souhaitez un intervalle [ l , u ] tel qu'avec (au moins) 95% de confiance, (au moins) 90 % de toutes les personnes qui voient le film le noteront pas plus bas que l et pas plus haut que u . [ l , u ] est alors un intervalle de tolérance (bilatéral) avec une confiance de 95% et une couverture de 90%. (Pour être clair, une confiance de 95% implique que si vous répétiez cette procédure plusieurs fois, 95% des intervalles produits auraient une couverture de la population d'au moins 90%.) Bien sûr, nous voulons généralement que [ l , u ] soit aussi étroit que possible tout en répondant à nos exigences.

J'ai vu diverses méthodes non paramétriques pour construire des intervalles de tolérance pour des variables aléatoires continues. J'ai également vu des méthodes pour construire des intervalles de tolérance pour les variables binomiales et de Poisson. (Le package R toleranceimplémente plusieurs de ces méthodes; Young, 2010.) Mais qu'en est-il des variables discrètes lorsque la distribution est inconnue? C'est généralement le cas pour des échelles de notation comme celle de mon exemple, et en supposant qu'une distribution binomiale ne semble pas sûre car les données réelles de l'échelle de notation présentent souvent des étrangetés telles que la multimodalité.

Serait-il sensé de se rabattre sur les méthodes non paramétriques pour les variables continues? Sinon, qu'en est-il d'une méthode Monte Carlo telle que la génération de 1000 répliques bootstrap de l'échantillon et la recherche d'un intervalle qui capture au moins 90% de l'échantillon dans au moins 950 des répliques?

Young, DS (2010). tolérance: Un package R pour estimer les intervalles de tolérance. Journal of Statistical Software, 36 (5), 1–39. Récupéré de http://www.jstatsoft.org/v36/i05

Kodiologue
la source
voulez-vous dire binomial ou multinomial? multinomial permettrait un comportement multimodal?
seanv507
Je veux dire binôme. Dans le cas d'une échelle de notation, par exemple, vous définiriez le nombre d'essais de Bernoulli sur le nombre de points d'échelle. Les intervalles entre les catégories d'une distribution multinomiale n'auraient pas beaucoup de sens, je pense, car les catégories ne sont pas ordonnées.
Kodiologist
@Kodiologist votre variable de résultat est une « échelle discrète de 1 à 10 » mais cela signifie qu'il est une réponse ordonnée multinomial. (Ou est-ce que je ne reçois pas quelque chose?)
Jim
@Jim "Multinomial ordonné" est un peu un oxymore. Dans une distribution multinomiale, l'ordre des catégories est arbitraire.
Kodiologist

Réponses:

1

La variable d'intérêt est distribuée multinomialement avec des probabilités de classe (cellule): . De plus, les classes sont dotées d'un ordre naturel.p1,p2,...,p10

Première tentative: le plus petit "intervalle prédictif" contenant90%

p     = [p1, ..., p10] # empirical proportions summing to 1
l     = 1
u     = length(p)
cover = 0.9

pmass = sum(p)

while (pmass - p[l] >= cover) OR (pmass - p[u] >= cover)
    if p[l] <= p[u]
       pmass = pmass - p[l]
       l     = l + 1  
    else # p[l] > p[u]
       pmass = pmass - p[u]
       u     = u - 1
    end        
end

Une mesure non paramétrique de l'incertitude (par exemple, variance, confiance) dans les estimations quantiles pourrait en effet être obtenue par des méthodes de bootstrap standard .l,u

Deuxième approche: "recherche bootstrap" directe

Ci-dessous, je fournis du code Matlab exécutable qui aborde la question directement du point de vue du bootstrap (le code n'est pas vectorisé de manière optimale).

%% set DGP parameters:
p = [0.35, 0.8, 3.5, 2.2, 0.3, 2.9, 4.3, 2.1, 0.4, 0.2];
p = p./sum(p); % true probabilities

ncat = numel(p);
cats = 1:ncat;

% draw a sample:
rng(1703) % set seed
nsamp = 10^3; 
samp  = datasample(1:10, nsamp, 'Weights', p, 'Replace', true);

Vérifiez que cela a du sens.

psamp = mean(bsxfun(@eq, samp', cats)); % sample probabilities
bar([p(:), psamp(:)])

entrez la description de l'image ici

Exécutez la simulation d'amorçage.

%% bootstrap simulation:
rng(240947)

nboots = 2*10^3;
cover  = 0.9;
conf   = 0.95;    

tic
Pmat = nan(nboots, ncat, ncat);
for b = 1:nboots

    boot  = datasample(samp, nsamp, 'Replace', true); % draw bootstrap sample    
    pboot = mean(bsxfun(@eq, boot', cats));      

    for l = 1:ncat
        for u = l:ncat
            Pmat(b, l, u) = sum(pboot(l:u));   
        end
    end

end
toc % Elapsed time is 0.442703 seconds.

Le filtre de chaque bootstrap reproduit les intervalles, , qui contiennent au moins masse de probabilité et calcule une estimation de confiance (fréquentiste) de ces intervalles.[l,u]90%

conf_mat = squeeze(mean(Pmat >= cover, 1))

     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    0.3360    0.9770    1.0000
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0

Sélectionnez ceux qui satisfont le desideratum de confiance.

[L, U] = find(conf_mat >= conf);
[L, U]

 1     8
 2     8
 1     9
 2     9
 3     9
 1    10
 2    10
 3    10

Vous convaincre que la méthode d'amorçage ci-dessus est valide

Les échantillons de bootstrap sont destinés à remplacer quelque chose que nous aimerions avoir, mais pas, c'est-à-dire: de nouveaux tirages indépendants de la véritable population sous-jacente (bref: nouvelles données).

Dans l'exemple que j'ai donné, nous connaissons le processus de génération de données (DGP), donc nous pourrions "tricher" et remplacer les lignes de code relatives aux rééchantillons de bootstrap par de nouveaux tirages indépendants du DGP réel.

newsamp = datasample(cats, nsamp, 'Weights', p, 'Replace', true);
pnew    = mean(bsxfun(@eq, newsamp', cats));

Ensuite, nous pouvons valider l'approche bootstrap en la comparant à l'idéal. Voici les résultats.

La matrice de confiance à partir de nouvelles données indépendantes tire:

     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    0.4075    0.9925    1.0000
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0

Les limites inférieures et supérieures de confiance à correspondantes :95%

 1     8
 2     8
 1     9
 2     9
 3     9
 1    10
 2    10
 3    10

On constate que les matrices de confiance s'accordent étroitement et que les bornes sont identiques ... Validant ainsi l'approche bootstrap.

Jim
la source
2
Les intervalles de tolérance et les intervalles de confiance sont des choses différentes. En fait, ce que vous avez décrit n'est pas un intervalle de confiance mais un intervalle de prédiction, qui est encore un autre type d'intervalle distinct.
Kodiologist
1
Votre modification semble être une implémentation de ce que je voulais dire quand j'ai écrit "une méthode Monte Carlo telle que la génération de 1000 répliques bootstrap de l'échantillon et la recherche d'un intervalle qui capture au moins 90% de l'échantillon dans au moins 950 des répliques". Bien qu'intuitif, je ne suis pas sûr que cela fonctionne ou ait du sens, c'est pourquoi j'ai créé cette question.
Kodiologist
@Kodiologist La réponse contient maintenant une section validant l'approche bootstrap. Bien sûr, cela pourrait être poussé plus loin, par exemple imbriqué dans des boucles sur la taille de l'échantillon et les probabilités de classe.
Jim
Montrer que la méthode bootstrap est correcte pour ce problème dans son intégralité signifie montrer qu'elle a la bonne confiance et la bonne couverture indépendamment de la distribution préalable des paramètres (c'est, après tout, une méthode fréquentiste). Pour cela, je pense que la simulation ne suffirait pas; vous auriez besoin d'une preuve mathématique. Mais vous avez été remarquablement persistant et montré que le bootstrap fonctionne au moins parfois, vous méritez donc un A pour l'effort.
Kodiologist