Say you’ve got a positive dataset and you want to calculate the variance. However, the numbers in your dataset are huge, so huge you need to represent them in the log-domain. How do you compute the log-variance without things blowing up?
I ran into the problem today. To my surprise, I couldn’t find a standard solution.
The bad solution
Suppose that your data is 
This of course is a terrible idea: When 
The mediocre solution
The first idea I had for this problem was relatively elegant. We can of course represent the variance as
Instead of calculating 
To do this, we can introduce a “log domain mean” operator, a close relative of the good-old scipy.special.logsumexp
def log_domain_mean(logx):
"np.log(np.mean(np.exp(x))) but more stable"
n = len(logx)
damax = np.max(logx)
return np.log(np.sum(np.exp(logx-damax))) \
+ damax-np.log(n)
Next, introduce a “log-sub-add” operator. (A variant of np.logaddexp)
def logsubadd(a,b):
"np.log(np.exp(a)-np.exp(b)) but more stable"
return a + np.log(1-np.exp(b-a))
Then, we can compute the log-variance as
def log_domain_var(logx):
a = log_domain_mean(2*logx)
b = log_domain_mean(logx)*2
c = logsubadd(a,b)
return c
Here a is 
b is 
Nice, right? Well, it’s much better then the first solution. But it isn’t that good. The problem is that when the variance is small, a and b are close. When they are both close and large, logsubadd runs into numerical problems. It’s not clear that there is a way to fix this problem with logsubadd.
To solve this, abandon elegance!
The good solution
For the good solution, the math is a series of not-too-intuitive transformations. (I put them at the end.) These start with
and end with
Why this form? Well, we’ve reduced to things we can do relatively stably: Compute the log-mean, and do a (small variant of) log-sum-exp.
def log_domain_var(logx):
"""like np.log(np.var(np.exp(logx)))
except more stable"""
n = len(logx)
log_xmean = log_domain_mean(logx)
return np.log(np.sum( np.expm1(logx-log_xmean)**2))\
+ 2*log_xmean - np.log(n)
This uses the log_domain_mean implementation from above, and also np.expm1 to compute 
a is close to zero.
Why is this stable? Is it really stable? Well, umm, I’m not sure. I derived transformations that “looked stable” to me, but there’s no proof that this is best. I’d be surprised if a better solution wasn’t possible. (I’d also be surprised if there isn’t a paper from 25+ years ago that describes that better solution.)
In any case, I’ve experimentally found that this function will (while working in single precision) happily compute the variance even when logx is in the range of 
As always, failure cases are probably out there. Numerical instability always wins when it can be bothered to make an effort.
. Assume we observe some fixed 
and want to sample from 
Assume 
, return 
. Thus
So this gives exact samples.
and (ii) to be able to simulate 
. You don’t even need to be able to evaluate either of these!
will typically be small unless you tend to get the same data regardless of the value of the latent variable. (In which case is the problem even interesting?)
you need a proposal distribution 
you can sample from and you need to know a constant such that 
The above algorithm is just using
is equivalent to 
which is always true.
. In the ABC setting we can’t calculate either of these, but we exploit that we can calculate the ratio 
for some small 
and accept if 
rather than if 
as before.
in any one iteration is![\displaystyle \begin{array}{rcl} \text{probability of returning }z^{*} & = & p(z^{*})\sum_{x}p(x\vert z^{*})I[s(x)=s(x^{*})].\\ & = & p(z^{*})p(s(x)\vert z^{*})\\ & = & p(s(x))p(z^{*}\vert s(x)). \end{array}](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle++%5Cbegin%7Barray%7D%7Brcl%7D++%5Ctext%7Bprobability+of+returning+%7Dz%5E%7B%2A%7D+%26+%3D+%26+p%28z%5E%7B%2A%7D%29%5Csum_%7Bx%7Dp%28x%5Cvert+z%5E%7B%2A%7D%29I%5Bs%28x%29%3Ds%28x%5E%7B%2A%7D%29%5D.%5C%5C+%26+%3D+%26+p%28z%5E%7B%2A%7D%29p%28s%28x%29%5Cvert+z%5E%7B%2A%7D%29%5C%5C+%26+%3D+%26+p%28s%28x%29%29p%28z%5E%7B%2A%7D%5Cvert+s%28x%29%29.+%5Cend%7Barray%7D+&bg=ffffff&fg=000000&s=0 "\displaystyle \begin{array}{rcl} \text{probability of returning }z^{*} & = & p(z^{*})\sum_{x}p(x\vert z^{*})I[s(x)=s(x^{*})].\\ & = & p(z^{*})p(s(x)\vert z^{*})\\ & = & p(s(x))p(z^{*}\vert s(x)). \end{array}")
![{p(s'\vert z)=\sum_{x}p(x\vert z)I[s'=s(x)]}](https://s0.wp.com/latex.php?latex=%7Bp%28s%27%5Cvert+z%29%3D%5Csum_%7Bx%7Dp%28x%5Cvert+z%29I%5Bs%27%3Ds%28x%29%5D%7D&bg=ffffff&fg=000000&s=0 "{p(s'\vert z)=\sum_{x}p(x\vert z)I[s'=s(x)]}")
and ![{p(s')=\sum_{x}p(x)I[s(x)=s'].}](https://s0.wp.com/latex.php?latex=%7Bp%28s%27%29%3D%5Csum_%7Bx%7Dp%28x%29I%5Bs%28x%29%3Ds%27%5D.%7D&bg=ffffff&fg=000000&s=0 "{p(s')=\sum_{x}p(x)I[s(x)=s'].}")
This could be much higher if there are many different datasets that yield the same summary statistics.
is in the 
. This is very nice.
. If our proposal distribution was 
then we’d do
from some proposal distribution
.
then 
.
instead. We will suggest the following algorithm instead, with changes shown in blue.
.
![{{\displaystyle I[\tilde{x}^{*}=x^{*}]}}](https://s0.wp.com/latex.php?latex=%7B%7B%5Cdisplaystyle+I%5B%5Ctilde%7Bx%7D%5E%7B%2A%7D%3Dx%5E%7B%2A%7D%5D%7D%7D&bg=ffffff&fg=0000ff&s=0 "{{\displaystyle I[\tilde{x}^{*}=x^{*}]}}")
.
, we generate a synthetic dataset. We can accept only if the synthetic dataset is the same as the observed one.
may be a terrible proposal distribution for the posterior 
only yields nearby points, than once the typical set has been reached, the probability of propsing a “good” ![\displaystyle p(z,x^{*},x)=p(z,x^{*})I[x^{*}=x].](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle++p%28z%2Cx%5E%7B%2A%7D%2Cx%29%3Dp%28z%2Cx%5E%7B%2A%7D%29I%5Bx%5E%7B%2A%7D%3Dx%5D.+&bg=ffffff&fg=000000&s=0 "\displaystyle p(z,x^{*},x)=p(z,x^{*})I[x^{*}=x].")
using Metropolis-Hastings. We choose the proposal distribution
![\displaystyle \begin{array}{rcl} \alpha & = & \min\left(1,\frac{p(\tilde{z},\tilde{x}^{*},x)}{p(z,x^{*},x)}\frac{q(z,x^{*}\vert\tilde{z},\tilde{x}^{*})}{q(\tilde{z},\tilde{x}^{*}\vert z,x^{*})}\right)\\ & = & \min\left(1,\frac{p(\tilde{z},\tilde{x}^{*})I[\tilde{x}^{*}=x]}{p(z,x^{*})I[x^{*}=x]}\frac{q(z\vert\tilde{z})p(x^{*}\vert z)}{q(\tilde{z}\vert z)p(\tilde{x}^{*}\vert\tilde{z})}\right) \end{array} .](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle++%5Cbegin%7Barray%7D%7Brcl%7D++%5Calpha+%26+%3D+%26+%5Cmin%5Cleft%281%2C%5Cfrac%7Bp%28%5Ctilde%7Bz%7D%2C%5Ctilde%7Bx%7D%5E%7B%2A%7D%2Cx%29%7D%7Bp%28z%2Cx%5E%7B%2A%7D%2Cx%29%7D%5Cfrac%7Bq%28z%2Cx%5E%7B%2A%7D%5Cvert%5Ctilde%7Bz%7D%2C%5Ctilde%7Bx%7D%5E%7B%2A%7D%29%7D%7Bq%28%5Ctilde%7Bz%7D%2C%5Ctilde%7Bx%7D%5E%7B%2A%7D%5Cvert+z%2Cx%5E%7B%2A%7D%29%7D%5Cright%29%5C%5C+%26+%3D+%26+%5Cmin%5Cleft%281%2C%5Cfrac%7Bp%28%5Ctilde%7Bz%7D%2C%5Ctilde%7Bx%7D%5E%7B%2A%7D%29I%5B%5Ctilde%7Bx%7D%5E%7B%2A%7D%3Dx%5D%7D%7Bp%28z%2Cx%5E%7B%2A%7D%29I%5Bx%5E%7B%2A%7D%3Dx%5D%7D%5Cfrac%7Bq%28z%5Cvert%5Ctilde%7Bz%7D%29p%28x%5E%7B%2A%7D%5Cvert+z%29%7D%7Bq%28%5Ctilde%7Bz%7D%5Cvert+z%29p%28%5Ctilde%7Bx%7D%5E%7B%2A%7D%5Cvert%5Ctilde%7Bz%7D%29%7D%5Cright%29+%5Cend%7Barray%7D+.&bg=ffffff&fg=000000&s=0 "\displaystyle \begin{array}{rcl} \alpha & = & \min\left(1,\frac{p(\tilde{z},\tilde{x}^{*},x)}{p(z,x^{*},x)}\frac{q(z,x^{*}\vert\tilde{z},\tilde{x}^{*})}{q(\tilde{z},\tilde{x}^{*}\vert z,x^{*})}\right)\\ & = & \min\left(1,\frac{p(\tilde{z},\tilde{x}^{*})I[\tilde{x}^{*}=x]}{p(z,x^{*})I[x^{*}=x]}\frac{q(z\vert\tilde{z})p(x^{*}\vert z)}{q(\tilde{z}\vert z)p(\tilde{x}^{*}\vert\tilde{z})}\right) \end{array} .")
was accepted, it must be true that ![\displaystyle \begin{array}{rcl} \alpha & = & \min\left(1,\frac{p(\tilde{z})}{p(z)}\frac{q(z\vert\tilde{z})}{q(\tilde{z}\vert z)}\right)I[\tilde{x}^{*}=x]. \end{array}](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle++%5Cbegin%7Barray%7D%7Brcl%7D++%5Calpha+%26+%3D+%26+%5Cmin%5Cleft%281%2C%5Cfrac%7Bp%28%5Ctilde%7Bz%7D%29%7D%7Bp%28z%29%7D%5Cfrac%7Bq%28z%5Cvert%5Ctilde%7Bz%7D%29%7D%7Bq%28%5Ctilde%7Bz%7D%5Cvert+z%29%7D%5Cright%29I%5B%5Ctilde%7Bx%7D%5E%7B%2A%7D%3Dx%5D.+%5Cend%7Barray%7D+&bg=ffffff&fg=000000&s=0 "\displaystyle \begin{array}{rcl} \alpha & = & \min\left(1,\frac{p(\tilde{z})}{p(z)}\frac{q(z\vert\tilde{z})}{q(\tilde{z}\vert z)}\right)I[\tilde{x}^{*}=x]. \end{array}")
![{I[\tilde{x}^{*}=x]}](https://s0.wp.com/latex.php?latex=%7BI%5B%5Ctilde%7Bx%7D%5E%7B%2A%7D%3Dx%5D%7D&bg=ffffff&fg=000000&s=0 "{I[\tilde{x}^{*}=x]}")
into ![{I[s(\tilde{x}^{*})=s(x)].}](https://s0.wp.com/latex.php?latex=%7BI%5Bs%28%5Ctilde%7Bx%7D%5E%7B%2A%7D%29%3Ds%28x%29%5D.%7D&bg=ffffff&fg=000000&s=0 "{I[s(\tilde{x}^{*})=s(x)].}")
You can also add a slop factor if you want.
can be something interesting like a Multivariate Gaussian. The acceptance probability then instead becomes
defined in 2.10, which itself involves other assumptions tracing backwards in the paper.
)
this becomes
and 
with 
and you need to choose 
so as to ensure that 
.
, then 
, but this is awkward with multiple lines.
. You don’t know 
times.
heads and 
tails.
: the probability of seeing exactly 
.
and 
.
, with 
flips, you tend to see around 
heads. Here’s an example with 
. Here, the most common value is 
, close to 
.
.
is defined comparing 
to 
, and so lose 
is a decision rule– a precise recipe of what decision to make in each possible situation.
such decision rules, corresponding to 
, (always predict heads), 
(predict heads if you see at least one heads), up to 
(always predict tails). These are shown here:
and 
, you don’t lose anything by restricting to this kind of decision rule. However, there are losses for which these decision rules are not enough. (Imagine you lose more when your guess is correct.)
) with loss 
in the bottom plot is the total area of the green bars in the middle. You can do the same sweep for 
, which you can is pictured here:
?
over 
:
are shown in light grey. Naturally enough, for each value of 
. For other losses, these simple ones might not be enough, and there are a lot of decision rules. (Even for this toy problem there are 
, since you can output heads or tails for each of 
, 
, …, 
.)
:
defined above, but we didn’t use this assumption in deriving the rules.)
)) then the optimal Bayesian decision rule is determined.
. However, on any particular day, some other decision rule could certainly be better.
, while 
. That is, the Bayesian method predicts heads when it observes 2 or more, while the minimax rule predicts heads if it observes even one. This makes sense intuitively: The minimax decision rule proceeds as if the “worst” w (a small number) is fixed, whereas the Bayesian decision rule less pessimistically averages over all w.
is often computationally hard, since you need to integrate over all possible observations. In general, finding exact Bayes or minimax optimal decision rules is a huge computational challenge, and at least some degree of approximation is required. That’s the game, that’s why statistics is hard. Still, even for complex situations the rules are the same– you win by finding a decision rule with low risk.
for some simple class of distributions 
parameterized by 
is a mixture of two Gaussians (in red), and 
is a single Gaussian (in blue), the VI optimization in one dimension would arrive at the solution below.
. The random walk is carefully constructed so that if you run it a long time, the probability it will end up in any given state is proportional to 
and a target distribution 
(so we can do a random walk) so that
.
. That’s difficult since 
is defined by marginalizing 
and then you have the bound
for any 
in the [0,1] interval.
depends on the choice of 
where 
is a normalizing constant. Here, you can think of 
as something akin to a “base measure”. 
is restricted to beconstant over 
![s(w) = \log r(w) - \log r_z + E_{q(Z|w)} [\beta^{-1} \log p(Z) + (1-\beta^{-1}) \log q(Z|w)]](https://s0.wp.com/latex.php?latex=s%28w%29+%3D+%5Clog+r%28w%29+-+%5Clog+r_z+%2B+E_%7Bq%28Z%7Cw%29%7D+%5B%5Cbeta%5E%7B-1%7D+%5Clog+p%28Z%29+%2B+%281-%5Cbeta%5E%7B-1%7D%29+%5Clog+q%28Z%7Cw%29%5D&bg=ffffff&fg=000000&s=0 "s(w) = \log r(w) - \log r_z + E_{q(Z|w)} [\beta^{-1} \log p(Z) + (1-\beta^{-1}) \log q(Z|w)]")
up to a constant, i.e.
.
is noise from a standard Gaussian and 
is a step-size.
derived above. Then, after a re-scaling, this becomes the new hybrid algorithm
is the entropy of 
. This clearly becomes the previous VI algorithm in the limit of 
. It also essentially becomes Langevin with 
. That’s because the distribution 
(not yet defined!) will prefer 
is highly concentrated. Thus, only the mean parameters of 
and 
. Second, Tthe gradient of the objective involves a KL-divergence. Exactly computing this is intractable, but can be approximated with standard tricks from stochastic VI, namely data subsampling and the “reparameterization trick”. Third, 
where 
is a universal constant chosen for each 
is pretty simplistic. (Recall, 
. Using the fact that 
you can show that it’s possible to find derivatives of 
online, and thus tighten the bound while the algorithm is running. (I didn’t want to do this in this paper since neither VI nor MCMC do this, and it complicates the interpretation of the experiments.)
. So, if you want to estimate a derivative, an easy way to do so would be to just pick some small 
, using a range of different 
, we get:
?
?
) and Betty (
) live in a parallel universe with two kinds of blocks, green (
) and yellow (
). Alfred likes green blocks, and Betty prefers the yellow blocks. One day, a Wizard decides to give one of the monkeys the magical power to send one block over to our universe each day.
. That is, the Wizard chooses Alfred with probability 
and Betty with probability 
. Both Alfred and Betty send their favorite colored block with probability 
We can mechanically apply 
that we haven’t explicitly calculated) we can normalize to obtain that
?
or rainy 
. Both the monkeys prefer clear weather. When the weather is clear, they send their preferred block with probability 
That is, we have
. Now, what is the probability that Alfred is the magical monkey?
is constant between the two, we can normalize it out. Thus, we get that
)
)
)
indexes the different days of observation.
. So, we simply sum over the different models
where 
is a noisy unbiased estimate of the gradient at 
and 
is some sequence of step-sizes that obeys the Robbins-Monro conditions [1]. However, in practice, the difference between, say, 
, and 
can be enormous, and finding these constants is a bit of a dark art.
and 
.)
is not linear in 
complexity to 
. (Calculate 
such that 
.) This is particularly true when the lower-complexity algorithm has higher constant factors. As an example, I’ve always found the 
ball
