Why Bayesian Stats Needs Monte-Carlo Methods
countbayesie.com
countbayesie.com
Analytic methods are really useful when they can be done, but often have to be worked out for each specific model.
1. Why do we need Bayesian estimation? Can't we just all use MLE, and live in harmony? If MLE was good for Gauss and Laplace and Fisher, why isn't it good for me? Answer: for many problems the likelihood function is unbound. In fact the cases where the likelihood function are bounded, so the MLE makes sense, are few and far between. But if you think of the likelihood as a function with a hump somewhere, you have several choices for a single number that best describes the function: mode, mean, median for example. MLE corresponds to mode, can we use the mean instead? Also, does it not even make better sense to use mean rather than mode? Sure it does, and if you use mean you become a Bayesian. But can the likelihood be unbounded and have a finite mean, you ask? Yes, it happens quite often, sufficiently often that in practice people don't worry about that.
Once you start using mean instead of mode (and therefore you can proudly call yourself a Bayesian), you get a number of side benefits: you get parameter uncertainty estimation for free, and you can start multiplying your likelihoods with various functions that make your life easier. You can interpret those functions as regularization functions, or prior, you can even start talking about expert knowledge, and delve into philosophical meanderings. No need to do that , there are enough people on the internet who do that.
2. Why Monte Carlo? Because calculating means of high-dimensional distributions is not easy. In most cases you are hit with the curse of dimensionality. Separately, the distributions that you work with are not known only up to a multiplicative constant. So the problem is how do you estimate the mean of a vector (v1, v2, ...., vn) when you know its distribution, but only up to a constant? The answer is Markov Chain Monte Carlo (MCMC) and if you learn that, you get special powers.
So why do Bayesians need Monte Carlo? Because they like to be like super-heroes and have special powers.
- Why Engineering needs PDE solvers.
- Why Computer Science needs compilers.
- Why Geography needs digital maps.
PS: I think it would be great to have a clickbait title builder based on top HN article
And one for generating (unrelated) comments to a title: https://github.com/leod/hncynic
The only sense I can make of this is that:
- Predictions about stochastic dynamic systems are hard, especially when there can be these exogenous variables intervening from out of nowhere.
- Competitive situations are hard to predict, especially if they change based on the measurements produced about them (i.e. the campaigns strategize based on polling info). This effectively makes those measurements less predictive.
To me, this suggests that all the polling and analysis are a waste of resources and attention. The expected value of information is extremely low. We should kick off campaigns on Halloween, talk about them for a couple days, and vote immediately. All the pollsters and analysts can be more usefully deployed towards studying other systems.
It seems more likely that the media companies profit enormously from all of this coverage in spite of (and perhaps even because of) the low certainty provided by this polling and analysis.
In industry, you just won’t find this. Statistics is a small domain specialist component of a much larger software ecosystem. You don’t require advanced statistics, you just might want to try it or derive a special benefit from it, and none of that will be worth it unless it efficiently connects to the regular software development life cycle, architecture patterns, testing, and maintenance processes that the larger system requires.
Someone who can do “a little software development” would be a liability in that scenario, no matter what their other domain specializations are.
That’s why you’ll find most data scientists and ML engineers professionally are highly skilled at software engineering. They approach statistical modeling tasks from the point of view of standard software development life cycle tasks where the components just happen to involve statistical models or techniques.
People that have this skillset (knowledge of Monte Carlo methods for statistical analysis/verification) are useful in this sphere.
An example problem we deal with is as follows:
We have some arbitrary histogram (list of slot game awards and probabilities), and we wish to map this to a particular game of bingo (N ball calls out of M balls in total, 5x5 bingo card). The constraint here is that when you come up with a list of bingo patterns for awards, they are ordered in the sense that when evaluating the set of bingo patterns, from top-down, the first bingo pattern that matches is awarded and the evaluation stops.
Getting the optimal fit theoretically in the general case requires exponential time, as the chain of P(this pattern given (not-that pattern and not-other pattern and not-...)) grows with an exponential number of operations in the theoretical expansion. Each time a new pattern is fit to an element of the histogram, the probabilities of all remaining possible patterns to fit changes.
There are a few ways of tackling this problem in a more optimal way: If you can solve this with good enough accuracy in linear time with a good math/cs background, you’d be useful to any company playing in this space.
> this is an awesome example. is there an easy way (as in non-brute force) to finding the beta parameters that'll match a probability?
It makes sense, I'm just surprised the blog never comes out and says "you have to do a brute force search to get the parameters for the likelihoods A ~ Beta(2,13) and B ~ Beta(3,11) for P(B > A) to have about the right value".
Further what do you mean by
> “ the blog never comes out and says "you have to do a brute force search to get the parameters for the likelihoods A ~ Beta(2,13) and B ~ Beta(3,11) for P(B > A) to have about the right value".”
What do you mean by “brute force search” here?
He also explains that easy analytical solutions often don't exists to seemingly easy questions, such as this one. That's why Bayesian statistics needs Monte-Carlo.
What I'm asking is: how do you find parameters for A and B given "P(B > A) = 0.71"?
About the actual question, I don't think that's very well specified in this form. You could have A~Beta0 where Beta0 is very concentrated around let's say 0.5, and then find a B~Beta1 where Beta1 has 71% of its weight above 0.5. That would be a pretty good approximation for it, in fact as close as you want if you increase sample size of A. I can't verify this, but I think finding Beta1 should be easy. Obv this is not what you are looking for though, because this would require A to have a much larger sample size than B. I guess we should add some restriction, like similar sample sizes.
from scipy.stats import beta
import numpy as np
def betamax(params0, params1, integration_points=5000):
u0 = params0[0] / sum(params0)
u1 = params1[0] / sum(params1)
flip = u0 < u1
if flip:
params0, params1 = params1, params0
q = beta(*params0).ppf(np.linspace(0, 1, integration_points + 2)[1:-1])
c = beta(*params1).cdf(q)
p = c.mean()
if flip:
return 1 - p
return p >>> betamax((3, 14), (4, 12))
0.29448379324802676
>>> betamaxinv(0.29, allowed_error=0.01)
[((3, 14), (4, 12)), ...]
I agree that 71% is a low amount of evidence, but election forecasts really make my brain work to rationalize them. I think it's helpful to remember that forecasts are produced _after_ weighting them by states' electoral college shares.
So when someone analysing an election comes up with a 70/30 split, the number seems incredibly high - there is now way we are getting this one wrong. It takes at least a second look to see that this guy is talking about something different and that the race actually is quite close this time.
I think the confusing thing here is that at first glance it looks like "Joe Biden will get 71% of the votes" which is a landslide victory. If it said "there is as much chance of Joe Biden winning as rolling a 1-7 on a 10 sided dice", which means almost exactly the same thing, no one would get surprised if Trump won, because we are all familiar with dice rolls with probabilities in the 0.1-0.9 range.
https://math.stackexchange.com/questions/1635949/is-there-a-...
The use of 'would' implies that automatic differentiation has not been established. Also, if something is 40+ years old, then it follows that its foundation is even older.
Integration is well behaved numerically, but poorly behaved algebraically.
Differentiation is poorly behaved numerically, but well behaved algebraically.
Therefore if one wants to do "automatic integration", one has to approach it in a radically different way to autodiff. And arguably, such a method does already exist; it's just very inefficient.
> Differentiation is poorly behaved numerically, but well behaved algebraically.
This is so pithy and true, am stealing that.
on why integration is a fundamentally harder problem than differentiation. New techniques would be required to programatically analytically solve integrals as well as current differentiation programs, but it would be exciting to see.
One way to think of this is to realize how easy it is to multiply numbers but how much more work it takes to divide numbers.
For something like automatic differentiation, you're essentially applying the chain rule for partial derivatives repeatedly. This is analytically pretty straightforward to do for most applications. All you need is an analytical derivative for the simple functions your more-complex function is comprised of (e.g. a neural network).
For integration, the analogue of the chain rule is integration by substitution [1]. The toolbox for solving integration problems is more limited than for differentiation. You run into issues where the answer cannot even be expressed using standard mathematical notation [2]. Sometimes you get lucky and the answer can be expressed via an alternating Taylor series so you can estimate the answer within some margin of error [3].
Stan is a piece of software that runs state-of-the-art MCMC methods to basically just compute fancy integrals. A Stan model will take an order of magnitude more time to run than a simple neural network via something like PyTorch on the same dataset. But they answer different questions.
[1] https://math.stackexchange.com/questions/1635949/is-there-a-...
[2] https://math.stackexchange.com/questions/1397132/why-cant-so...
[3] https://math.stackexchange.com/questions/145087/how-to-calcu...
When you input the variable values into a symbolic derivative you just get a value at the end. d/dx x^2 = 2x. If x = 0.5 then d/dx = 1. The same is true for symbolic integrals. For most practical applications, we don't really care about the full symbolic expression. We just want the answer, or at least a good approximation. This post uses a specific example of the difference between two Beta distributions. We want to get that 0.71. It is very hard to "automatically" make that happen.
Say we want the derivative of f(x) = x^3 - 2x^2 + 5. That becomes:
(x+e)^3 - 2(x+e)^2 + 5
= x^3 + 3x^2e - 2x^2 -4xe + 5
= (x^3 - 2x^2 + 5) + (3x^2 - 4x)e
The term in front of 'e' is "3x^2 - 4x", which is f'(x).
f'(x) = lim_{e->0} (f(x+e) - f(x))/e
If you're able to express f(x+e) on the form (f(x) + y e) then it follows that y is the derivative.
It also should be noted that auto-diff doesn't let you skip the rules for derivation and you're using the same calculation as you would to show that e.g. f'(x^n)=n*x^{n-1}.
But what about something like cos(x). Either you can lazily evaluate the power series or you know its sin(x).
What's the integral of exp( sqrt ( 1 + (tan^(3/2) X)2 ) ) ) ?
We only know a handful of forms that can be integrated in closed form and its down to our creativity to discover new forms that can be integrated (same deal with solving differential equations and the reasons are the same).
The forms that we know how to integrate can be done by a computer. CAS tools will do that for you. For example Mathematica.
Finding an "auto-integrate method" would probably involve finding a way of calculating the integral in a decomposable way, and that indeed would be amazing, but I don't really see that happening any time soon.
But then again, it's not clear what the OP meant by "automatic integration" anyway.
[0] https://newbooksnetwork.com/david-bressoud-calculus-reordere...
The problem is that it's unknown if a symbolic equality algorithm exists for elementary functions and what those functions should be.
It is known that the general case is undecidable: https://en.wikipedia.org/wiki/Richardson%27s_theorem
For integration we need to prove that such a object exists for each expression we want to integrate before we try integrating. Which as far as I'm aware is a problem equivalent to that solved by the Risch algorithm.
Can't we just integrate away and check if it's correct by differentiating it back?
EDIT: oh, I guess that's why we need to be able to check for function equality :)
[1] - https://en.wikipedia.org/wiki/Computable_analysis#Basic_resu...
The problem is that integration suffers from the curse of dimensionality much worse than differentiation.
Differentiation is a local activity. You approximate a derivative near any point. Integration is non-local, you must aggregate a function over a span of points in the domain space.
As the dimensions get higher, differentiation merely must apply local approximation to each successive dimension of a point in a grid. It only adds the computation cost.
Integration requires exponentially more data (for a given level of accuracy) because adding a new dimension enlarges the entire integration domain space that must be ranged over, and worse you generally need all that new “volume” of the enlarged space to be covered evenly to avoid bias.
So there are some fundamental asymmetries between differentiation and integration that play a role in why “automatic integration” isn’t as straightforwardly feasible as for differentiation.
Usually one cannot afford to visit all of the space. One visits high density regions preferentially.