Документация

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

Только, чтобы ввести некоторые концепции, давайте рассмотрим нереалистично простой пример, где мы хотим протестировать среднее значение, и мы знаем стандартное отклонение. Наши данные непрерывны, и могут быть смоделированы с нормальным распределением. Мы должны определить объем выборки 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 в области значений демонстрационных размеров от 1 470 до 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 and 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 могут использоваться, чтобы определить объем выборки, требуемый достигнуть желаемого уровня степени в тесте гипотезы. В некоторых проблемах объем выборки может быть, вычисляют непосредственно; в других необходимо искать в области значений объемов выборки, пока правильное значение не найдено. Генераторы случайных чисел могут помочь проверить, что желаемой степени соответствуют и можно также использовать, чтобы изучить степень определенного теста при альтернативных условиях.