Проверка нормального среднего с известным стандартным отклонением, односторонним

Просто, чтобы ввести некоторые концепции, давайте рассмотрим нереально простой пример, где мы хотим проверить среднее и мы знаем стандартное отклонение. Наши данные непрерывны и могут быть смоделированы с помощью нормального распределения. Нам нужно определить размер выборки N, чтобы мы могли различать среднее значение 100 и среднее 110. Мы знаем, что стандартное отклонение составляет 20.

Когда мы проводим статистический тест, мы обычно проверяем нулевую гипотезу против альтернативной гипотезы. Мы находим тестовую статистику T и смотрим на ее распределение под нулевой гипотезой. Если мы наблюдаем необычное значение, скажем, то, которое имеет менее 5% шансов произойти, если нулевая гипотеза верна, то мы отвергаем нулевую гипотезу в пользу альтернативы. (Вероятность 5% известна как уровень значимости теста.) Если значение не является необычным, мы не отвергаем нулевую гипотезу.

В этом случае тестовая статистика T является средним значением выборки. При нулевой гипотезе имеет среднее значение 100 и стандартное отклонение 20/sqrt (N). Сначала рассмотрим фиксированный размер выборки, скажем N = 16. Мы отвергнем нулевую гипотезу, если T находится в затененной области, которая является верхним хвостом его распределения. Это делает этот односторонний тест, поскольку мы не отвергнем, если T находится в нижнем хвосте. Срез для этой затененной области составляет 1,6 стандартных отклонений выше среднего.

rng(0,'twister');
mu0 = 100;
sig = 20;
N = 16;
alpha = 0.05;
conf = 1-alpha;
cutoff = norminv(conf, mu0, sig/sqrt(N));
x = [linspace(90,cutoff), linspace(cutoff,127)];
y = normpdf(x,mu0,sig/sqrt(N));
h1 = plot(x,y);
xhi = [cutoff, x(x>=cutoff)];
yhi = [0, y(x>=cutoff)];
patch(xhi,yhi,'b');
title('Distribution of sample mean, N=16');
xlabel('Sample mean');
ylabel('Density');
text(96,.01,sprintf('Reject if mean>%.4g\nProb = 0.05',cutoff),'Color','b');

Это то, как T ведет себя при нулевой гипотезе, но как насчет альтернативы? Наше альтернативное распределение имеет среднее значение 110, как представлено красной кривой.

mu1 = 110;
y2 = normpdf(x,mu1,sig/sqrt(N));
h2 = line(x,y2,'Color','r');
yhi = [0, y2(x>=cutoff)];
patch(xhi,yhi,'r','FaceAlpha',0.25);
P = 1 - normcdf(cutoff,mu1,sig/sqrt(N));
text(115,.06,sprintf('Reject if T>%.4g\nProb = %.2g',cutoff,P),'Color',[1 0 0]);
legend([h1 h2],'Null hypothesis','Alternative hypothesis');

Существует большая вероятность отклонить нулевую гипотезу, если альтернатива верна. Это то, чего мы хотим. Легче визуализировать, если посмотреть на совокупную функцию распределения (cdf) вместо плотности (pdf). Мы можем считать вероятности непосредственно из этого графика, вместо того, чтобы вычислять области.

ynull = normcdf(x,mu0,sig/sqrt(N));
yalt = normcdf(x,mu1,sig/sqrt(N));
h12 = plot(x,ynull,'b-',x,yalt,'r-');
zval = norminv(conf);
cutoff = mu0 + zval * sig / sqrt(N);
line([90,cutoff,cutoff],[conf, conf, 0],'LineStyle',':');
msg = sprintf(' Cutoff = 100 + %.2g \\times 20 / \\surd{n}',zval);
text(cutoff,.15,msg,'Color','b');
text(min(x),conf,sprintf('   %g%% test',100*alpha),'Color','b',...
                         'verticalalignment','top')
palt = normcdf(cutoff,mu1,sig/sqrt(N));
line([90,cutoff],[palt,palt],'Color','r','LineStyle',':');
text(91,palt+.02,sprintf(' Power is 1-%.2g = %.2g',palt,1-palt),'Color',[1 0 0]);
legend(h12,'Null hypothesis','Alternative hypothesis');

Этот график показывает вероятность получения значительной статистики (отклонение нулевой гипотезы) для двух различных значений mu при N = 16.

Степень определяется как вероятность отклонения нулевой гипотезы, когда альтернатива верна. Это зависит от значения альтернативы и от размера выборки. Мы нарисуем степень (которая является одним минус cdf) как функцию от N, фиксируя альтернативу в 110. Мы выберем N, чтобы достичь степени 80%. График показывает, что нам нужно около N = 25.

DesiredPower = 0.80;
Nvec = 1:30;
cutoff = mu0 + norminv(conf)*sig./sqrt(Nvec);
power = 1 - normcdf(cutoff, mu1, sig./sqrt(Nvec));
plot(Nvec,power,'bo-',[0 30],[DesiredPower DesiredPower],'k:');
xlabel('N = sample size')
ylabel('Power')
title('Power function for the alternative: \mu = 110')

В этом слишком простом примере существует формула, чтобы вычислить необходимое значение непосредственно, чтобы получить степень 80%:

mudiff = (mu1 - mu0) / sig;
N80 = ceil(((norminv(1-DesiredPower)-norminv(conf)) / mudiff)^2)

N80 =

    25

Чтобы убедиться, что это работает, давайте сделаем симуляцию Монте-Карло и сгенерируем 400 выборки размера 25 как при нулевой гипотезе со средним 100, так и при альтернативной гипотезе со средним 110. Если мы проверяем каждую выборку, чтобы увидеть, является ли ее среднее значение 100, мы должны ожидать, что около 5% первой группы и 80% второй группы будут значительными.

nsamples = 400;
samplenum = 1:nsamples;
N = 25;
h0 = zeros(1,nsamples);
h1 = h0;
for j = 1:nsamples
    Z0 = normrnd(mu0,sig,N,1);
    h0(j) = ztest(Z0,mu0,sig,alpha,'right');
    Z1 = normrnd(mu1,sig,N,1);
    h1(j) = ztest(Z1,mu0,sig,alpha,'right');
end
p0 = cumsum(h0) ./ samplenum;
p1 = cumsum(h1) ./ samplenum;
plot(samplenum,p0,'b-',samplenum,p1,'r-')
xlabel('Sample number')
ylabel('Proportion significant')
title('Verification of power computation')
legend('Null hypothesis','Alternative hypothesis','Location','East')

Проверка нормального среднего с неизвестным стандартным отклонением, двусторонним

Теперь предположим, что мы не знаем стандартного отклонения, и мы хотим выполнить двусторонний тест, то есть тот, который отвергает нулевую гипотезу, является ли среднее значение выборки слишком высоким или слишком низким.

Тестовая статистическая величина является статистической величиной t, которая является различием между средней выборкой и средним тестируемым значением, разделенной на стандартную ошибку среднего значения. Согласно нулевой гипотезе, это имеет t- распределения Стьюдента с N-1 степенями свободы. Согласно альтернативной гипотезе, распределение является нецентральным t-распределением с параметром нецентральности, равным нормализованному различию между истинным средним и проверяемым средним.

Для этого двустороннего теста мы должны выделить 5% вероятности ошибки при нулевой гипотезе в равной степени обоим хвостам и отклонить, если тестовая статистика слишком экстремальна в любом направлении. Мы также должны рассмотреть оба хвоста при любой альтернативе.

N = 16;
df = N-1;
alpha = 0.05;
conf = 1-alpha;
cutoff1 = tinv(alpha/2,df);
cutoff2 = tinv(1-alpha/2,df);
x = [linspace(-5,cutoff1), linspace(cutoff1,cutoff2),linspace(cutoff2,5)];
y = tpdf(x,df);
h1 = plot(x,y);
xlo = [x(x<=cutoff1),cutoff1];
ylo = [y(x<=cutoff1),0];
xhi = [cutoff2,x(x>=cutoff2)];
yhi = [0, y(x>=cutoff2)];
patch(xlo,ylo,'b');
patch(xhi,yhi,'b');
title('Distribution of t statistic, N=16');
xlabel('t');
ylabel('Density');
text(2.5,.05,sprintf('Reject if t>%.4g\nProb = 0.025',cutoff2),'Color','b');
text(-4.5,.05,sprintf('Reject if t<%.4g\nProb = 0.025',cutoff1),'Color','b');

Вместо исследования степени как раз при нулевой гипотезе и единственном альтернативном значении mu, мы можем рассматривать ее как функцию mu. Степень увеличивается, когда mu отходит от значения нулевой гипотезы в любом направлении. Мы можем использовать sampsizepwr функция для вычисления степени. Для вычисления степени нам нужно задать значение для стандартного отклонения, которое, как мы подозреваем, будет примерно 20. Вот изображение функции степени для размера выборки N = 16.

N = 16;
x = linspace(90,127);
power = sampsizepwr('t',[100 20],x,[],N);
plot(x,power);
xlabel('True mean')
ylabel('Power')
title('Power function for N=16')

Нам нужна степень 80%, когда среднее значение 110. Согласно этому графику, наша степень менее 50% с размером выборки N = 16. Какой размер образца даст степень, которую мы хотим?

N = sampsizepwr('t',[100 20],110,0.8)

N =

    34

Нам нужен размер выборки около 34. По сравнению с предыдущим примером нам нужно взять девять дополнительных наблюдений, чтобы позволить двухсторонний тест и восполнить незнание истинного стандартного отклонения.

Мы можем сделать график функции степени для различных значений N.

Nvec = 2:40;
power = sampsizepwr('t',[100 20],110,[],Nvec);
plot(Nvec,power,'bo-',[0 40],[DesiredPower DesiredPower],'k:');
xlabel('N = sample size')
ylabel('Power')
title('Power function for the alternative: \mu = 110')

И мы можем сделать симуляцию, подобную предыдущей, чтобы убедиться, что мы получаем степень.

nsamples = 400;
samplenum = 1:nsamples;
N = 34;
h0 = zeros(1,nsamples);
h1 = h0;
for j = 1:nsamples
    Z0 = normrnd(mu0,sig,N,1);
    h0(j) = ttest(Z0,mu0,alpha);
    Z1 = normrnd(mu1,sig,N,1);
    h1(j) = ttest(Z1,mu0,alpha);
end
p0 = cumsum(h0) ./ samplenum;
p1 = cumsum(h1) ./ samplenum;
plot(samplenum,p0,'b-',samplenum,p1,'r-')
xlabel('Sample number')
ylabel('Proportion significant')
title('Verification of power computation')
legend('Null hypothesis','Alternative hypothesis','Location','East')

Предположим, мы можем позволить себе взять выборку размером 50. Предположительно, наша степень для обнаружения альтернативного значения mu = 110 будет больше 80%. Если мы поддерживаем степень на уровне 80%, какую альтернативу мы могли бы обнаружить?

mu1 = sampsizepwr('t',[100 20],[],.8,50)

mu1 =

  108.0837

Проверка пропорции

Теперь давайте обратимся к проблеме определения размера выборки, необходимой для различения двух пропорций. Предположим, что мы отбираем население, в которой около 30% отдают предпочтение некоторому кандидату, и мы хотим отобрать достаточное количество людей, чтобы мы могли отличить это значение от 33%.

Идея здесь такая же, как и раньше. Здесь мы можем использовать отсчет выборки как нашу тестовую статистику. Это количество имеет биномиальное распределение. Для любого размера выборки N мы можем вычислить срез для отклонения нулевой гипотезы P = 0,30. Для N = 100, например, мы бы отвергли нулевую гипотезу, если количество выборок больше, чем значение отсечения, вычисленное следующим образом:

N = 100;
alpha = 0.05;
p0 = 0.30;
p1 = 0.33;
cutoff = binoinv(1-alpha, N, p0)

cutoff =

    38

Осложнение с биномиальным распределением возникает, потому что оно дискретно. Вероятность превышения значения отсечения не совсем 5%:

1 - binocdf(cutoff, N, p0)

ans =

    0.0340

Еще раз, давайте вычислим степень относительно альтернативы P = 0,33 для области значений размеров выборки. Мы будем использовать односторонний (правохвостый) тест, потому что нас интересуют только альтернативные значения, превышающие 30%.

Nvec = 50:50:2000;
power = sampsizepwr('p',p0,p1,[],Nvec,'tail','right');
plot(Nvec,power,'bo-',[0 2000],[DesiredPower DesiredPower],'k:');
xlabel('N = sample size')
ylabel('Power')
title('Power function for the alternative: p = 0.33')

Мы можем использовать sampsizepwr функция, чтобы запросить размер выборки, необходимый для степени 80%.

approxN = sampsizepwr('p',p0,p1,0.80,[],'tail','right')

Warning: Values N>200 are approximate.  Plotting the power as a function
of N may reveal lower N values that have the required power. 

approxN =

        1500

Предупреждающее сообщение сообщает нам, что ответ просто приблизительный. Если мы рассмотрим функцию степени для различных размеров выборки, мы видим, что функция обычно увеличивается, но нерегулярна, потому что биномиальное распределение дискретно. Рассмотрим вероятность отклонения нулевой гипотезы для как p = 0,30, так и p = 0,33 в области значений размеров выборок от 1470 до 1480.

subplot(3,1,1);
Nvec = 1470:1480;
power = sampsizepwr('p',p0,p1,[],Nvec,'tail','right');
plot(Nvec,power,'ro-',[min(Nvec),max(Nvec)],[DesiredPower DesiredPower],'k:');
ylabel(sprintf('Prob[T>cutoff]\nif p=0.33'))
h_gca = gca;
h_gca.XTickLabel = '';
ylim([.78, .82]);

subplot(3,1,2);
alf = sampsizepwr('p',p0,p0,[],Nvec,'tail','right');
plot(Nvec,alf,'bo-',[min(Nvec),max(Nvec)],[alpha alpha],'k:');
ylabel(sprintf('Prob[T>cutoff]\nif p=0.30'))
h_gca = gca;
h_gca.XTickLabel = '';
ylim([.04, .06]);

subplot(3,1,3);
cutoff = binoinv(1-alpha, Nvec, p0);
plot(Nvec,cutoff,'go-');
xlabel('N = sample size')
ylabel('Cutoff')

Этот график показывает, что кривая функции степени (верхний график) не только неправильна, но и уменьшается при некоторых размерах выборки. Это размеры выборки, для которых необходимо увеличить значение среза (нижний график), порядок сохранить уровень значимости (средний график) не более 5%. Мы можем найти меньший размер выборки в этой области значений, который дает желаемую степень 80%:

min(Nvec(power>=0.80))

ans =

        1478

Проверка корреляции

В примерах, которые мы рассмотрели до сих пор, мы смогли вычислить отсечение для тестовой статистики, чтобы достичь определенного уровня значимости и вычислить вероятность превышения этого отсечения в соответствии с альтернативной гипотезой. Для нашего последнего примера давайте рассмотрим проблему, где это не так просто.

Представьте, что мы можем взять выборки из двух переменных X и Y, и мы хотим знать, какой размер выборки нам нужно было бы проверить, являются ли они некоррелированными по сравнению с альтернативой, что их корреляция достигает 0,4. Хотя возможно проработать распределение корреляции выборки путем преобразования его в t-распределение, давайте будем использовать метод, который мы можем использовать даже в задачах, где мы не можем проработать распределение тестовой статистики.

Для заданного размера выборки мы можем использовать симуляцию Монте-Карло, чтобы определить приблизительное значение среза для теста корреляции. Давайте сделаем большой запуск симуляции, чтобы мы могли получить это значение точно. Начнем с выборки размером 25.

nsamples = 10000;
N = 25;
alpha = 0.05;
conf = 1-alpha;
r = zeros(1,nsamples);
for j = 1:nsamples
    xy = normrnd(0,1,N,2);
    r(j) = corr(xy(:,1),xy(:,2));
end
cutoff = quantile(r,conf)

cutoff =

    0.3372

Затем мы можем сгенерировать выборки под альтернативной гипотезой и оценить степень теста.

nsamples = 1000;
mu = [0; 0];
sig = [1 0.4; 0.4 1];
r = zeros(1,nsamples);
for j = 1:nsamples
    xy = mvnrnd(mu,sig,N);
    r(j) = corr(xy(:,1),xy(:,2));
end
[power,powerci] = binofit(sum(r>cutoff),nsamples)

power =

    0.6470


powerci =

    0.6165    0.6767

Мы оцениваем степень в 65%, и у нас есть 95% доверия, что истинное значение находится между 62% и 68%. Чтобы получить степень 80%, нам нужен больший размер выборки. Мы можем попробовать увеличить N до 50, оценить значение среза для этого размера выборки и повторить симуляцию степени.

nsamples = 10000;
N = 50;
alpha = 0.05;
conf = 1-alpha;
r = zeros(1,nsamples);
for j = 1:nsamples
    xy = normrnd(0,1,N,2);
    r(j) = corr(xy(:,1),xy(:,2));
end
cutoff = quantile(r,conf)

nsamples = 1000;
mu = [0; 0];
sig = [1 0.4; 0.4 1];
r = zeros(1,nsamples);
for j = 1:nsamples
    xy = mvnrnd(mu,sig,N);
    r(j) = corr(xy(:,1),xy(:,2));
end
[power,powerci] = binofit(sum(r>cutoff),nsamples)

cutoff =

    0.2315


power =

    0.8990


powerci =

    0.8786    0.9170

Этот размер выборки дает степень лучше, чем наша цель 80%. Мы могли бы продолжить экспериментировать таким образом, пытаясь найти размер выборки менее 50, который соответствовал бы нашим требованиям.

Заключение

Функции вероятности в Statistics and Machine Learning Toolbox могут использоваться, чтобы определить размер выборки, необходимый для достижения желаемого уровня степени в тесте гипотезы. В некоторых задачах размер выборки может быть вычислен непосредственно; в других необходимо искать по области значений размеров выборки до тех пор, пока не будет найдено правильное значение. Генераторы случайных чисел могут помочь проверить, что требуемая степень достигнута, а также могут использоваться, чтобы изучить степень определенного теста в альтернативных условиях.