DZone
Thanks for visiting DZone today,
Edit Profile
  • Manage Email Subscriptions
  • How to Post to DZone
  • Article Submission Guidelines
Sign Out View Profile
  • Post an Article
  • Manage My Drafts
Over 2 million developers have joined DZone.
Log In / Join
Refcards Trend Reports Events Over 2 million developers have joined DZone. Join Today! Thanks for visiting DZone today,
Edit Profile Manage Email Subscriptions Moderation Admin Console How to Post to DZone Article Submission Guidelines
View Profile
Sign Out
Refcards
Trend Reports
Events
Zones
Culture and Methodologies Agile Career Development Methodologies Team Management
Data Engineering AI/ML Big Data Data Databases IoT
Software Design and Architecture Cloud Architecture Containers Integration Microservices Performance Security
Coding Frameworks Java JavaScript Languages Tools
Testing, Deployment, and Maintenance Deployment DevOps and CI/CD Maintenance Monitoring and Observability Testing, Tools, and Frameworks
Culture and Methodologies
Agile Career Development Methodologies Team Management
Data Engineering
AI/ML Big Data Data Databases IoT
Software Design and Architecture
Cloud Architecture Containers Integration Microservices Performance Security
Coding
Frameworks Java JavaScript Languages Tools
Testing, Deployment, and Maintenance
Deployment DevOps and CI/CD Maintenance Monitoring and Observability Testing, Tools, and Frameworks
What's in store for DevOps in 2023? Hear from the experts in our "DZone 2023 Preview: DevOps Edition" on Fri, Jan 27!
Save your seat
  1. DZone
  2. Data Engineering
  3. Databases
  4. Pills, Half Pills, and Probabilities

Pills, Half Pills, and Probabilities

Arthur Charpentier user avatar by
Arthur Charpentier
·
Feb. 12, 13 · Interview
Like (0)
Save
Tweet
Share
2.29K Views

Join the DZone community and get the full member experience.

Join For Free

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 , 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…

Matrix (protocol) POST (HTTP) R (programming language) Column (database) Links IT Row (database) Distribution (differential geometry)

Published at DZone with permission of Arthur Charpentier, DZone MVB. See the original article here.

Opinions expressed by DZone contributors are their own.

Popular on DZone

  • What Java Version Are You Running? Let’s Take a Look Under the Hood of the JDK!
  • Simulate Network Latency and Packet Drop In Linux
  • Spring Cloud: How To Deal With Microservice Configuration (Part 1)
  • Unleashing the Power of JavaScript Modules: A Beginner’s Guide

Comments

Partner Resources

X

ABOUT US

  • About DZone
  • Send feedback
  • Careers
  • Sitemap

ADVERTISE

  • Advertise with DZone

CONTRIBUTE ON DZONE

  • Article Submission Guidelines
  • Become a Contributor
  • Visit the Writers' Zone

LEGAL

  • Terms of Service
  • Privacy Policy

CONTACT US

  • 600 Park Offices Drive
  • Suite 300
  • Durham, NC 27709
  • support@dzone.com
  • +1 (919) 678-0300

Let's be friends: