# Coin with an unknown weight

#### AngleWyrm

##### Active Member
Was watching this video, supposedly part two of a three-part series, but the author never did get around to making the final part which he advertised as answering the question set up in part one. The question asked "2 defects found in 100 samples; what can we say about the probability of a defect?"

He got as far as "What's the probability density function that describes the value h after seeing a few outcomes?" but that was two years ago.

#### Dason

Was he a frequentist? If so he probably was going to make one of the many different options for confidence intervals.

Was he Bayesian? If so they probably just looked at the posterior distribution to get an interval.

#### AngleWyrm

##### Active Member
This is the screen presented in the video:

His claim is there's a way to find the PDF from the sample set, and it can answer questions such as what's the probability that the P(heads) is in this range.

Does that make him a Frequentist or a Bayesian?

#### Dason

If they're putting a statement like P(0.6 <= P(H) <= 0.8) into the video then they're definitely either a Bayesian or a Frequentist that is very very sloppy or undereducated.

My guess is Bayesian though. They're probably just using a uniform prior so the posterior will just be a beta distribution.

#### hlsmith

##### Less is more. Stay pure. Stay poor.
So of all the realizations that could have happened, how many times were these within the below range / sum of all possibilities.

Edit: I got 57% of values landing in this range, but I could have buggered this a little.

Code:
# define grid
p_grid <- seq( from=0 , to=1 , length.out=100 )
# define prior
prior <- rep( 1 , 100 )
# compute likelihood at each value in grid
likelihood <- dbinom( 7 , size=10 , prob=p_grid )
# compute product of likelihood and prior
unstd.posterior <- likelihood * prior
# standardize the posterior, so it sums to 1
posterior <- unstd.posterior / sum(unstd.posterior)
plot( p_grid , posterior , type="b" ,
xlab="probability of Heads" , ylab="posterior probability")
abline(v = c(0.6,0.8), col = "red", lty = 3)
# Calculate the probability
sum(posterior[60:80])

#### Dason

So of all the realizations that could have happened, how many times were these within the below range / sum of all possibilities.
View attachment 3956
Edit: I got 57% of values landing in this range, but I could have buggered this a little.

Code:
# define grid
p_grid <- seq( from=0 , to=1 , length.out=100 )
# define prior
prior <- rep( 1 , 100 )
# compute likelihood at each value in grid
likelihood <- dbinom( 7 , size=10 , prob=p_grid )
# compute product of likelihood and prior
unstd.posterior <- likelihood * prior
# standardize the posterior, so it sums to 1
posterior <- unstd.posterior / sum(unstd.posterior)
plot( p_grid , posterior , type="b" ,
xlab="probability of Heads" , ylab="posterior probability")
abline(v = c(0.6,0.8), col = "red", lty = 3)
# Calculate the probability
sum(posterior[60:80])
A grid approach isn't awful but in this particular case we do have the exact posterior distribution known.
Code:
> diff(pbeta(c(.6, .8), 1+7, 1+3))
[1] 0.5425765
Which is somewhat close to what you had. But your code a few issues - 1) why only 100 sampling points? They weren't even on the 100ths spots. The other is that you used indices 60-80 but those weren't exactly what you wantedv since that actually gave you some spots less than .6
Code:
> p_grid <- seq( from=0 , to=1 , length.out=100 )
> p_grid[60:80]
[1] 0.5959596 0.6060606 0.6161616 0.6262626 0.6363636 0.6464646 0.6565657
[8] 0.6666667 0.6767677 0.6868687 0.6969697 0.7070707 0.7171717 0.7272727
[15] 0.7373737 0.7474747 0.7575758 0.7676768 0.7777778 0.7878788 0.7979798
Let's increase the amount of sampling points and fix the indexing
Code:
p_grid <- seq( from=0 , to=1 , length.out=1000 )
# define prior
prior <- rep( 1 , 1000 )
# compute likelihood at each value in grid
likelihood <- dbinom( 7 , size=10 , prob=p_grid )
# compute product of likelihood and prior
unstd.posterior <- likelihood * prior
# standardize the posterior, so it sums to 1
posterior <- unstd.posterior / sum(unstd.posterior)
plot( p_grid , posterior , type="b" ,
xlab="probability of Heads" , ylab="posterior probability")
abline(v = c(0.6,0.8), col = "red", lty = 3)
# Calculate the probability
sum(posterior[which(p_grid >= .6 & p_grid <= .8)])
which gives
Code:
> sum(posterior[which(p_grid >= .6 & p_grid <= .8)])
[1] 0.5430052
Which is much closer to the analytic solution.

#### hlsmith

##### Less is more. Stay pure. Stay poor.
I had thought about increasing the sample to something larger. Which I would normally do in a similar setting or for bootstrapping to get at percentage values like 97.5%, which you can't when only using 100 values.

Good point on my use of sum(posterior[60:80]), I had wanted to do it your way but was too lazy and didn't think hard enough to realize a couple values outside the bounds could have snuck in.

#### hlsmith

##### Less is more. Stay pure. Stay poor.
Side note, people like Sander Greenland have been pushing for pvalue functions, which can be used to get a density looking plot like used in Bayesian posteriors, but the issue is their interpretation would still be that of a null hypothesis.

#### AngleWyrm

##### Active Member
Well here's the scenario I'm working with: There's a gambler's box with a lever, pull the lever and you could win a prize. The creator of the box has been messing around with many cogs and gears and pulleys and such, but says "Trust me, it's 50% chance to win!"

What I want to know is at what point can I reasonably declare the box is fair?
For any given number of pulls on the lever how confident can I be the result is drawn from the stated population of 1/2 win 1/2 lose?

#### hlsmith

##### Less is more. Stay pure. Stay poor.
Well you have to define a threshold for 'reasonably declare the box is fair". This is likely subjective or based on an acceptable loss level.

#### AngleWyrm

##### Active Member
See part two of the question:
For any given number of pulls on the lever how confident can I be the result is drawn from the stated population of 1/2 win 1/2 lose?