# Tests, Power and Significance

In the mathematical statistics course today, we started talking about tests, and decision rules. To illustrate all the concepts introduced today, we considered the case where we have a sample $\boldsymbol{x}=\{x_1,\cdots,x_n\}$ with $U_i\sim\mathcal{U}([0,\theta])$. And we want to test

$H_0:\theta\leq \theta_0$  against $H_1:\theta> \theta_0$

In the course, we’ve seen that we could use a test based on the order statistics $x_{n:n}=\max\{x_i\}$.  The test would be

$\psi(\boldsymbol{x})=\boldsymbol{1}_{(c\ ;\ +\infty)}(x_{n:n})$

i.e. if $\psi(\boldsymbol{x})=1$ we choose $H_1$, and if $\psi(\boldsymbol{x})=0$, we choose $H_0$.

From the definition of the first order risk,

$\alpha=\sup_{\theta\in\Theta_0}\left\lbrace{\mathbb E}_{\theta}\lbrack \psi(\boldsymbol{X})\rbrack\right\rbrace={\mathbb E}_{\theta_0}\lbrack \psi(\boldsymbol{X})\rbrack$

we can easily get that

$c=\theta_0\cdot(1-\alpha)^{\frac{1}{n}}$

Thus, the power is then

$p(\theta)=\left(1-\left(\frac{\theta_0}{\theta}\right)^n(1-\alpha)\right)\boldsymbol{1}_{(c\ ;\ +\infty)}(\theta)$

To visualize it, use the following parameters

n=5
alpha=.1
theta0=1

Then

C1=theta0*(1-alpha)^(1/n)
theta=seq(0,2,by=.01)
P1=(1-(theta0/theta)^n*(1-alpha))*(theta>C1)
plot(theta,P1,type="l",lwd=2,col="blue",xlab="",ylab="Power")

Note that, so far, we did never consider the maximum of our sample. Assume that the maximum is $x_{n:n}$, then we can compute the $p$-value,

$p=\mathbb{P}(X_{n:n}>x_{n:n})=1-\left(\frac{x_{n:n}}{\theta_0}\right)^n$

Here it is

PV=(1-theta^n)*(theta<=1)
plot(theta,PV,type="l",lwd=2,col="blue",xlab="",ylab="p-value")

Now, why not consider another test, based on the minimum (since we have the distribution of the minimum of a sample from a uniform distribution). The test is the same as before

$\psi(\boldsymbol{x})=\boldsymbol{1}_{(c\ ;\ +\infty)}(x_{1:n})$

but here, the threshold is

$c=\theta_0\cdot (1-\alpha^{\frac{1}{n}})$

The power of the test is here

$\quad p(\theta)=\left(1-\frac{\theta_0}{\theta}(1-\alpha^{\frac{1}{n}})\right)^n \boldsymbol{1}_{( c\ ;\ +\infty)}(\theta)$

This test has the same significance level (by construction), but the power of the test is clearly lower than the one we got using the maximum of our sample, when $\theta\in\Theta_1$

C2=theta0*(1-alpha^(1/n))
P2=(1-(theta0/theta)*(1-alpha^(1/n)))^n*(theta>C2)
lines(theta,P2,type="l",lwd=2,col="red")

Why not consider a test based on $\overline{x}$? The problem is that we need the distribution (more specifically the survival distribution) of $\overline{X}$. We can compute it, numerically. But that might be painful. An alternative is to consider some approximation, based on the central limit theorem, i.e.

$2\overline{X}\sim\mathcal{N}\left(\theta,2^2 \frac{\theta^2}{12n}\right)$

Our test is based on $\psi(\boldsymbol{x})=\boldsymbol{1}_{(c\ ;\ +\infty)}(2\overline{x})$, and to get the same significance as before, use

$c=\Phi_{\mu,\sigma^2}^{-1}(1-\alpha)=\mu+\sigma \Phi^{-1}(1-\alpha)$

The power of the test is then

$p(\theta)=1-\Phi_{\theta,\sigma^2}(c)\cdot \boldsymbol{1}_{(c,+\infty)}(\theta)$

Here it is

mu=2*(theta0/2)
s2=2^2*(theta0^2/12)/n
C3=qnorm(1-alpha,mu,sqrt(s2))
(P=1-pnorm(C3,theta,sqrt(s2)))*(theta>C3)
lines(theta,P)

Observe here that the test based on the maximum is not more powerful than the one based on the average (I just wonder if it could be due to the Gaussian approximation…).

# Playing cards in Vegas?

In a previous post, a few weeks ago, I mentioned that I will be in Las Vegas by the end of July. And I took the opportunity to write a post on roulette(s). Since some colleagues told me I should take some time to play poker there, I guess I have to understand how to play poker… so I went back to basics on cards, and shuffling techniques.

Now, I have to confess that I have been surprised, while I was looking for mathematical models for shuffling, to find so many deterministic techniques (and results related to algebra, and cycles).

On http://mathworld.wolfram.com/ for instance, one can find nice articles on so-called in-shuffle or out-shuffle techniques. There is also a great article, Golomb (1961), but mainly on algebraic properties of permutations by cutting and shuffling, as well as Diaconis, Kantor and Graham’s The Mathematics of Perfect Shuffle And if you look at Monge’s shuffle, you can find a deterministic recursive relationship. As a statistician (or applied probabilist), I should confess that I did not find answer to the question I wanted to ask : how long should we shuffle before getting cards randomly sorted in ours hands ?

• Randomness (from a statistician perspective)

First, I need to define (as properly as possible) a notion of “cards randomly sorted“. Consider a game with 32 cards. Why 32 ? Mathematicians will tell you that 32 is a great number, since it is a power of 2, so there might be interesting (algebraic) properties when shuffling. From a computational point of view, 32 is smaller than 52, so my random generations will run faster. This is basically why I used 32. 10 would have been better, but not realistic with cards.

So, our 32 cards can be seen as a vector, or a list, of 32 items, say

$(a_1,a_2,\cdots, a_{31},a_{32})$

In order to assess if my cards are randomly sorted, let us get back to number properties (real valued numbers). If there were 10 cards, the list can be seen as an element of the following set

$(a_1,a_2,\cdots, a_{9},a_{10})\in\{0,1,\cdots,8,9\}^{10}$

(or to be more specific, a subset of that set, since numbers have to be different – it has to be a permutation – we cannot have duplicates, we’ll get back to that point in a few seconds). Let us see this list as a decimal number, with 10 digits. More precisely,

$u=\frac{a_1}{10}+\frac{a_2}{10^2}+\cdots+\frac{a_9}{10^9}+\frac{a_{10}}{10^{10}}\in[0,1]$

Now,  it is natural to say that cards are randomly sorted is this number is uniformly distributed on the unit interval, isn’t it ? (if we use the same shuffle many times, with the same starting point)

Well, if we think about it twice, uniform on the unit interval is probably not the proper distribution, since (as mentioned above) all digits have to be different. For instance, the smallest number would be $0.0123456789$ and the largest $0.9876543210$ . But as we will see, it this uniform assumption might not be too strong, actually.

And if we want to get back to our initial problem, with 32 cards, we simply have to use a decomposition in the 32-basis.

$u=\frac{a_1}{32}+\frac{a_2}{32^2}+\cdots+\frac{a_{31}}{32^{31}}+\frac{a_{32}}{32^{32}}\in[0,1]$

So if we have an algorithm to shuffle cards, we just have to run it several times (with the same starting value) and see when $u$ starts to be uniformly distributed. We start with a Dirac distribution, we have some kind of transition matrix, we expect our limiting distribution to be uniform and we wonder when the limiting distribution is reached… And from a statistical point of view, that should not be that difficult to assess, since we do have several goodness of fit tests that can be used.

Actually, it is possible to check if our technique passes the test of a uniform distribution, when digit are randomly generated (without replacement). The code to generate $u$ is

> j = 32
> X3 = (0:(j-1))[sample(1:j)]
> x3 = sum(j^(-(1:j))*X3)

If we run it a few times, and check if the assumption of a uniform distribution is valid (on samples with, say, 500 observations),

> P3=NULL
> for(i in 1:10000){
+   U3=NULL
+   for(s in 1:500){
+     X3 =(0:(j-1))[sample(1:j)]
+     x3 =sum(j^(-(1:j))*X3)
+     U3 =c(U3,x3)}
+   P3 =c(P3,ks.test(U3,punif)$p.value) + } in 95% of the scenarios, the $p$-value exceeds 5% > mean(P3>.05) [1] 0.9529 (which is something we should have under the null), More precisely, we can check that the $p$-value is uniformly distributed on the unit interval. > hist(P3,freq=FALSE) So assuming that our number is uniform on the unit interval might be a good notion for “cards are randomly sorted“. What we need now is some shuffling algorithms. Or to be more specific, some feasible shuffling algorithm. I mean here that I just start playing with cards, so it should be some techniques that I should be able to perform, to understand how it works…. So you will have to wait a few weeks before I start talking about the riffle or dovetail shuffle (you know the kind of shuffle in which half of the deck is held in each hand, and then cards are released by the thumbs so that they fall to the table interleaved… like in the movies) ! • Top in at random shuffle, and related (simple) algorithm My first algorithm is simple: the top-in at random shuffle. We start with the following ordering  N=1:m There are $m$ cards, and n denote the place where the card on top will go.  n=sample(2:m,size=1) if(n<m) N=c(N[2:n],N[1],N[(n+1):m]) if(n==m) N=c(N[2:n],N[1]) Then, we repeat that transfer of the card on top several times. schuffle1=function(m,ns=10){ N=1:m for(i in 1:ns) { n=sample(2:m,size=1) if(n<m) N=c(N[2:n],N[1],N[(n+1):m]) if(n==m) N=c(N[2:n],N[1]) } return(N)} Now, it is also possible to consider a bottom-in at random shuffle. The idea is the same, the only difference it that you start from the card at the bottom of the deck. But that would be the same as the one before (in terms of time before reaching randomness)  n=sample(1:(m-1),size=1) if(n>1) N=c(N[1:(n-1)],N[m],N[n:(m-1)]) if(n==1) N=c(N[m],N[1:(m-1)]) Why not mixing ? Randomly. Call it randomly mixed top-bottom in at random shuffle. You start either with the card on top, or at bottom (with identical probability), of the deck and then move the card somewhere,  card=sample(c("top","bottom"),size=1) if(card=="top"){ n=sample(2:m,size=1) if(n<m) N=c(N[2:n],N[1],N[(n+1):m]) if(n==m) N=c(N[2:n],N[1])} if(card=="bottom"){ n=sample(1:(m-1),size=1) if(n>1) N=c(N[1:(n-1)],N[m],N[n:(m-1)]) if(n==1) N=c(N[m],N[1:(m-1)])} All those codes can be together (within the same function), schuffle1=function(m,ns=10,which="top"){ N=1:m if(which=="top"){ for(i in 1:ns) { n=sample(2:m,size=1) if(n<m) N=c(N[2:n],N[1],N[(n+1):m]) if(n==m) N=c(N[2:n],N[1]) }} if(which=="bottom"){ for(i in 1:ns) { n=sample(1:(m-1),size=1) if(n>1) N=c(N[1:(n-1)],N[m],N[n:(m-1)]) if(n==1) N=c(N[m],N[1:(m-1)]) }} if(which=="mixed"){ for(i in 1:ns) {card=sample(c("top","bottom"),size=1) if(card=="top"){ n=sample(2:m,size=1) if(n<m) N=c(N[2:n],N[1],N[(n+1):m]) if(n==m) N=c(N[2:n],N[1]) } if(card=="bottom"){ n=sample(1:(m-1),size=1) if(n>1) N=c(N[1:(n-1)],N[m],N[n:(m-1)]) if(n==1) N=c(N[m],N[1:(m-1)]) } }} return(N)} But why do we take only one card ? It won’t be more complex to take 2. Or 3. Or more. • Tops in at random shuffle, and related (mixed) algorithm Yes, I used tops to say that we would take several cards on top of the deck. Say a random number of cards. And then, the strategy is the same, so the previous code is (slightly) adapted, as follows  k=sample(1:(m-1),size=1) n=sample((k+1):m,size=1); if(k==m-1) n=m if(n<m) N=c(N[(k+1):n],N[1:k],N[(n+1):m]) if(n==m) N=c(N[(k+1):n],N[1:k]) The idea is the following, here As earlier, it is possible to take cards at the bottom of the deck, or, one more time, to use a mixed strategy. The codes would be  card=sample(c("top","bottom"),size=1) if(card=="top"){ k=sample(1:(m-1),size=1) n=sample((k+1):m,size=1); if(k==m-1) n=m if(n<m) N=c(N[(k+1):n],N[1:k],N[(n+1):m]) if(n==m) N=c(N[(k+1):n],N[1:k])} if(card=="bottom"){ k=sample(2:m,size=1) n=sample(1:(k-1),size=1); if(k==1) n=1 if(n>1) N=c(N[1:(n-1)],N[k:m],N[n:(k-1)]) if(n==1) N=c(N[k:m],N[n:(k-1)])} Again, it is possible to have all those codes in the same function, schuffle2=function(m,ns=10,which="top"){ N=1:m if(which=="top"){ for(i in 1:ns) { k=sample(1:(m-1),size=1) n=sample((k+1):m,size=1); if(k==m-1) n=m if(n<m) N=c(N[(k+1):n],N[1:k],N[(n+1):m]) if(n==m) N=c(N[(k+1):n],N[1:k]) }} if(which=="bottom"){ for(i in 1:ns) { k=sample(2:m,size=1) n=sample(1:(k-1),size=1); if(k==1) n=1 if(n>1) N=c(N[1:(n-1)],N[k:m],N[n:(k-1)]) if(n==1) N=c(N[k:m],N[n:(k-1)]) }} if(which=="mixed"){ for(i in 1:ns) {card=sample(c("top","bottom"),size=1) if(card=="top"){ k=sample(1:(m-1),size=1) n=sample((k+1):m,size=1); if(k==m-1) n=m if(n<m) N=c(N[(k+1):n],N[1:k],N[(n+1):m]) if(n==m) N=c(N[(k+1):n],N[1:k]) } if(card=="bottom"){ k=sample(2:m,size=1) n=sample(1:(k-1),size=1); if(k==1) n=1 if(n>1) N=c(N[1:(n-1)],N[k:m],N[n:(k-1)]) if(n==1) N=c(N[k:m],N[n:(k-1)]) } }} return(N)} • How long should we shuffle before having cards randomly sorted ? With the codes mentioned above, it is possible to run generations of shuffles, distu=function(k=100,j=32){ U1B=U1T=U1M=U2B=U2T=U2M=U3=NULL for(s in 1:100){ X1T=(0:(j-1))[schuffle1(j,k,"top")] X1B=(0:(j-1))[schuffle1(j,k,"bottom")] X1M=(0:(j-1))[schuffle1(j,k,"mixed")] X2T=(0:(j-1))[schuffle2(j,k,"top")] X2B=(0:(j-1))[schuffle2(j,k,"bottom")] X2M=(0:(j-1))[schuffle2(j,k,"mixed")] X3 =(0:(j-1))[sample(1:j)] x1T=sum(j^(-(1:j))*X1T) x1B=sum(j^(-(1:j))*X1B) x1M=sum(j^(-(1:j))*X1M) x2T=sum(j^(-(1:j))*X2T) x2B=sum(j^(-(1:j))*X2B) x2M=sum(j^(-(1:j))*X2M) x3 =sum(j^(-(1:j))*X3) U1T=c(U1T,x1T) U1B=c(U1B,x1B) U1M=c(U1M,x1M) U2T=c(U2T,x2T) U2B=c(U2B,x2B) U2M=c(U2M,x2M) U3 =c(U3,x3) } B=list(U1T=U1T,...) } and then, we run tests to see if the samples can be assumed to be uniformly distributed on the unit interval, e.g. for the very first kind first shuffle describe above, it would be ks.test(B$U1T,punif)$p.value More precisely, we use the following function, to estimate to proportion of scenarios where the $p$-value exceeds 5%, PV=function(k){ P1B=P1T=P1M=P2B=P2T=P2M=P3=NULL for(i in 1:10000){ B=dist(k,j=32) P1T=c(P1T,ks.test(B$U1T,punif)$p.value) P1M=c(P1M,ks.test(B$U1M,punif)$p.value) P2T=c(P2T,ks.test(B$U2T,punif)$p.value) P2M=c(P2M,ks.test(B$U2M,punif)$p.value) P3 =c(P3,ks.test(B$U3,punif)\$p.value)}
return(list(
p1T=mean(P1T>.05),
p1M=mean(P1M>.05),
p2T=mean(P2T>.05),
p2M=mean(P2M>.05),
p3=mean(P3>.05)))}

If we plot the results, we have

K=1:100
MP=Vectorize(PV)(K)
plot(K,MP[1,],col="red",type="b",ylim=0:1,pch=1)
lines(K,MP[2,],type="b",pch=19,col="red")
lines(K,MP[3,],col="blue",type="b",pch=1)
lines(K,MP[4,],type="b",pch=19,col="blue")
lines(K,MP[5,],type="b",pch=3,col="black")

Here, we look at the proportion of $p$-values that exceed 5%. We can pretend that we have a uniform distribution if that proportion is close to 95%. So basically, we just have to see when we reached for the first time the 95% region.If we zoom in the upper part of the graph, we get

With 32 cards,

• with a top in at random, we have to shuffle about 70 or 80 cards before having a randomly sorted set of cards. Which is large, but which is quite intuitive. One can imagine that it might take a while before getting the cards at bottom much higher in the pack,
• with a randomly mixed top in at random strategy, it is faster, slightly (we do not have that problem with cards at bottom that stay at bottom), since it takes about 60 or 70 cards.
• with a tops in a random, it goes again faster, with about 35 rounds,
• with a randomly mixed tops-bottoms in at random, it takes about 10 to 15 rounds.

Those results were obtained on tests on samples of size 100. The same code ran on a server during the week-end, with samples of size 500. Note that the output is rather close,

Note that those algorithm were mentioned because they were feasible, not only from a computational point of view, but when playing with real cards, in paper. Like with kids. I can actually ask my kids to perform those shuffle techniques next time we play with cards. The good thing is that randomly mixed tops-bottoms in at random shuffle technique: kids can do it 10 times, and cards should be randomly ordered in the deck…

Now, for those willing to see more algorithms, there are the so-called Fisher-Yates (also Knuth) shuffle. But may I keep that for another post ?

# Fisher-Tippett theorem and limiting distribution for the maximum

Tomorrow, we will discuss Fisher-Tippett theorem. The idea is that there are only three possible limiting distributions for normalized versions of the maxima of i.i.d. samples . For bounded distribution, consider e.g. the uniform distribution on the unit interval, i.e. on the unit interval. Let and . Then, for all and ,

i.e. the limiting distribution of the maximum is Weibull’s.

set.seed(1)
s=1000000
n=100
M=matrix(runif(s),n,s/n)
V=apply(M,2,max)
bn=1
an=1/n
U=(V-bn)/an
hist(U,probability=TRUE,,col="light green",
xlim=c(-7,1),main="",breaks=seq(-20,10,by=.25))
u=seq(-10,0,by=.1)
v=exp(u)
lines(u,v,lwd=3,col="red")

For heavy tailed distribution, or Pareto-type tails, consider Pareto samples, with distribution function . Let and , then

which means that the limiting distribution is Fréchet’s.

set.seed(1)
s=1000000
n=100
M=matrix((runif(s))^(-1/2),n,s/n)
V=apply(M,2,max)
bn=0
an=n^(1/2)
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(0,7),main="",breaks=seq(0,max(U)+1,by=.25))
u=seq(0,10,by=.1)
v=dfrechet(u,shape=2)
lines(u,v,lwd=3,col="red")

For light tailed distribution, or exponential tails, consider e.g. a sample of exponentially distribution variates, with common distribution function . Let and , then

i.e. the limiting distribution for the maximum is Gumbel’s distribution.

library(evd)
set.seed(1)
s=1000000
n=100
M=matrix(rexp(s,1),n,s/n)
V=apply(M,2,max)
(bn=qexp(1-1/n))
log(n)
an=1
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Consider now a Gaussian sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)

as . Let and . Then we can get

as . I.e. the limiting distribution of the maximum of a Gaussian sample is Gumbel’s. But what we do not see here is that for a Gaussian sample, the convergence is extremely slow, i.e., with 100 observations, we are still far away from Gumbel distribution,

and it is only slightly better with 1,000 observations,

set.seed(1)
s=10000000
n=1000
M=matrix(rnorm(s,0,1),n,s/n)
V=apply(M,2,max)
(bn=qnorm(1-1/n,0,1))
an=1/bn
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Even worst, consider lognormal observations. In that case, recall that if we consider (increasing) transformation of variates, we are in the same domain of attraction. Hence, since , if

then

i.e. using Taylor’s approximation on the right term,

This gives us normalizing coefficients we should use here.

set.seed(1)
s=10000000
n=1000
M=matrix(rlnorm(s,0,1),n,s/n)
V=apply(M,2,max)
bn=exp(qnorm(1-1/n,0,1))
an=exp(qnorm(1-1/n,0,1))/(qnorm(1-1/n,0,1))
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,40,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")