Prueba de independencia de Matlab

Question

Para 1,000,000 de observaciones, observé un evento discreto, X, 3 veces para el grupo de control y 10 veces para el grupo de prueba.

Necesito realizar una prueba de independencia de Chi cuadrado en Matlab. Esta es la forma en que lo haría en r:

función

m <- rbind(c(3, 1000000-3), c(10, 1000000-10)) 
#  [,1] [,2] 
# [1,] 3 999997 
# [2,] 10 999990 
chisq.test(m) 

El R Devuelve-chi cuadrado = 2,7692, df = 1, p-valor = 0,0961.

¿Qué función de Matlab debo usar o crear para hacer esto?

Answer 1

Aquí es mi propia aplicación que utilizo:

function [hNull pValue X2] = ChiSquareTest(o, alpha) 
    %# CHISQUARETEST Pearson's Chi-Square test of independence 
    %# 
    %# @param o   The Contignecy Table of the joint frequencies 
    %#      of the two events (attributes) 
    %# @param alpha  Significance level for the test 
    %# @return hNull  hNull = 1: null hypothesis accepted (independent) 
    %#      hNull = 0: null hypothesis rejected (dependent) 
    %# @return pValue The p-value of the test (the prob of obtaining 
    %#      the observed frequencies under hNull) 
    %# @return X2  The value for the chi square statistic 
    %# 

    %# o:  observed frequency 
    %# e:  expected frequency 
    %# dof: degree of freedom 

    [r c] = size(o); 
    dof = (r-1)*(c-1); 

    %# e = (count(A=ai)*count(B=bi))/N 
    e = sum(o,2)*sum(o,1)/sum(o(:)); 

    %# [ sum_r [ sum_c ((o_ij-e_ij)^2/e_ij) ] ] 
    X2 = sum(sum((o-e).^2 ./ e)); 

    %# p-value needed to reject hNull at the significance level with dof 
    pValue = 1 - chi2cdf(X2, dof); 
    hNull = (pValue > alpha); 

    %# X2 value needed to reject hNull at the significance level with dof 
    %#X2table = chi2inv(1-alpha, dof); 
    %#hNull = (X2table > X2); 

end 

Y un ejemplo para ilustrar:

t = [3 999997 ; 10 999990] 
[hNull pVal X2] = ChiSquareTest(t, 0.05) 

hNull = 
    1 
pVal = 
    0.052203 
X2 = 
     3.7693 

Tenga en cuenta que los resultados son diferentes a la suya, porque chisq.test realiza una corrección por defecto, de acuerdo a ?chisq.test

correcto: a logical indicando si para aplicar la corrección de continuidad al calcular el estadístico de prueba para tablas 2x2: la mitad es restado de todos | O - E | diferencias

---

Alternativamente, si usted tiene las observaciones reales de los dos eventos en cuestión, puede utilizar la función CROSSTAB que calcula la tabla de contingencia y devolver los ji 2 y p-valor medidas:

X = randi([1 2],[1000 1]); 
Y = randi([1 2],[1000 1]); 
[t X2 pVal] = crosstab(X,Y) 

t = 
    229 247 
    257 267 
X2 = 
    0.087581 
pVal = 
     0.76728 

el equivalente en I sería:

chisq.test(X, Y, correct = FALSE) 

Nota: Tanto (MATLAB) se aproxima por encima de requerir las estadísticas Caja de herramientas

Answer 2

Esta función probará la independencia mediante el estadístico de chi-cuadrado de Pearson y la estadística de razón de verosimilitud, junto con el cálculo de los residuos. Sé que esto se puede vectorizar más, pero estoy tratando de mostrar las matemáticas para cada paso.

function independenceTest(data) 
df = (size(data,1)-1)*(size(data,2)-1); % Mean Degrees of Freedom 
sd = sqrt(2*df);      % Standard Deviation 

u   = nan(size(data)); % Estimated expected frequencies 
p   = nan(size(data)); % Values used to calculate chi-square 
lr  = nan(size(data)); % Values used to calculate likelihood-ratio 
residuals = nan(size(data)); % Residuals 

rowTotals = sum(data,1); 
colTotals = sum(data,2); 
overallTotal = sum(rowTotals); 

%% Calculate estimated expected frequencies 
for r=1:1:size(data,1) 
    for c=1:1:size(data,2) 
     u(r,c) = (rowTotals(c) * colTotals(r))/overallTotal; 
    end 
end 

%% Calculate chi-squared statistic 
for r=1:1:size(data,1) 
    for c=1:1:size(data,2) 
     p(r,c) = (data(r,c) - u(r,c))^2/u(r,c); 
    end 
end 
chi = sum(sum(p)); % Chi-square statistic 

%% Calculate likelihood-ratio statistic 
for r=1:1:size(data,1) 
    for c=1:1:size(data,2) 
     lr(r,c) = data(r,c) * log(data(r,c)/u(r,c)); 
    end 
end 
G = 2 * sum(sum(lr)); % Likelihood-Ratio statisitc 

%% Calculate residuals 
for r=1:1:size(data,1) 
    for c=1:1:size(data,2) 
     numerator = data(r,c) - u(r,c); 
     denominator = sqrt(u(r,c) * (1 - colTotals(r)/overallTotal) * (1 - rowTotals(c)/overallTotal)); 
     residuals(r,c) = numerator/denominator; 
    end 
end