Merge branch 'fix_marginalised_likelihood' into 'master'

Fix marginalised likelihood

See merge request Monash/tupak!45

examples/injection_examples/marginalized_likelihood.py

@@ -20,10 +20,9 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
                             waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
 
 # Create the waveform_generator using a LAL BinaryBlackHole source function
-waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
-                                                                sampling_frequency=sampling_frequency,
-                                                                frequency_domain_source_model=tupak.source.lal_binary_black_hole,
-                                                                parameters=injection_parameters)
+waveform_generator = tupak.waveform_generator.WaveformGenerator(
+    time_duration=time_duration, sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=tupak.source.lal_binary_black_hole, parameters=injection_parameters)
 hf_signal = waveform_generator.frequency_domain_strain()
 
 # Set up interferometers.
@@ -39,8 +38,9 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota
 
 # Initialise GravitationalWaveTransient
 # Note that we now need to pass the: priors and flags for each thing that's being marginalised.
-likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator, prior=priors,
-                                         distance_marginalization=True, phase_marginalization=True)
+likelihood = tupak.likelihood.GravitationalWaveTransient(
+    interferometers=IFOs, waveform_generator=waveform_generator, prior=priors,
+    distance_marginalization=True, phase_marginalization=True)
 
 # Run sampler
 result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',

tupak/likelihood.py

@@ -14,10 +14,10 @@ import logging
 class Likelihood(object):
     """ Empty likelihood class to be subclassed by other likelihoods """
 
-    def log_likelihood():
+    def log_likelihood(self):
         return np.nan
 
-    def noise_log_likelihood():
+    def noise_log_likelihood(self):
         return np.nan
 
     def log_likelihood_ratio(self):
@@ -119,11 +119,20 @@ class GravitationalWaveTransient(Likelihood):
             matched_filter_snr_squared_array = matched_filter_snr_squared * \
                 self.waveform_generator.parameters['luminosity_distance'] / self.distance_array
 
+            if self.phase_marginalization:
+                matched_filter_snr_squared_array = self.bessel_function_interped(abs(matched_filter_snr_squared_array))
+            else:
+                matched_filter_snr_squared_array = np.real(matched_filter_snr_squared_array)
+
             log_l = logsumexp(matched_filter_snr_squared_array - optimal_snr_squared_array / 2,
                               b=self.distance_prior_array * self.delta_distance)
-        else:
+        elif self.phase_marginalization:
+            matched_filter_snr_squared = self.bessel_function_interped(abs(matched_filter_snr_squared))
             log_l = matched_filter_snr_squared - optimal_snr_squared / 2
 
+        else:
+            log_l = matched_filter_snr_squared.real - optimal_snr_squared / 2
+
         return log_l.real
 
     def log_likelihood(self):