Nichtparametrische Toleranzintervalle für diskrete Variablen

8

Angenommen, Sie haben eine Gruppe von Leuten, die bewerten, wie sehr sie einen Film auf einer diskreten Skala von 1 bis 10 mochten, und Sie möchten ein Intervall [ l , u ], das mit (mindestens) 95% Vertrauen (mindestens) 90 beträgt % aller Personen, die den Film sehen, bewerten ihn nicht niedriger als l und nicht höher als u . [ l , u ] ist dann ein (zweiseitiges) Toleranzintervall mit 95% Konfidenz und 90% Abdeckung. (Um klar zu sein, bedeutet 95% iges Vertrauen, dass bei mehrmaliger Wiederholung dieses Verfahrens 95% der erzeugten Intervalle eine Bevölkerungsabdeckung von mindestens 90% erhalten würden.) Natürlich möchten wir im Allgemeinen, dass [ l , u ] so eng wie möglich ist möglich, während unsere Anforderungen noch erfüllt werden.

Ich habe verschiedene nichtparametrische Methoden zum Erstellen von Toleranzintervallen für kontinuierliche Zufallsvariablen gesehen. Ich habe auch Methoden zum Erstellen von Toleranzintervallen für Binomial- und Poisson-Variablen gesehen. (Das R-Paket toleranceimplementiert mehrere dieser Methoden; Young, 2010.) Aber was ist mit diskreten Variablen, wenn die Verteilung unbekannt ist? Dies ist im Allgemeinen bei Bewertungsskalen wie der in meinem Beispiel der Fall, und die Annahme einer Binomialverteilung scheint nicht sicher zu sein, da echte Bewertungsskalendaten häufig Verrücktheiten wie Multimodalität aufweisen.

Wäre es sinnvoll, auf die nichtparametrischen Methoden für kontinuierliche Variablen zurückzugreifen? Was ist alternativ mit einer Monte-Carlo-Methode wie dem Generieren von 1.000 Bootstrap-Replikaten der Stichprobe und dem Finden eines Intervalls, das mindestens 90% der Stichprobe in mindestens 950 der Replikate erfasst?

Young, DS (2010). Toleranz: Ein R-Paket zur Schätzung von Toleranzintervallen. Journal of Statistical Software, 36 (5), 1–39. Abgerufen von http://www.jstatsoft.org/v36/i05

Kodiologe
quelle
meinst du binomial oder multinomial? Multinomial würde multimodales Verhalten ermöglichen?
Seanv507
Ich meine Binomial. Im Fall einer Bewertungsskala würden Sie beispielsweise die Anzahl der Bernoulli-Versuche auf die Anzahl der Skalenpunkte einstellen. Intervalle über die Kategorien einer multinomialen Verteilung wären meiner Meinung nach nicht sehr sinnvoll, da die Kategorien ungeordnet sind.
Kodiologe
@Kodiologist Ihre Ergebnisvariable ist eine "diskrete Skala von 1 bis 10", aber das bedeutet, dass es sich um eine geordnete multinomiale Antwort handelt. (Oder bekomme ich nichts?)
Jim
@ Jim "Ordered Multinomial" ist ein bisschen wie ein Oxymoron. In einer multinomialen Verteilung ist die Reihenfolge der Kategorien beliebig.
Kodiologe

Antworten:

1

Die interessierende Variable ist multinomial verteilt mit Klassenwahrscheinlichkeiten (Zellenwahrscheinlichkeiten): . Ferner sind die Klassen mit einer natürlichen Ordnung ausgestattet.p1,p2,...,p10

Erster Versuch: kleinstes "Vorhersageintervall" mit90%

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

Ein nichtparametrisches Maß für die Unsicherheit (z. B. Varianz, Vertrauen) in den -quantilen Schätzungen könnte tatsächlich durch Standard-Bootstrap-Methoden erhalten werden .l,u

Zweiter Ansatz: direkte "Bootstrap-Suche"

Im Folgenden stelle ich ausführbaren Matlab-Code zur Verfügung, der sich der Frage direkt aus einer Bootstrap-Perspektive nähert (der Code ist nicht optimal vektorisiert).

%% 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);

Überprüfen Sie, ob dies sinnvoll ist.

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

Geben Sie hier die Bildbeschreibung ein

Führen Sie die Bootstrap-Simulation aus.

%% 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.

Filter aus jedem Bootstrap replizieren die Intervalle , die eine Wahrscheinlichkeitsmasse von mindestens enthalten, und berechnen eine (häufig auftretende) Konfidenzschätzung dieser Intervalle.90 %[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

Wählen Sie diejenigen aus, die das Vertrauens-Desiderat erfüllen.

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

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

Überzeugen Sie sich davon, dass die oben genannte Bootstrap-Methode gültig ist

Bootstrap-Beispiele sollen als Ersatz für etwas dienen, das wir gerne hätten, aber nicht, dh: neue, unabhängige Draws aus der tatsächlichen zugrunde liegenden Population (kurz: neue Daten).

In dem Beispiel, das ich gegeben habe, kennen wir den Datengenerierungsprozess (DGP), daher könnten wir die Codezeilen für Bootstrap-Neuabtastungen durch neue, unabhängige Draws aus dem tatsächlichen DGP "betrügen" und ersetzen.

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

Dann können wir den Bootstrap-Ansatz validieren, indem wir ihn mit dem Ideal vergleichen. Unten sind die Ergebnisse.

Die Konfidenzmatrix aus neuen, unabhängigen Daten zeichnet:

     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

Die entsprechenden -Konfidenz-Unter- und Obergrenzen:95%

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

Wir stellen fest, dass die Konfidenzmatrizen eng übereinstimmen und dass die Grenzen identisch sind ... Dies bestätigt den Bootstrap-Ansatz.

Jim
quelle
2
Toleranzintervalle und Konfidenzintervalle sind verschiedene Dinge. Was Sie beschrieben haben, ist eigentlich kein Konfidenzintervall, sondern ein Vorhersageintervall, das eine weitere Art von Intervall darstellt.
Kodiologe
1
Ihre Bearbeitung scheint eine Implementierung dessen zu sein, was ich meinte, als ich "eine Monte-Carlo-Methode schrieb, wie das Generieren von 1.000 Bootstrap-Replikaten des Samples und das Finden eines Intervalls, das mindestens 90% des Samples in mindestens 950 der Replikate erfasst". Obwohl intuitiv, bin ich mir nicht sicher, ob dies tatsächlich funktioniert oder Sinn macht, weshalb ich diese Frage erstellt habe.
Kodiologe
@Kodiologist Die Antwort enthält jetzt einen Abschnitt, in dem der Bootstrap-Ansatz überprüft wird. Dies könnte natürlich weiter vorangetrieben werden, z. B. verschachtelt in Schleifen über Stichprobengröße und Klassenwahrscheinlichkeiten.
Jim
Um zu zeigen, dass die Bootstrap-Methode in ihrer Gesamtheit für dieses Problem korrekt ist, muss gezeigt werden, dass sie unabhängig von der vorherigen Verteilung der Parameter über das richtige Vertrauen und die richtige Abdeckung verfügt (es handelt sich schließlich um eine häufig verwendete Methode). Dafür würde die Simulation meiner Meinung nach nicht ausreichen. Sie würden mathematische Beweise benötigen. Aber Sie waren bemerkenswert hartnäckig und haben gezeigt, dass der Bootstrap zumindest manchmal funktioniert, sodass Sie ein A für Mühe verdienen.
Kodiologe