sharnett

Interview with Lugg (YC S15) with some details on how they did it: https://www.ycombinator.com/blog/interview-with-lugg-founder...

Very interesting! I explored this method and compared it to Newton's method, you may like it: https://github.com/sharnett/point_to_ellipse/blob/master/poi...

"x is with a likelihood of 95% between a and b" -- that is not what a confidence interval means. Confidence intervals are about as confusing as p-values and don't really solve the problem.

Very nice, much faster than simulating it. I guess if you're using an informative prior, instead of adding 1 to the number of successes and failures, you add the corresponding parameters of your beta prior?

A pretty good shortcut (if you don't have a log-beta function, for example) is to approximate both A and B with the normal distribution. Then their difference is also normal, so you can just check the probability that a normal random variable is greater than zero.

Specifically, the mean μ of a beta distribution is α/(α+β) and the variance σ^2 is αβ/((α+β)^2(α+β+1)). Use these as the parameters for your normal approximation, and we have the difference D ~ N(μ_A-μ_B, σ_A^2+σ_B^2). The probability that B beats A is just the CDF of D evaluated at 0.

In Python:

  from scipy.stats import norm as norm

  def beta_mean(a, b):
      return a/(a+b)
  def beta_var(a, b):
      return a*b/((a+b)**2*(a+b+1))
  def probability_B_beats_A(α_A, β_A, α_B, β_B):
      mu = beta_mean(α_A, β_A) - beta_mean(α_B, β_B)
      sigma = (beta_var(α_A, β_A) + beta_var(α_B, β_B))**.5
      return norm.cdf(0, mu, sigma)

A nit pick: I believe your 'readable' interpretation of the confidence interval isn't quite right: http://en.wikipedia.org/wiki/Confidence_interval#Meaning_and...