In the software i'm contributing to, there's a procedure to take samples from a certain 2D distribution (if you're curious, it's for sampling a BSDF, used in rendering) and a procedure to return the corresponding pdf (over the hemisphere in terms of solid angle measure).

In order to test the implementation of the sample() and pdf() functions, i want to test if the binned sample counts agree with the bin probabilities $P_i$ i get from the pdf function (by integrating pdf over each bin). Basically i want to compare the two distributions in a significance test way and based on the result, let the unit test for the software pass or fail.

In the literature and similar software projects i mainly see people using Chi-squared tests, which can work to a degree, but produces problems with bins with zero counts or in general when the bin counts vary over several orders of magnitude.

Therefore i started thinking how to tackle this problem from a Bayesian standpoint and came up with the following:

  • I take n samples and bin them into N bins, interpreting them as counts drawn from a multinomial distribution.
  • I construct the posterior for the bin probabilities using a Jeffrey's prior, giving me a Dirichlet distribution: $p(\theta|\text{bincounts})=\text{Dirichlet}(\alpha_1+\frac{1}{2},\cdots,\alpha_N+\frac{1}{2}$)

Now my first idea was to compute a credible interval (volume in this case) from the Dirichlet distribution and check if the bin probabilities that i computed from pdf() are inside this region to decide on the test result.

But since i have on the order of a thousand bins, directly computing this region for a given significance level is not feasible.

So instead i had the idea to invert the idea and basically compute the significance level (or credibility) given our bin probabilities and Dirichlet distribution by doing the following:

  • calculate $\text{threshold}=p(\text{pdf}(P_1,\cdots,P_N)|\text{bincounts})$ by evaluating our posterior with the bin probabilities we got from pdf().
  • integrate the indicator function (showing which regions have lower posterior density than our threshold)$$ I(\theta) = \left\{ \begin{array}{lr} 1 :& 0 \lt p(\theta) < \text{threshold} \\ 0 :& \text{otherwise} \end{array} \right. $$ over the whole support of the Dirichlet posterior density.
  • to make this numerically tractable, estimate the integral by using importance sampling over the posterior density itself

This should then result in a number between 0 and 1 showing how much posterior 'volume' has lower density than our threshold giving us a kind of credibility that we can then use to decide if the test fails.

I implemented this in Mathematica and there my algorithm looks like this:

DirichletCredibility[observedcounts_, probs_, n_] := Module[
  {dist, mcsamples, pdf, threshold},
  dist = DirichletDistribution[observedcounts + 1/2];
  threshold = PDF[dist, Most[probs]];
  mcsamples = RandomVariate[dist, n];
  N@Mean[Composition[Boole[# <= threshold] &, PDF[dist]] /@ mcsamples]

I hope i explained well enough what i'm trying to achieve, if not i'm happy to clarify and give more details.

Now what happens when i use this is twofold: - if i test DirichletCredibility with bin counts drawn from a multinomial distribution and bin probabilities which don't match the multinomial probabilities the resulting credibility goes to zero with a high enough number of total samples, which is great! - if on the other hand i test DirichletCredibility with perfectly matching multinomial bin probabilities the result doesn't approach one, but instead tends to a uniform distribution (by repeating the credibility computation with different samples from the multinomial distribution).

My questions are: - Do you think my methodology can work in principle or does it have major flaws that i should address or consider? - If this is not feasible how would you tackle the problem in a Bayesian way to get a useful significance test for comparing sampled bin counts with reference probabilities?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.