Normalize the sky marginalization weights over the series, not the draws

The sky and time marginalized likelihood grew with the number of points
it was told to use: over a range of sixty four in the count it climbed by
four and a half in the log likelihood ratio, many times its own scatter,
so no count gave the same answer as any other.

The times are drawn from each detector's signal to noise series in
proportion to the likelihood along it, and what each carries is the
probability it was drawn with. That was being worked out over the drawn
times rather than over the series they came from, which makes it one in
however many were drawn, and leaves the answer growing by the log of the
count. The steps matched: 1.39, 1.69 and 1.41 against log 4 of 1.386.

Normalize over the series each detector drew from, as the time-only
marginalization already does. The same range of counts now moves the
answer by 0.36, within the scatter of the counts themselves, while the
error on it still falls with the count as it should.

Author: ahnitz

pycbc/inference/models/tools.py

@@ -526,6 +526,13 @@ def make_init():
 
             w = snr.squared_norm().numpy() / 2.0
             i = draw_sample(w, size=vsamples)
+            # the times were drawn from this detector's series in
+            # proportion to its likelihood, so normalizing over that
+            # series is what turns it into the probability they were
+            # drawn with. Normalizing over the drawn times instead would
+            # make it a probability of one in however many were drawn,
+            # and the answer would grow with the number of them.
+            w -= logsumexp(w)
 
             if sref is not None:
                 mcweight += w[i]
@@ -538,7 +545,6 @@ def make_init():
                 mcweight = w[i]
 
             idx.append(i)
-        mcweight -= logsumexp(mcweight)
 
         # check if delay is in dict, if not, throw out
         ti = []

test/test_marg_normalization.py

@@ -29,7 +29,13 @@
 from utils import simple_exit
 
 from pycbc.detector import Detector
-from pycbc.distributions import JointDistribution, SinAngle, Uniform
+from pycbc.distributions import (
+    CosAngle,
+    JointDistribution,
+    SinAngle,
+    Uniform,
+    UniformAngle,
+)
 from pycbc.inference import models
 from pycbc.noise import noise_from_psd
 from pycbc.psd import aLIGOZeroDetHighPower
@@ -81,6 +87,55 @@ def prior(self, halfwidth):
             Uniform(distance=(10, 200)),
             Uniform(tc=(TC - halfwidth, TC + halfwidth)))
 
+    def sky_marginalized(self, npoint, nseed=6):
+        """Mean and spread of the sky and time marginalized likelihood."""
+        variable = ['distance', 'inclination', 'tc', 'polarization',
+                    'ra', 'dec']
+        static = {k: v for k, v in self.static.items()
+                  if k not in ('ra', 'dec', 'polarization')}
+        dists = [SinAngle(inclination=None), Uniform(distance=(10, 200)),
+                 Uniform(tc=(TC - 0.1, TC + 0.1)),
+                 UniformAngle(polarization=None), UniformAngle(ra=None),
+                 CosAngle(dec=None)]
+        values = []
+        for s in range(nseed):
+            numpy.random.seed(700 + s)
+            model = models.MarginalizedTime(
+                list(variable), copy.deepcopy(self.data),
+                low_frequency_cutoff=self.flow, psds=self.psds,
+                static_params=static,
+                prior=JointDistribution(list(variable), *dists),
+                marginalize_phase=True,
+                marginalize_vector_params='tc,ra,dec,polarization',
+                marginalize_vector_samples=npoint, sample_rate=4096,
+                marginalize_sky_initial_samples=1e6)
+            model.update(**self.point)
+            values.append(model.loglr)
+        return numpy.mean(values), numpy.std(values) / nseed ** 0.5
+
+    def test_sky_marginalization_does_not_depend_on_the_point_count(self):
+        """The number of points must buy precision, not a bigger answer.
+
+        The times are drawn from each detector's signal to noise series,
+        so what they carry is their probability under that series. Taking
+        it over the drawn times rather than over the series made it one in
+        however many were drawn, and the answer grew by the log of that:
+        over a range of sixty four in the count it climbed by four and a
+        half, far outside its own scatter.
+        """
+        seen = [(n,) + self.sky_marginalized(n) for n in (128, 2048)]
+        (_, coarse, coarse_err), (_, fine, fine_err) = seen
+        self.assertLess(abs(coarse - fine),
+                        4 * (coarse_err ** 2 + fine_err ** 2) ** 0.5,
+                        "%s" % seen)
+
+    def test_more_sky_points_are_more_precise(self):
+        """What the count does buy is precision."""
+        _, coarse = self.sky_marginalized(128)
+        _, fine = self.sky_marginalized(2048)
+        self.assertLess(fine, coarse,
+                        "error on the mean went %.4f to %.4f" % (coarse, fine))
+
     def integrated(self, halfwidth, npoint=4001):
         """The marginal computed by summing the likelihood over the prior."""
         variable, prior = self.prior(halfwidth)