How to verify Bayes theorem output?

Question

I am looking for if I can do some simulation in R for verifying the output of Bayes theorem , if you have some other technique I am open to know that also. Anything that can help me in verification would be useful.

Actually my problem is something like this:-

I have probabilities of dropping a call based on the duration of time customer has waited for before dropping.

Probability of dropping during 0-20sec: $0.23$

Probability of dropping during 20-40sec: $0.11$

Probability of dropping during 40-50sec: $0.04$

Probability of dropping during 50-60sec: $0.61$

Now the probability of dropping the call during 50-60sec given that the call has not been dropped during 0-20sec is denoted as P(50-60sec dropped/Not dropped 0-20sec).

And P(50-60sec dropped/Not dropped 0-20sec) is found to be $\frac{0.61}{1-0.23}=0.79 $.

So now I want to verify if 0.79 is correct or not by some simulation or any other technique.

EDIT: The sum of probabilities adds up to 0.99 because those are calculated from a dataset and rounded off up to 2 digits.

Answer 1

Your first problem is that your probabilities don't sum to 1.

Here's some Python that generates 100,000 samples. I'll leave it to you to write it in R.

import numpy as np

n50_60 = n_not0_20 = 0

for i in range(100_000):
    u = np.random.uniform()
    if u <= 0.23:
        pass
    elif u <= 0.23 + 0.11:
        n_not0_20 += 1
    elif u <= 0.23 + 0.11 + 0.04:
        n_not0_20 += 1
    else:
        n_not0_20 += 1
        n50_60 += 1

print(n50_60 / n_not0_20)
```

Answer 2

Here is the simulation in R which tries to verify the Probability of dropping during 50-60sec by using P(50-60sec dropped/Not dropped 0-20sec).

#This will have ones with the probability of 0.77
not_20_sec<-rbinom(n=10000, size=1, prob=1-0.23)
#Below code will put 'ones'1' in 'y' with the probabilty of 0.7951 
#wherever it was a '1' in  not_20_sec else it would add '0'.
y=c()
for (i in not_20_sec){
  if (i==1){
    y<-c(y,rbinom(n=10, size=1, prob=0.7951))
  }
  else{
    y<-c(y,replicate(10,0))
  }
}
#This is the probability of dropping during 50-60sec.
print(sum(y==1)/length(y)) # 0.609