Wie spezifiziert man eine a priori Korrelation zwischen zufällig gezogenen Stichproben aus zwei multinomialen Verteilungen?

Question

Betrachten Sie das folgende Spiel: in jedem Versuch werden Sie mit x roten und y blauen Punkten dargestellt. Sie müssen entscheiden, ob es mehr rote als blaue Punkte gibt. Für jeden Versuch ist die minimale Anzahl von Punkten in einer gegebenen Farbe 10, das Maximum ist 50. Rote und blaue Punkte folgen einer identischen multinomialen Verteilung (Betrachten wir zur Vereinfachung, dass die Wahrscheinlichkeit des Auftretens jeder ganzen Zahl zwischen 10 und 50 ähnlich ist).

Ich möchte 300 Versuche bauen. Um dies zu tun, zeichne ich 300 Proben von jeder multinomialen Verteilung. Wichtig ist, dass ich (a priori) die Korrelation zwischen den 300 Proben aus der ersten Verteilung und den 300 Proben aus der zweiten Verteilung angeben möchte. Ich würde gerne eine Korrelation von -0,8, -0,5, 0, 0,5 und 0,8 in fünf Paaren von Stichprobensätzen haben.

Vorzugsweise möchte ich auch diese Sätze so, dass in jedem Satz (X, Y) mit einer der angegebenen Korrelationen, die Hälfte der X-Proben größer als Y (x(i) > y(i)), und die andere Hälfte wird kleiner sein als Y (x(i) < y(i)).

Wie kann ich das in Python, R oder MATLAB?

Answer 1

Grundsätzlich fragen Sie, wie create 2 vectors with a specified correlation, so dass es mehr Statistiken als programing Frage ist, aber es kann auf folgende Weise erfolgen:

Schritt 1 - Schaffung von zwei Vektor mit der gewünschten Korrelation

r = 0.75;    % r is the desired correlation 
M = rand(10000,2);  % two vectors from uniform distribution between 0 to 1 
R = [1 r; r 1]; 
L = chol(R);    % this is Cholesky decomposition of R 
M = M*L;     % when multiplied by M it gives the wanted correlation 
M = (M+abs(min(M(:)))); % shift the vector to only positive values 
M = M./max(M(:));  % normalize the vector... 
M = round(40*M)+10;  % ...to values between 10 to 50 
disp([min(M(:)) max(M(:))]) 
first_r = corr(M(:,1), M(:,2));  % and check the resulted correlation 

Die rand-Funktion könnte in jede beliebige Zufallszahlenfunktion wie randi oder randn geändert werden, und wenn eine bestimmte Verteilung gewünscht wird, könnte sie using the it's cdf erhalten werden.

Schritt 2 - Probenahme diese Vektoren für zwei Sätze von Proben, eine mit x> y und eine mit y> x

x = M(:,1); 
y = M(:,2); 
Xy = x>y;    % logical index for all x > y 
Yx = y>x;    % logical index for all y > x 
xy1 = datasample([x(Xy) y(Xy)],150,'Replace',false); % make a 1/2 sample like Xy 
xy2 = datasample([x(Yx) y(Yx)],150,'Replace',false); % make a 1/2 sample like Yx 
x = [xy1(:,1);xy2(:,1)];   % concat the smaples back to x 
y = [xy1(:,2);xy2(:,2)];   % concat the smaples back to y 
checkx = sum(x>y)     % how many times x is bigger than y 
checky = sum(y>x)     % how many times y is bigger than x 
final_r = corr(x,y)    % and check the new correlation 

Schritt 3 - Korrigieren des Korrelations

Wie Sie sehe die final_r ist nicht wie die gewünschte r, so dass Sie die erste r durch ihre Entfernung von der final_r zu verschieben haben. Hier ist ein Beispiel - zunächst der Ausgang, wenn r = 0.75:

10 50 
checkx = 
    150 
checky = 
    150 
final_r = 
     0.67511 

wir sehen, dass die final_r durch 0,074886 nach unten verschoben, so wollen wir die ursprünglichen r bis um diesen Wert zu verschieben, um unsere final_r korrekt zu erhalten. Also, wenn wir es wieder mit r = 0.75+0.074886 laufen, erhalten wir:

10 50 
checkx = 
    150 
checky = 
    150 
final_r = 
     0.76379 

die auf die gewünschte r ziemlich nahe ist. Ich würde eine Schleife über den Prozess für sagen sagen, 1000 Mal, um die nächste r zu dem gewünschten zu finden, oder einfach einen Schwellenwert setzen, die weiter suchen, bis die final_r nahe genug ist, was Sie wollen.