Yesterday, I was uploading some old posts to complete the migration (I get back to my old posts, one by one, to check links of pictures, reformating R codes, etc). And I re-discovered a post published amost 2 years ago, on nuns and Hell’s Angels in an airplaine.

It reminded me an old probability problem (that might be known as one on Feymann’s problems):* suppose that you have a prescription to take half pills for 6 days. Unfortunately the pharmacist was a bit lazy (or just wanted to help me to write a mathematical problem), and he gives 3 (full) pills in a small box. Day 1, you take a pill, break it in two parts, eat one, and return the other half in the box. Day 2, you draw randomly ‘something’ from the box, i.e. either half a pill, or a pill. If it’s a half one, then you eat it. If it is a fill one, you break it in two, eat one half, and return the other half in the box. Etc.On Day 6, if my story was well explained, you should know that there can only be one half pill. So far, so good. But what about Day 5 ? There were either two half pills, or one full pill. But what was the probability that there was a fill pill in the box on Day 5 ?*

Nice problem, isn’t it ?

The good thing is that it can be modeled as a Markovian model. Assume that we do have pills. After days, the box will be empty. Consider the pair denoting the number of half pills, and complete pills. can take all values, from 0 to , and will be positive, with . Thus, the number of states – possible pairs from Day 1 till Day – will be , i.e. . More precisely, define those states in a dataframe,

> n=3 > COMPLETE=HALF=NULL > for(i in n:0){ + HALF=c(0:(n-i),HALF) + COMPLETE=c(rep(i,length(0:(n-i))),COMPLETE) + } > k=length(COMPLETE) > state=data.frame(s=1:k,nc=rev(COMPLETE),nh=rev(HALF)) > state s nc nh 1 1 3 0 2 2 2 1 3 3 2 0 4 4 1 2 5 5 1 1 6 6 1 0 7 7 0 3 8 8 0 2 9 9 0 1 10 10 0 0

Now, we can play to derive the transition matrix of the Markov chain.

> attach(state) > P=matrix(0,k,k) > for(i in 1:k){ + C=state$nc[i] + H=state$nh[i] + if((C>0)&(H>0)){ + P[i,state[(nc==C-1)&(nh==H+1),"s"]]= C/(C+H) + P[i,state[(nc==C)&(nh==H-1),"s"]]= H/(C+H)} + if((C>0)&(H==0)){ + P[i,state[(nc==C-1)&(nh==H+1),"s"]]=1} + if((C==0)&(H>0)){ + P[i,state[(nc==C)&(nh==H-1),"s"]]=1} + if((C==0)&(H==0)){ + P[i,state[(nc==C)&(nh==H),"s"]]=1} + }

We do have a transition matrix (or a probability matrix) since all elements are positive, and the sum per line is 1,

> apply(P,1,sum) [1] 1 1 1 1 1 1 1 1 1 1

Here, the transition matrix is the following

> P [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 1 0.00 0.00 0.00 0.0 0.00 0.0 0 0 [2,] 0 0 0.33 0.66 0.00 0.0 0.00 0.0 0 0 [3,] 0 0 0.00 0.00 1.00 0.0 0.00 0.0 0 0 [4,] 0 0 0.00 0.00 0.66 0.0 0.33 0.0 0 0 [5,] 0 0 0.00 0.00 0.00 0.5 0.00 0.5 0 0 [6,] 0 0 0.00 0.00 0.00 0.0 0.00 0.0 1 0 [7,] 0 0 0.00 0.00 0.00 0.0 0.00 1.0 0 0 [8,] 0 0 0.00 0.00 0.00 0.0 0.00 0.0 1 0 [9,] 0 0 0.00 0.00 0.00 0.0 0.00 0.0 0 1 [10,] 0 0 0.00 0.00 0.00 0.0 0.00 0.0 0 1

In order to get our probability, let us start from state 1 – or – with probability 1, and let us look at the distribution at different periods,

> dist=c(1,rep(0,k-1)) > MatDist=matrix(NA,2*n+1,k) > MatDist[1,]=dist > for(i in 1:(2*n)){dist=as.vector(t(dist)%*%P) + MatDist[i+1,]=dist + }

(one can check that after days, the box is empty). The probability is given in row , and we just have to check which column corresponds to the pair ,

> vs=state[which(MatDist[2*n-1,]>0),] > proba=MatDist[2*n-1,vs[vs$nc==1,"s"]] > proba [1] 0.3888889

Here the probability of having a full pair on Day 5 is 38.89%.

Actually, it is possible to study the evolution of this probability as a function of ,

> computeproba=function(n=3){ + COMPLETE=HALF=NULL + for(i in n:0){ + HALF=c(0:(n-i),HALF) + COMPLETE=c(rep(i,length(0:(n-i))),COMPLETE) + } + k=length(COMPLETE) + state=data.frame(s=1:k,nc=rev(COMPLETE),nh=rev(HALF)) + P=matrix(0,k,k) + for(i in 1:k){ + C=state$nc[i] + H=state$nh[i] + if((C>0)&(H>0)){ + P[i,state[(state$nc==C-1)&(state$nh==H+1),"s"]]= C/(C+H) + P[i,state[(state$nc==C)&(state$nh==H-1),"s"]]= H/(C+H)} + if((C>0)&(H==0)){ + P[i,state[(state$nc==C-1)&(state$nh==H+1),"s"]]=1} + if((C==0)&(H>0)){ + P[i,state[(state$nc==C)&(state$nh==H-1),"s"]]=1} + if((C==0)&(H==0)){ + P[i,state[(state$nc==C)&(state$nh==H),"s"]]=1} + } + dist=c(1,rep(0,k-1)) + MatDist=matrix(NA,2*n+1,k) + MatDist[1,]=dist + for(i in 1:(2*n)){dist=as.vector(t(dist)%*%P) + MatDist[i+1,]=dist + } + vs=state[which(MatDist[2*n-1,]>0),] + proba=MatDist[2*n-1,vs[vs$nc==1,"s"]] + return(proba) + }

If we plot the probability as a function of , we get

> P=Vectorize(computeproba)(2:40) > plot(2:40,P,ylim=c(0,.5))

One can observe that the probability is decreasing. But slowly, extremely slowly. With a log scale on the y-axis, we have

> plot(2:40,P,ylim=c(0,.5),log="y")

If we look for ‘high’ values, we can get

> computeproba(100) [1] 0.14218

I do not know if this limit goes to 0 as goes to infinity. Actually, since we do have to compute a matrix with entries i.e. roughly , cannot be that large… Too bad. If anyone knows how this probability behaves as a function of , when is large, I’d be glad to know…

For N pills that probability that there is a full pill in the penultimate day is 1 / (2*N – 1).

What the distribution of full pills to half pills looks like over the intermediate day is certainly interesting though

Doggone it! Christian beat me to it . The approach you took(I think) became difficult because it solves for the available “choices” at every day rather than just at the penultimate day.

Hi Arthur!

Think of this as coloured half-pills. The probably of having one intact pill in the box on day 5 is the same as that of having two half-pills of the same colour. Thinking about it this way, you don’t have to worry about when a pill is split in two halves. Then you can solve it using the multinomial theorem (with the colours as “boxes”).