- Ranks. Behavior under the null and under a shift alternative.
- The Wilcoxon rank-sum test against a shift alternative. Connection to permutation tests including Fisher's exact test.
- The Siegel-Tukey test for an alternative that the dispersions differ.
- The Smirnov test against the omnibus alternative.

References: Lehmann, E.L., 1998. *Nonparametrics: Statistical
Methods Based on Ranks*. Upper Saddle River, N.J.: Prentice Hall;
SticiGui
Chapter 19.

The *rank* of an observation is its position in the list of observations
after sorting the observations from smallest to largest.
The *j*th smallest observation is the observation with rank *j*.
If the data are
{ *t*_{1}, … *t*_{N} },
then *t*_{(j)} is the observation with rank *j*,
for *j*=1, 2, … , *N*.
Thus

*t*_{(1)} ≤ *t*_{(2)} ≤ … ≤
*t*_{(N)}.

Ties can be broken arbitrarily.
The smallest datum is *t*_{(1)} and the largest observation is
*t*_{(N)}.

A large class of nonparametric methods is based on working with the ranks of
the data instead of the original values of the data.
That is, *t*_{(1)} is replaced by 1, *t*_{(2)}
is replaced by 2, and so on.

In some problems, only the ranks of the data are observed. This often happens in situations where it is possible to make comparative judgements about outcomes, but not absolute judgements. For example, a doctor might be able to rank a collection of patients according to the severity with which they suffer from some disease they share, but might not be able to rate the severity on an objective quantitative scale—and certainly not a scale for which 5 is "better than" 4 by the same amount that 2 is "better than" 1. Similarly, a consumer might be able to order a collection of items by preference: she prefers A to B, B to C and C to D, for example. Yet she might not be able to quantify the strength of her preferences. Having techniques for dealing with ranks is then useful.

Working with ranks is also gives a useful replacement for the permutation test
we derived in the previous chapter:
The null distribution of *Y* for the permutation test
depended on the observed responses, but the ranks of *N*
distinct observations are always the integers from 1 to *N*.
Thus if we replace each observation by its rank and then perform a
permutation test, the null distribution of the new test statistic depends only on
*N* and *n*.

In the following exposition, the notation and most of the derivations follow
Lehmann, E.L. (1998. *Nonparametrics: Statistical
Methods Based on Ranks*. Upper Saddle River, N.J.: Prentice Hall).

In the randomization model we have been discussing, let
*S*_{1}, … *S*_{n} be
the ranks of the *n* responses of the treatment group among the
*N* responses of all subjects.
That is, *S*_{1} is the rank of the response of the first treated subject,
*S*_{2} is the rank of the response of the second treated subject,
and so on.
Each *S*_{j} takes a value between 1 and *N*, inclusive.
Assume for the moment that no two responses are equal, so that the ranks are not
ambiguous.
Define

*W*_{s} ≡ *S*_{1} + *S*_{2} +
*S*_{3} +
… + *S*_{n}.

The statistic *W*_{s} is the *Wilcoxon rank-sum*.
We can base a test of the strong null hypothesis on *W*_{s}.
To test against the alternative that treatment tends to increase responses,
we would reject for large values of *W*_{s}.
To test against the alternative that treatment tends to decrease responses,
we would reject for small values of *W*_{s}.
The critical value of the test is set using the probability distribution of
*W*_{s} on the assumption that the strong null hypothesis is
true.
For a level-alpha test against the alternative that treatment increases
responses, we would find the smallest *c* such that, if the
strong null is true, P(*W*_{s} ≥ *c*) ≤ α.
We would then reject the strong null if the observed value of
*W*_{s} is *c* or greater.
Recall that

1 + 2 + … + *k* = *k*(*k*+1)/2.

If the treated subjects have the smallest possible ranks, 1 to *n*,
then

*W*_{s}= 1 + 2 + … + *n* = *n*(*n*+1)/2.

If the treated subjects have the largest possible ranks,
*N*−*n*+1 to *N*, then

*W*_{s}= (*N*−*n*+1) + (*N*−*n*+2) +
… + *N*

= (*N*−*n*) + 1 + (*N*−*n*) + 2 + … +
(*N*−*n*) + *n*

= *n*(*N*−*n*) + (1 + 2 + … + *n*)

= *n*(*N*−*n*) + *n*(*n*+1)/2.

All the integers between *n*(*n*+1)/2 and
*n*(*N*−*n*) + *n*(*n*+1)/2 are possible values
of *W*_{s}.
The null distribution of *W*_{s} under the strong null
hypothesis is symmetric about
*n*(*N*+1)/2.
*W*_{s} under the null hypothesis
is *n*×(*N*+1)/2.
(This also follows from the fact that the expected value of the sample sum of
a simple random sample of size *n* from a box of *N* numbers is
equal to *n* times the mean of the *N* numbers; the mean of the
integers 1 to *N* is (*N*+1)/2.)
The variance of *W*_{s} under the strong null hypothesis
is *m**n*(*N*+1)/12.

Define *W*_{r} to be the sum of the control ranks.
Because the treatment ranks and control ranks together comprise all the ranks,

*W*_{r} + *W*_{s} = 1 + 2 + … + *N*
= *N*(*N*+1)/2.

Thus,

*W*_{r} = *N*(*N*+1)/2 − *W*_{s}.

It follows that the (null) expected value of *W*_{r} is

*N*(*N*+1)/2 − *n*(*N*+1)/2 = *m*(*N*+1)/2,

and that the (null) variance of *W*_{r} is also
*m**n*(*N*+1)/12.
When *n* and *m* are both large, the normal approximations to the
null distributions of
*W*_{r} and to *W*_{s} tend to be accurate.
We can check this by simulation using the sampling applet; see
*N* and *n*.
To check the approximation, first draw 10,000 samples to get an approximation to the
null probability distribution of *W*_{s}.
Then highlight various intervals using the scrollbars in the figure and compare the
area under the histogram with the area under the normal curve.

The Wilcoxon rank-sum test is exactly what you get if you replace each observation by its
rank, then do a permutation test using the sum of the (ranks of the) responses in the
treatment group.
Thus, the code in the previous chapter can be used to simulate the null distribution,
by replacing the raw data with their ranks.
However, the normal approximation to the null distribution
of *W*_{s} can be better than the normal approximation to the null
distribution of the sample sum of the responses to treatment, because the ranks
are evenly spread out, whereas the raw responses can be highly skewed or multimodal.

Define *W*_{XY} ≡ *W*_{s} − *n*(*n*+1)/2.
This subtracts from *W*_{s} its minimum possible value.
Let *W*_{YX} ≡ *W*_{r} − *m*(*m*+1)/2,
*W*_{r} minus its minimum possible value.
Under the strong null hypothesis, the probability distribution of *W*_{XY}
is the same as the probability distribution of *W*_{YX},
a consequence of the symmetry of the probability distribution of *W*_{s}.
The statistics *W*_{XY} and *W*_{YX} are called
the *Mann-Whitney* statistics.

Let { *X*_{1}, … , *X*_{m} } denote
the control responses and let
{ *Y*_{1}, … , *Y*_{n} }
denote the treatment responses as before.
Consider

#{ (*i*, *j*) : 1 ≤ *i*≤ *m*,
1 ≤ *j*≤ *n*, and
*X*_{i} < *Y*_{j} }.

This is the number of (control, treatment) pairs such that the control response is less than the treatment response.

Let { *S*_{j}: *j* = 1, … , *n* }
be the ranks of the treatment responses as
before, and let
{ *R*_{j}: *j* = 1, … , *m* }
be the control responses.
Let { *S*_{(j)}: *j* = 1, … , *n* }
and { *R*_{(j)}: *j* = 1, … , *m* }
be the corresponding ordered ranks.
Use the value of the treatment response to partition the set of (control, treatment)
pairs for which the control response is less than the treatment response: The total number of
such pairs for which the control response is less than the treatment response
is the number of pairs where the control response is less than the smallest treatment
response, plus the number of pairs where the control response is less than the second-smallest
treatment response, and so on.
This will help us count the pairs.

How many response values are less than *S*_{(1)}, the rank of the smallest
treatment response?
By definition, there are *S*_{(1)}−1 of them—all of which are control
responses.
The number of response values that are less than *S*_{(2)} is
*S*_{(2)}−1, one of which is *S*_{(1)}, so the
total number of control responses that are less than *S*_{(2)} is
*S*_{(2)}−2.
The total number of control responses that are less than *S*_{(j)}
is *S*_{(j)} − *j*, so the total number of (control, treatment)
pairs with the control response less than the treatment response is

*S*_{(1)}−1 + *S*_{(2)}−2 + … +
*S*_{(n)} − *n* = *W*_{s} −
(1 + 2 + … + *n*) =
*W*_{s} − *n*(*n*+1)/2 = *W*_{XY}.

Thus

*W*_{XY} =
#{ (*i*, *j*) : 1 ≤ *i*≤ *m*,
1 ≤ *j*≤ *n*, and
*X*_{i} < *Y*_{j} }.

The Mann-Whitney statistic *W*_{XY} (and the Wilcoxon rank sum
*W*_{s}, up to an additive constant)
measures the number of (control, treatment) pairs for which the treatment response
is at least as large as the control response.
The larger the positive effect of treatment, the larger the Mann-Whitney and Wilcoxon
rank sum statistics tend to be.

If the control and treatment responses are random samples from continuous distributions,
the chance of ties among the data is zero.
However, ties can and do occur in practice, if only because of limits on the
precision with which data can be recorded.
When there are tied observations, the ranks are not uniquely defined.
We can patch the rank-sum approach by assigning each set of tied observations
the *mid-rank* of the set.
For example, if the sorted data are {1, 2, 3, 4, 4} as above, the mid-rank of the
last two observations (tied at 4) would be 4.5.
If the sorted data were {1, 2, 3, 3, 3, 4, 4, 5}, the corresponding mid-ranks would
be 1, 2, 4, 4, 4, 6.5, 6.5, and 8.
The Wilcoxon rank-sum statistic generalized to use mid-ranks is denoted
*W*_{s}^{*}.
Once there are ties, the null distribution of the rank sum depends on the observations:
which mid-ranks are represented.
When there are ties, the normal approximation to the null probability distribution
of the rank-sum tends to be worse, because there are fewer possible values of
the rank sum (the null probability distribution is "chunkier").

*N*=5, *n*=2.
Draw 10,000 samples to see the simulated null distribution of the rank sum.
Note that the distribution of the rank-sum is skewed in this case.
Compare the simulated null distribution of the rank-sums with and without
midranks by replacing the contents of the box by {1, 2, 3, 4, 5}, which
would be the ranks if there were no ties in the data.
Note that the distribution is then symmetric.
Replace the contents of the box with the mid-ranks
1, 2, 4, 4, 4, 6.5, 6.5, and 8 and change the sample size (*n*) to 4.
Compare the simulated sampling distribution of the rank sum with that
for ranks 1, 2, 3, 4, 5, 6, 7, 8,
and check the accuracy of the normal approximation in both cases.

Calculating midranks seems to involve a number of steps: sort the observations,
identify ties, and average the ranks of each group of ties.
An alternative point of view makes the computation much simpler.
If an observation is not tied, its rank is equal to the number of
observations that are less than or equal to it.
That is, suppose the data are *z*_{1}, … , *z*_{N}.
If the multiplicity of the value *z*_{k} is one,
then

rank(*z*_{k}) = #{ *j*: *z*_{j} ≤
*z*_{k} }.

On the other hand, if there are *M*(*z*_{k})
observations in all whose values are
equal to *z*_{k}, then the midrank of *z*_{k}
is

midrank(*z*_{k}) =
#{ *j*: *z*_{j} ≤ *z*_{k} } −
(*M*(*z*_{k})−1)/2.

We can calculate *M*(*z*_{k}) by

*M*(*z*_{k}) = #{ *j*: *z*_{j} ≤
*z*_{k} } −
#{ *j*: *z*_{j} < *z*_{k} }.

Thus

midrank(*z*_{k}) = (
#{ *j*: *z*_{j} ≤ *z*_{k} } +
#{ *j*: *z*_{j} < *z*_{k} } +
1 )/2.

Here is a Matlab function to compute the midranks of a list:

function r = midRanks(x) % function r = midRanks(x) % P.B. Stark, statistics.berkeley.edu/~stark 2/17/2003 % vector of midranks of a vector x. % midrank is (#<=) - ((#=) - 1)/2 = ((#<=) + (#<) + 1)/2 for j = 1:length(x), r(j) = sum(x <= x(j)) - (sum(x == x(j)) - 1)/2; end; return;

Here is an R function to do the same thing:

midRanks <- function(x) { # P.B. Stark, statistics.berkeley.edu/~stark 2/5/2008 # vector of midranks of a vector x. mr <- array(0,length(x)); for (j in 1:length(x)) { mr[j] <- sum(x <= x[j]) - (sum(x == x[j]) - 1)/2; } mr }

Here is a Matlab function to simulate the distribution of the sample
sum of a simple random sample of size *n* from a population of size *N*.

function dist = simSampleSum(z, n, iter) % function dist = simSampleSum(z, n, iter) % P.B. Stark, statistics.berkeley.edu/~stark 2/17/2003 % Simulate the sampling distribution of the sum of n of the elements of z dist = zeros(iter,1); N = length(z); if (n > N) disp(['error in simSampleSum: sample size exceeds population size']); return; elseif (n == N) dist = dist + sum(z); % constant random variable elseif (n == 0) return; % zero else for i=1:iter zp = z(randperm(length(z))); % random permutation of data dist(i) = sum(zp(1:n)); % add the first n end; end; return;

Here is a Matlab function to simulate the null distribution of the Wilcoxon rank sum using midranks.

function dist = simWilcox(x, y, iter) % function dist = simWilcox(x, y, iter) % P.B. Stark, statistics.berkeley.edu/~stark 2/17/2003 % simulates the null distribution of the Wilcoxon rank sum using midranks, % for data x (control) and y (treatment) using iter pseudorandom samples z = midRanks([x, y]); n = length(y); dist = simSampleSum(z, n, iter); return;

Here's an example of using the Matlab functions just defined to approximate
the 1-sided (upper tail) *P*-value for the Wilcoxon rank sum statistic
with midranks:

iter = 10000; zr = midRanks([x y]); % the control responses are x; the treatment % responses are y. zr now has their midranks. dist = simWilcox(x, y, iter); pVal = sum(dist >= sum(zr[(length(x)+1):length(zr)]) )/iter;

When the alternative is that the dispersion of the control responses differs
from the dispersion of the treatment responses, the Wilcoxon rank-sum statistic
has poor power.
We could imagine a permutation test based on the inter-quartile range of the
responses in the treatment group or the standard deviation of the responses
of the treatment group; under the strong null hypothesis, we could find
the distribution of that test statistic by calculation or simulation.
The randomization, in which every subset of *n* of the *N*
responses equally is likely to be the treatment responses, makes such calculations
and tests straightforward, at least conceptually.

There is a clever way to relabel the data that allows us to use the
calculations and tables for the Wilcoxon rank-sum statistic to test against
the dispersion alternative—we just redefine what we mean by "rank."
Suppose that the alternative is that the treatment decreases the dispersion of
the responses around some common measure of location (such as the median).
Then responses far from the middle rank should be less likely to occur in the
treatment group than in the control group.
Label the smallest observation "1", the largest "2",
the second-smallest "3", the second-largest "4", and so on,
ignoring the possibility of ties for the moment.
Consider the sum *T* of the *n* labels of the treated subjects.
The null distribution of *T* is the same as the null distribution of
the Wilcoxon rank-sum statistic *W*_{s}, and under the alternative
that treatment decreases dispersion, *T* tends to be large.
Therefore, rejecting the strong null hypothesis when

*T* > *c*

where *c* is the appropriate quantile of the null distribution of the
Wilcoxon rank-sum statistic, is a reasonable test
against this alternative.
This is the *Siegel-Tukey* test.

However, labeling the largest observation "1", the smallest "2",
the second-largest "3", the second-smallest "4", and so on,
leads to a test statistic with the same distribution under the strong null hypothesis,
but that typically gives a different *P*-value for a given set of data.
A different test statistic that is more symmetric is to assign the smallest and largest
observations the label 1.5, the second-smallest and second-largest the label 3.5,
the third-smallest and third-largest 5.5, and so on.
Ties can be treated as they are in the generalization of the Wilcoxon rank-sum
statistic *W*_{s}^{*}, by using "mid-labels"
(just like mid-ranks).
To find significance levels or *P*-values when *N* is large,
it is sometimes necessary to resort to simulation and settle for an approximation.

The Siegel-Tukey test tacitly assumes that the treatment responses and the control responses are scattered about a common typical value—that is, that there is no shift difference between the two groups. If there is a known shift between control and treatment, it could be subtracted before ranking the data. If the shift is unknown, it could be estimated from the data, but then the nominal significance level of the test will not be its real level. For large samples, the difference might be unimportant.

This material is from Lehmann (1998).

We have considered shift and dispersion alternatives; now we move to the
omnibus alternative.
The Smirnov test is based on the difference between the
empirical cumulative distribution function (cdf) of the treatment responses
and the empirical cdf of the control responses.
It has some power against all kinds of violations of the strong null hypothesis.
However, it has less power than *W*_{s} against
the alternatives for which the Wilcoxon rank-sum test is designed.
Let *F*_{Y,n}
denote the empirical cdf of the treatment responses:

*F*_{Y,n}(*y*)
≡ #{ *y*_{j}:
*y*_{j} ≤ *y* }/*n*,

and define *F*_{X,m} analogously as the empirical cdf of the
control responses.
If there are no ties within the treatment group, *F*_{Y,n}
jumps by 1/*n* at each treatment response value.
If there are no ties within the control group,
*F*_{X,m} jumps by 1/*m* at each control response value.
If *k* of the treatment responses are tied and
equal to *y*_{0}, then *F*_{Y,n} jumps by *k*/*n*
at *y*_{0}.
The Smirnov test statistic is

*D*_{m,n} ≡ sup |*F*_{Y,n}(*y*)
− *F*_{X,m}(*y*) |,

where the supremum is over all real values of *y*.
It is easy to see that the supremum is attained at one of the data values.
We can also see that
the supremum depends only on the ranks of the data, because the
order of the jumps matters, but the precise values of *y* at which the jumps occur
do not matter.
Therefore, the test

Reject if *D _{m}*

for an appropriately chosen value of *c*, is a nonparametric test of
the strong null hypothesis.

Let's calculate the null distribution of *D _{m}*

The null probability distribution of *D _{m}*

Thus a Smirnov test (against the omnibus alternative) with
significance level 0.2 would reject the strong null hypothesis when
*D*_{2,3}=1;
smaller significance levels are not attainable when *n* and *m* are so small.
Critical values can be calculated analytically fairly easily when
*n*=*m* (when the treatment and control groups are the same size)
and the significance level is an integer multiple *a* of 1/*n*.
Let *k* = ⌊ *n*/*a* ⌋ be the largest integer such that *n*
− *k*×*a* ≥ 0.
Then, under the strong null hypothesis, provided there are no ties,

P(*D _{m}*

When there are ties, *D _{m}*

Here is a Matlab m-file to calculate the Smirnov
statistic.
Let *x*=(*x*_{1}, … , *x*_{m})
denote the vector of control responses and let
*y*=(*y*_{1}, … , *y*_{n})
denote the vector of treatment responses.

function s = smirnov(x, y) % function s = smirnov(x, y) % P.B. Stark, statistics.berkeley.edu/~stark 2/17/2003 % calculates the Smirnov distance between two vectors of data m = length(x); n = length(y); z = [x y]; s = 0; for j=1:m+n, s = max(s, abs(sum(x <= z(j))/m - sum(y <= z(j))/n)); end; return;

The following Matlab function simulates the null distribution of the Smirnov statistic.
Again, the data are *x* and *y*.

function dist = simSmir(x, y, iter) % function dist = simSmir(x, y, iter) % P.B. Stark, statistics.berkeley.edu/~stark 2/17/2003 % simulates the null distribution of the Smirnov distance for data x and y % using iter pseudorandom samples dist = zeros(iter,1); m = length(x); n = length(y); z = [x y]; % pool the N = m+n observations for i=1:iter zp = z(randperm(length(z))); % random permutation of the data dist(i) = smirnov(zp(1:m),zp(m+1:m+n)); % first m control, last n treatment end; return;

Using these functions, you can find an approximate *P*-value for
the Smirnov test using the following script, which assumes that the data are
already in the vectors *x* and *y* and that *iter* has
been set (to a number like 10,000).

testVal = smirnov(x, y); simDist = simSmir(x, y, iter); pValue = sum(simDist >= testVal)/iter;

Here are versions of the functions in R:

smirnov <- function(x, y) { # P.B. Stark, statistics.berkeley.edu/~stark 9/8/2006 z <- c(x,y); s <- 0; for (j in 1:length(z)) { s <- max(s, abs(sum(x <= z[j])/length(x) - sum(y <= z[j])/length(y))); } s }

simSmir <- function(x, y, iter) { # P.B. Stark, statistics.berkeley.edu/~stark 9/8/2006 dist <- array(0, iter); z <- c(x,y); for (i in 1:iter) { zp = sample(z); dist[i] <- smirnov(zp[1:length(x)],zp[(length(x)+1):length(z)]); } dist }

testVal <- smirnov(x, y); simDist <- simSmir(x, y, iter); pValue <- sum(simDist >= testVal)/iter;

Here is a more cryptic "seven at one blow" version:

simSmirnovTest <- function(x, y, iter) { # P.B. Stark, statistics.berkeley.edu/~stark 9/8/2006 ts <- smirnov(x,y) z <- c(x, y) # pooled responses sum(replicate(iter, { zp <- sample(z); (smirnov(zp[1:length(x)],zp[(length(x)+1):length(z)]) >= ts) } ) )/iter }

We can also estimate the power of the Smirnov statistic against a shift alternative
by simulation.
Suppose that the effect of treatment is to increase each subject's response by *d*,
no matter what the response would have been in control.
(Remember, this is an absurdly restrictive assumption, but the results might be suggestive anyway.)
Then the responses *y*=(*y*_{1}, … *y*_{n})
of the treated subjects would have been

*y*−*d*=(*y*_{1} − *d*, …,
*y*_{n} − *d*)

had they been in the control group, and the responses *x* of the control group would
have been *x*+*d* had they been in the treatment group.
The following Matlab function simulates the distribution of the Smirnov statistic under
the shift alternative. It assumes that *x*, *y* and *d* are already
defined.

function altD = simSmirShift(x, y, d, iter) % function altD = simSmirShift(x, y, d, iter) % P.B. Stark, statistics.berkeley.edu/~stark 2/17/2003 % simulates the distribution of the Smirnov statistic under the % shift alternative. x and y are the data, d is the shift. % Uses iter pseudorandom samples. altD = zeros(iter,1); m = length(x); n = length(y); z = [x y-d]; % pool observations. Under the shift % alternative, treatment effect is d for i=1:iter zp = z(randperm(length(z))); % random permutation of data altD(i) = smirnov(zp(1:m),zp(m+1:m+n)+d); % first m control; last n treatment % treatment is shifted by d end; return;

Note that this reduces to the null distribution when *d*=0, so we really
don't need the function simSmir—we can just call simSmirShift with *d*=0.
To find the power of the test, we need to remember that the shift can change the number
of ties, depending on which subjects are assigned to treatment and which to control.
Thus, for each assignment, we need (in principle) to calculate a new critical value for
the test.
The following Matlab script does that.

function p = simSmirPower(x, y, d, alpha, iter) % function p = simSmirPower(x, y, d, alpha, iter) % % P.B. Stark, statistics.berkeley.edu/~stark 3/3/03 % finds the approximate power of a level alpha Smirnov % test against the alternative hypothesis that treatment % shifts the response by exactly d for all subjects. % % x is the vector of control observations % y is the vector of treatment observations % the lengths of x and y determine the sizes of the control and treatment % groups % alpha is the significance level % iter is the number of replications of each simulation--note: there are % two nested simulations. p = 0; z = [x y-d]; for i=1:iter dist = zeros(1,iter); zp = z(randperm(length(z))); xp = zp(1:length(x)); yp = zp(length(x)+1:length(x)+length(y)) + d; dist = simSmirShift(xp, yp, 0, iter); if (sum(dist >= smirnov(xp, yp))/iter <= alpha) p = p+1; end; end; p = p/iter; return;

The following Matlab script finds a critical value for a conservative level *a*;
test, the *P*-value, and the power of the test against a shift alternative with
shift *d*.
It assumes that *x*, *y*, *a*, *d* and *iter*
are already defined.
It uses the function upperCriticalVal(), whose definition follows.

testVal = smirnov(x, y); % value of the test statistic simDist = simSmirShift(x, y, 0, iter); % simulated null distribution critVal = upperCriticalVal(simDist, a); % find critical value for the test pValue = sum(simDist >= testVal)/iter; % approximate p-value power = simSmirPower(x, y, d, alpha, iter); % estimated power

function cv = upperCriticalVal(list, alpha) % function upperCriticalValue(list, alpha) % % P.B. Stark statistics.berkeley.edu/~stark 2/22/08 % Finds a value cv in the list such that at most a % fraction alpha of elements of the list are greater % than or equal to cv n = length(list); if (alpha <= sum(list == max(list))/n | alpha > 1) disp(['Error in upperCriticalVal: unattainable fraction ' num2str(alpha)]); cv = NaN; elseif (alpha == 1) cv = min(list); else list = sort(list); trialInx = ceil(n*(1-alpha)); % first guess at critical value, % but there could be duplicates. cv = list(trialInx); % first guess at critical value sigAttain = sum(list >= cv)/n; % attained estimated significance while (sigAttain > a & trialInx < n) % for conservative level-a test trialInx = trialInx + 1; cv = list(trialInx); sigAttain = sum(list >= cv)/n; end; end; return;