You're definitely messing something up somewhere in that analysis. Where's that numpy/betadistribution snippet coming from?
Back to basics: If you assume that each mile is uncorrelated, (on the large not a crazy assumption), then the overall number of accidents over N miles can be modeled by a binomial distribution with N trials and some crash probability p. You're trying to use bayesian analysis to estimate that p.
The core concept there is that the probability of the observation given a certain p is directly proportional to the probability of the p given a certain observation.
Given p, it's trivial to compute the chance of the observation; just plug it into the pmf. For exactly one crash (and a small p), you can even quite accurately approximate it by Np/exp(Np) due to the trivial binomial coefficient for k=1 and the taylor expansion of the base-near-1 power.
So regardless of the numerical analysis that lead you to believe otherwise: you can just plug in numbers can get a numeric approximation of the (scaled) distribution of the probability mass for p. When I do that, both with an exponential prior (i.e. p in 0.1-0.01 is as likely as p in 0.01-0.001) or a uniform prior (i.e. p in 0.1-0.2 is as likely as p in 0.2-0.3) I get the unsurprising result that the maximum likelihood estimation for p is 1/N, and that the median p is a little lower; at around 1 accident in 1.5e8 miles
In practice that means that common sense is correct: if you observe 1 accident in N miles, a reasonable approximation of the best guess of the true accident rate per mile is 1/N.
Now, you can do this analytically, and the conjugate prior of the binomial is indeed a beta distribution. But that forces you to pick parameters for the prior from that beta distribution, and it's easy to do that wrongly.
A reasonable prior might be Beta[1,3] which assumes p is likely closer to 0 than 1. Analytically [then](http://www.johndcook.com/blog/conjugate_prior_diagram/), the posterior distribution is Beta[1+1,3+1e8-1]; with a mean p of 2e-8 (but note that the details of the chosen parameters definitely impact the outcome). This is quite close to the numerically derived p, though probably less accurate since it is constrained to a unreasonable prior.
So I'm not sure exactly what your python script is computing, but commonsense, numerical analysis, and Bayesian analysis all arrive at roughly 1/N in my case - you probably made some mistake in your reasoning somewhere (if you elaborate in detail how you arrive at this numpy snippet, maybe we'll find the error).
Whoops, you are right, my posterior should have been Beta[2,100e6+1] (I chose a uniform prior). With that I get a median probability of 1.7e-8, and a mean probability of 2e-8. Good catch!
Note that this is not particularly sensitive to whether you choose beta[1,1] or Beta[1,3] as the prior.
Yeah, that makes virtually no difference. I'm not familiar with picking these priors, so I could imagine picking (say) Beta[0.5,1] or Beta[5,30] as prior too, and that does make some difference (as in the prior alpha has some impact).
Is there some principled reason to pick one over the other, or is it just "whatever looks reasonable"?
I also can't explain why plain numerical simulation doesn't come up with almost identical p values - I get 1/1.5e8 to 1/1.7e8, but that's 6.25e-9. I mean, it's the same order of magnitude, but it's nevertheless quite different.
Oh well, unless you've got some insight there, I guess I'm going to give up on that analysis - that kind of nitty gritty stuff is a timesink ;-).
Back to basics: If you assume that each mile is uncorrelated, (on the large not a crazy assumption), then the overall number of accidents over N miles can be modeled by a binomial distribution with N trials and some crash probability p. You're trying to use bayesian analysis to estimate that p.
The core concept there is that the probability of the observation given a certain p is directly proportional to the probability of the p given a certain observation.
Given p, it's trivial to compute the chance of the observation; just plug it into the pmf. For exactly one crash (and a small p), you can even quite accurately approximate it by Np/exp(Np) due to the trivial binomial coefficient for k=1 and the taylor expansion of the base-near-1 power.
So regardless of the numerical analysis that lead you to believe otherwise: you can just plug in numbers can get a numeric approximation of the (scaled) distribution of the probability mass for p. When I do that, both with an exponential prior (i.e. p in 0.1-0.01 is as likely as p in 0.01-0.001) or a uniform prior (i.e. p in 0.1-0.2 is as likely as p in 0.2-0.3) I get the unsurprising result that the maximum likelihood estimation for p is 1/N, and that the median p is a little lower; at around 1 accident in 1.5e8 miles
In practice that means that common sense is correct: if you observe 1 accident in N miles, a reasonable approximation of the best guess of the true accident rate per mile is 1/N.
Now, you can do this analytically, and the conjugate prior of the binomial is indeed a beta distribution. But that forces you to pick parameters for the prior from that beta distribution, and it's easy to do that wrongly.
A reasonable prior might be Beta[1,3] which assumes p is likely closer to 0 than 1. Analytically [then](http://www.johndcook.com/blog/conjugate_prior_diagram/), the posterior distribution is Beta[1+1,3+1e8-1]; with a mean p of 2e-8 (but note that the details of the chosen parameters definitely impact the outcome). This is quite close to the numerically derived p, though probably less accurate since it is constrained to a unreasonable prior.
So I'm not sure exactly what your python script is computing, but commonsense, numerical analysis, and Bayesian analysis all arrive at roughly 1/N in my case - you probably made some mistake in your reasoning somewhere (if you elaborate in detail how you arrive at this numpy snippet, maybe we'll find the error).