Log Gamma approximation

By the standard I set out yesterday, C#'s math library is less than decent. I was a bit surprised that even the stats library didn't have a function for the log of gamma. Fortunately, there's no secret to the approximation. Most implementations use a variant of the Stirling Approximation, which can be super accurate if you want it to be. I don't really need much more than a few decimal places, so I coded up the Nemes version which only uses a few terms and, according to my tests, gives 6-digits when z is greater than 1.5. The CISS confidence intervals are pretty meaningless without at least 1 block sampled, so the lower bound on the range I care about is 2.0. You can read the wiki page (or the source article, if you know German) for the mathematical details. I'll simply post the code for my C# brethren:

        // log gamma function using simplified Stirling's formula
        // Gergő Nemes: New asymptotic expansion for the Gamma function
        // Archiv der Mathematik, August 2010, Volume 95, Issue 2, pp 161-169
        public static double LogGamma(double z)
        {
            return (Math.Log(2.0 * Math.PI) - Math.Log(z)) / 2.0 +
                z * (Math.Log(z + 1.0 / (12.0 * z - 1.0 / (10.0 * z))) - 1.0);
        }

Closed form estimates

There's a reason Bayesian's have been using the Beta distribution for Binomial and Bernoulli priors for so long: it does tie up into a closed form rather nicely. It may turn out to be not the best choice, but if it is workable, I don't have to do any numerical estimation to get the variance on CISS results.

Recall from yesterday that, if X is a random variable representing the query sum from a block:

E(X|θ) = μ(θ) = θnb_k/2

E(X) = μ = ∫μ(θ)P(θ)dθ    where P(θ) is the posterior on θ

Factoring out the constant from μ(θ), this simply becomes:

μ = (nb_k/2)∫θP(θ)dθ = E(X|block contains data)E(θ)

Well, duh. As long as E(X|block contains data) is independent of θ, that's always going to be true. Unfortunately, that's a flimsy assumption. If the query attributes are highly correlated within a block, θ drops and E(X|block contains data) rises. But, we're deferring that for the moment. We will have to revisit it when we get to dynamic partitioning because then we'll be intentionally increasing the negative correlation between the two. (We may have to look at it sooner; I don't trust any assumptions until I run real data past them).

At any rate, as long as we're using a distribution with a closed-form mean (Beta certainly qualifies), we can move on to the variance:

Var(X|θ) = σ(θ)² = E(|X - μ|²) = (1-θ)μ² + θ[(nb_k)²/3 - nb_k μ + μ²]

Var(X) = σ² = ∫σ(θ)²P(θ)dθ

Before we go plowing ahead with the integration, let's simplify our symbols a bit by letting p = μ² and q = (nb_k)²/3 - nb_k μ + μ². This helps us see that these terms are just constant coefficients and can be moved in and out of the integral giving:

σ(θ)² = (1-θ)p + θq

σ² = ∫σ(θ)²P(θ) dθ = ∫[(1-θ)p + θq]P(θ) dθ

Substituting in the density for θ ~ Beta(a,b) gives

σ² = ∫[(1-θ)p + θq]θ^a-1(1 - θ)^b-1 dθ

   = ∫(1-θ)^bpθ^a-1 + θ^aq(1 - θ)^b-1 dθ

   = p∫θ^a-1(1-θ)^b dθ + q∫θ^a(1 - θ)^b-1 dθ

   = pΒ(a,b+1) + qB(a+1,b)   where B(a,b) is the Beta function.

Now, before you accuse me of cheating and simply burying the integrals inside a predefined function, recall that:

B(a,b) = Γ(a)Γ(b) / Γ(a+b)

and that the Γ function is simply an extension of factorial to the real number line:

Γ(i) = (i-1)!

So, if a and b are integer:

B(a,b) = (a-1)!(b-1)! / (a+b-1)!

and

σ² = p[(a-1)!b! / ((a+b)!] + q[a!(b-1)! / (a+b)!]

OK, factorial isn't technically closed form, but since we're dealing with block counts rather than row counts, we're not talking about massive looping to get answers. In fact, while it's no problem to set the prior so I only get integer parameters a and b on the posterior, any decent math software library has a closed form approximation of the Gamma function, so I could just call that.

E Pluribus Unum

Don't panic, I'm not switching this over to be a political blog. It did occur to me though, that the nature of my data correlation allows for a nice trick which may simplify matters considerably. To review:

• We want a weak prior on the proportion of values returned.
• We want a somewhat stronger prior on the magnitude of those values.
• Unless we adjust for correlation, the posterior on the proportion will be too tight because we'll give too much credence to the first data points.
• Correlation is strong within a block.
• Correlation is weak (nonexistent?) between blocks.
• All we really care about is the sum of values out of a block (actually, we only care about the sum for the whole stratum, but we can only read one block at a time so we have to put that sum somewhere).

So...

Why not simply ditch the individual rows altogether and just work with block sums? Those are very nicely distributed random variables:

P(X = x) = {1-θ for x = 0;  θU(0, nb_k) for x ≠ 0}
where U(0, nb_k) is uniform between zero and the stratum upper bound times the block size. (In the case of the negative strata, it will be U(nb_k, 0), which flips the sign on the expected value but doesn't change the variance.) Thus,

E(X|θ) = μ(θ) = θnb_k/2
E(X) = μ = ∫μ(θ)P(θ)dθ    where P(θ) is the posterior on θ

and

Var(X|θ) = σ(θ)² = E(|X - μ|²) = (1-θ)μ² + θ[(nb_k)²/3 - nb_k μ + μ²]
Var(X) = σ² = ∫σ(θ)²P(θ)dθ

If the second term in Var(X|θ) looks weird, that's because we're looking at the variability around the posterior mean, E(X), rather than the mean of the uniform random variable. If you substitute nb_k/2 for μ, you get the familiar uniform variance of (nb_k)²/12 multiplied by θ for that term, but that puts us right back where we were last Friday.

Even with something as computationally convenient as the Beta distribution for θ, those integrals are going to be messy. In the early going, I can get by with a rough approximation. As the confidence interval gets closer to the goal, I'll need to be more accurate. Fortunately, the distribution of θ will tighten up as well, so I may not have to integrate over the entire domain.

Also note that in the first integral, it's just a constant times the variable times the density. The second isn't much worse: μ is fixed, so we're just placing the density on top of the affine space of two constants; it's really not that bad.

In exchange for the additional computation, this allows me to intentionally increase the correlation within a block, which I certainly do intend to do when we get to dynamic partitioning. Being able to know that a block can be skipped because the partition contains no data in the query filter will more than pay back grinding out a numerical integral. As long as I'm just concerned with the sum, I can still compute a confidence interval without a bunch of correlation adjustments.

Next step is to get these ideas into CISS and test them out. Code modifications should be pretty straightforward. Really hoping to have a stable version this week so I can get to writing.

Prior resurrected

It was nice to take a couple days off for Easter. Classes are off all week for Spring Break, but I'm back on it.

I know why my confidence intervals are too tight. This should have been obvious from the get go, but sometimes you need to take a step back to see the obvious. The confidence intervals are based on the distribution of the values given θ and λ. I just took the point estimates for those two. However, they aren't constants, they are random variables themselves. So, I really need to compute the distribution given θ and λ and then integrate that over the joint distribution of the inputs. That's not terribly hard to do, but it does bring us back the problem that I dismissed too quickly: I don't even have good marginal distributions for θ and λ. I don't have a clue what the joint distribution would be.

So, the first thing I'm going to do is take lambda out of the equation. Rather than stratify on absolute value, I'm going to stratify on value. Thus, within a stratum, values will be considered uniform between the stratum bounds and the big question is how many get picked up by the query (θ).

Because this doubles the number of strata, I'm going to back off on pinning the magnitude to 2^k and instead allow for configurable bounds on each strata. I've got some ideas on the optimal way to do that, but for now, I'll just set them manually.

That leaves the prior on θ. If I'm really going to use this as a distribution for purposes of confidence intervals, I need to play by the rules and compute a real posterior based on likelihood. So, the fast-converging point estimate isn't going to cut it. I think I can still get around the problem by letting observations in adjacent strata inform the posterior. So, we start with a fairly neutral prior that reflects the variability due to the fact that we don't know the value of θ and the posterior gets updated anytime we sample a block from strata k-1, k, or k+1. This will speed the convergence, while still preserving a lot of variability, at least in the early going.

Although mathematical convenience makes it tempting to use the Beta distribution as the prior, I think it converges much too quickly, especially since we have correlation within blocks. I'm going to have to work something else out. This seems like a good application of Gibbs sampling. I can try a bunch of different distributions and see which ones behave right.

Overshot

Looks like I tightened up the confidence interval just a bit too much. Here's the convergence plot for the same query posted earlier in the week. This time the uncertainty bounds are "real" confidence intervals (except that they're wrong). Ideally, the true value would be between the blue lines 95% of the time.

The convergence path is just a bit slower due to the fact that the algorithm thinks it's doing better than it is so it's not picking strata as well. More importantly, the stopping rule is the width of the confidence interval, so the algorithm terminates before the desired confidence has been reached.  I had it set to $1000 to ensure it ran out nearly to the end (note, the scale on the left is in billions). For an all-expense query like this, the actuaries would probably want things to the nearest $50-million or so, which would end it at block 282, but it was still off by $225 million at that point.

Not sure which of my assumptions is being violated in a bad way, but I'm going to leave it at this for the week. Definitely good enough progress to turn my attention to far greater matters for a few days. Happy Easter; I won't be posting again until Monday.

Prior on hit rate

Despite a couple drawbacks, the Beta distribution is the obvious choice for my prior on the proportion of rows returned by a query. For those who haven't done much with Bayesian stats, the Beta is a favorite for proportions because it's a "conjugate prior", meaning that the posterior distribution is of the same form. This makes iterative application pretty easy.

It's more than just a mathematical convenience; it really does do a nice job of distributing uncertainty. Don't have a clue? Use Beta(1,1) and your posterior is simply the proportion in your data with confidence intervals that match "traditional" frequentist models. Sort of sure, but willing to have your mind changed? Use Beta(a, b) where a/(a+b) is your believed proportion and a+b indicates how sure you are of that. If your data size is greater than a+b, the posterior will be biased towards the data. If your data size is smaller, the prior will show through more strongly. The graph below illustrates:

Here, the prior indicates a belief that we should get 16 hits for every 36 misses. We sample the data and get 57 hits and 43 misses. The posterior now predicts (16+57) hits for every (36+43) misses. The confidence interval has been tightened appropriately. We could now apply more data using the posterior as the new prior and get the same final result as if we had applied all the data at once. Convenient, intuitive, and flexible.

But, the real world is often messier than that. Because these aren't truly random samples, but samples of blocks of correlated rows, I need the posterior to reflect a pessimistic view of the uncertainty. Specifically, I want it to assume that there's a lot of relevant data out there that we simply haven't gotten to yet (until the data we have gotten to overwhelmingly suggests otherwise). That would suggest using Beta(a, b) where a is very large and b is very small, even though I don't really believe that. Specifically, I'd like to assume that EVERY row is included. The problem with setting a prior to 0% or 100% is that it's no longer a distribution at all, simply a statement of fact. As such, any data that doesn't match that fact has a likelihood of zero and the posterior degenerates back to the prior. (There's a Bayesian "proof" of God's existence that counts on people missing this detail).

There is an out. I don't actually need a prior on θ for the first iteration. I only need the point estimate to compute the variance of the observations. I can just set it to 1 and compute the confidence interval for the stratum (which will obviously be quite wide). Then, after a block has been sampled, I can set the posterior to Beta(N+h, m) where h and m are the number of hits and misses from the first block and N is the total number of rows.

The rub is if the first block has no misses. Then I'm back to a degenerate prior. That's not really terrible. If I just continue to add hits and misses to the distribution, the function will work; it just doesn't make a lot of intuitive sense.

The real problem is that this isn't going to collapse fast enough. Remember, I don't really believe this prior. I'm just setting it that way to keep the algorithm sampling. After sampling all the data, I'll have a distribution on θ of Beta(N+H, M) where H and M are the total hits and misses for the stratum. That's clearly nonsense because everything has been sampled, so there's no uncertainty left. It should be just θ = H/(M + H). So, I'm going to use a little sleight of hand. Instead of treating the sampling as cumulative, I'm going to treat it as replacing. That is, a sampled block will now replace an unsampled block in the prior, not the posterior. Under that scheme, the prior becomes Beta(N-(h+m), 0) and the posterior is Beta(N-m, m). Repeating this over all subsequent samplings will result in a final state of Beta(H, M), which is exactly what we're looking for.

Mathematically, this is a bit bogus and I'd not feel good about it if I was giving a confidence interval on θ. However, I'm not. I'm just using a point estimate of θ to compute a confidence interval on the stratum sum. So, at each step, I'll just use θ = (N-M)/N where M is the cumulative misses so far. This gives a well-defined point estimate at each step and converges to the real answer if we end up sampling the entire stratum. Given that, I don't need to stress over the actual distribution.

Variance by strata

It's just symbol manipulation, but there's enough of it to be worth writing down. So, I'll write it here. This is the derivation of the variance for an individual observation from stratum k.

P(X = 2^k) = θλ
P(X = -2^k) = θ(1-λ)
P(X = 0) = 1-θ

stratum mean μ = ΣxP(x) = 2^kθλ - 2^kθ(1-λ) = 2^kθ(2λ - 1)

stratum variance σ² = Σ(μ-x)²P(x)

   = θλ(μ - 2^k)² + θ(1 - λ)(μ + 2^k)² + (1 - θ)μ²

   = θλμ² - 2θλμ2^k + θλ2^2k + θ(1-λ)μ² + 2θ(1-λ)μ2^k + θ(1-λ)2^2k + (1-θ)μ²

   = μ² + 2θ(1-2λ)μ2^k + θ2^2k

Substituting in the formula for the mean to get the variance just in terms of θ and λ is not enlightening and doesn't help the computational stability, so I'll leave it as an exercise to any reader that just really likes moving letters around.