Bug fixes: correct observable prior and log evidence calculation

docs/src/fit-rv-astrom.md

@@ -57,20 +57,35 @@ Generate a simulated RV curve from the same orbit:
 using Random
 Random.seed!(1)
 
-epochs = 58849 .+ (20:25:90)
+epochs = 58849 .+ range(0,step=1.5, length=20)
 planet_sim_mass = 0.001 # solar masses here
 
 
 rvlike = MarginalizedStarAbsoluteRVLikelihood(
     Table(
         epoch=epochs,
-        rv=radvel.(orb_template, epochs, planet_sim_mass) ,
+        rv=radvel.(orb_template, epochs, planet_sim_mass) .+ 150,
         σ_rv=fill(5.0, size(epochs)),
     ),
-    jitter=:jitter1
+    jitter=:jitter1,
+    instrument_name="inst1"
 )
 
-Makie.scatter(rvlike.table.epoch[:], rvlike.table.rv[:])
+epochs = 58949 .+ range(0,step=1.5, length=20)
+
+rvlike2 = MarginalizedStarAbsoluteRVLikelihood(
+    Table(
+        epoch=epochs,
+        rv=radvel.(orb_template, epochs, planet_sim_mass) .- 150,
+        σ_rv=fill(5.0, size(epochs)),
+    ),
+    jitter=:jitter1,
+    instrument_name="inst2"
+)
+
+fap = Makie.scatter(rvlike.table.epoch[:], rvlike.table.rv[:])
+Makie.scatter!(rvlike2.table.epoch[:], rvlike2.table.rv[:])
+fap
 ```
 
 
@@ -92,14 +107,14 @@ end astrom
     M ~ truncated(Normal(1, 0.04),lower=0.1) # (Baines & Armstrong 2011).
     plx = 100.0
     jitter1 ~ truncated(Normal(0,10),lower=0)
-end rvlike b
+end rvlike rvlike2 b
 
 model = Octofitter.LogDensityModel(test)
 
 using Random
 rng = Xoshiro(0) # seed the random number generator for reproducible results
 
-results = octofit(rng, model, max_depth=9, iterations=5000)
+results = octofit(rng, model, max_depth=9, adaptation=300, iterations=400)
 ```
 
 Display results as a corner plot:
@@ -109,7 +124,7 @@ octocorner(model,results, small=true)
 
 Plot RV curve, phase folded curve, and binned residuals:
 ```@example 1
-Octofitter.rvpostplot(model, results,)
+Octofitter.rvpostplot(model, results)
 ```
 
 Display RV, PMA, astrometry, relative separation, position angle, and 3D projected views:

docs/src/rel-astrom-obs.md

@@ -21,6 +21,7 @@ astrom_like = PlanetRelAstromLikelihood(
     (epoch = 50720, ra = -469.1, dec = 114.3, σ_ra =   5.8, σ_dec =   5.8, cor=  0.1),
     (epoch = 50840, ra = -458.3, dec = 145.3, σ_ra =   4.2, σ_dec =   4.2, cor= -0.2),
 )
+obs_pri = ObsPriorAstromONeil2019(astrom_like)
 @planet b Visual{KepOrbit} begin
     e ~ Uniform(0.0, 0.5)
     i ~ Sine()
@@ -31,7 +32,7 @@ astrom_like = PlanetRelAstromLikelihood(
     a = ∛(system.M * b.P^2)
     θ ~ UniformCircular()
     tp = θ_at_epoch_to_tperi(system,b,50420)
-end astrom_like ObsPriorAstromONeil2019(astrom_like)
+end astrom_like obs_pri
 
 @system TutoriaPrime begin
     M ~ truncated(Normal(1.2, 0.1), lower=0.1)
@@ -62,22 +63,26 @@ Compare this with the previous fit using uniform priors:
     e ~ Uniform(0.0, 0.5) # hide
     i ~ Sine() # hide
     ω ~ UniformCircular() # hide
-    Ω ~ UniformCircular() # hide
+    # Ω ~ UniformCircular() # hide
+    # P = √(b.a^3/system.M) # hide
+    # θ ~ UniformCircular() # hide
+    # ω ~ Uniform(0,2pi) # hide
+    Ω ~ Uniform(0,2pi) # hide
     P = √(b.a^3/system.M) # hide
-    θ ~ UniformCircular() # hide
+    θ ~ Uniform(0,2pi) # hide
     tp = θ_at_epoch_to_tperi(system,b,50420) # hide
 end astrom_like # hide
 @system Tutoria begin # hide
     M ~ truncated(Normal(1.2, 0.1), lower=0.1) # hide
     plx ~ truncated(Normal(50.0, 0.02), lower=0.1) # hide
 end b # hide
-model = Octofitter.LogDensityModel(Tutoria) # hide
+model_with_uniform_priors = Octofitter.LogDensityModel(Tutoria) # hide
 Random.seed!(0) # hide
-results_unif_pri = octofit(model,iterations=5000,verbosity=0) # hide
-octoplot(model, results_unif_pri)
+results_unif_pri = octofit(model_with_uniform_priors,iterations=5000,verbosity=0) # hide
+octoplot(model_with_uniform_priors, results_unif_pri)
 ```
 
 We can compare the results in a corner plot:
 ```@example 1
 octocorner(model,results_unif_pri,results_obspri,small=true)
-```
+```

docs/src/rv-multi-planet.md

@@ -239,4 +239,31 @@ Octofitter.rvpostplot_animated(model_2p_v2, results_2p_v2)
 
 ```@raw html
 <video src="assets/rv-posterior.mp4" autoplay loop width=300 height=300>
+```
+
+
+
+```@example 1
+using Dynesty, HypercubeTransform
+chn_dyn_2p_v2 = Dynesty.dysample(model_2p_v2, DynamicNestedSampler())
+
+Octofitter.rvpostplot(model_2p_v2, chn_dyn_2p_v2, 1)
+
+```
+
+
+## Note about the evidence ratio
+The pigeons method returns the log evidence ratio. If the priors are properly normalized, this is equal to the log evidence.
+
+In other cases (e.g. if using `ObsPriorAstromONeil2019` or `UniformCircular`) you may need to calculate the log_Z0 term yourself. This can be done as follows:
+```@example 1
+prior_model = Octofitter.LogDensityModel(Octofitter.prior_only_model(pt_1p.system))
+_, pt_prior = octofit_pigeons(prior_model, n_rounds=10) # should be very quick!
+log_Z0 = pt_prior.shared.reports.summary.stepping_stone[end]
+```
+
+Subtract this from the stepping stone value to get the true evidence:
+```@example 1
+log_Z1_over_Z0 = pt_1p.shared.reports.summary.stepping_stone[end]
+log_Z1 = log_Z1_over_Z0 - log_Z0
 ```

src/cross-validation.jl

@@ -38,21 +38,21 @@ end
 Creates a copy of a system model that is stripped of observations. The result is 
 a model that only samples from the priors. This can be used eg. for tempering.
 """
-function prior_only_model(system::System)
+function prior_only_model(system::System;exclude_all=false)
 
 
     # Generate new observations for each planet in the system
     newplanets = map(1:length(system.planets)) do i
         planet = system.planets[i]
-        newplanet_obs = []
+        newplanet_obs = exclude_all ? AbstractLikelihood[] : filter(_isprior, planet.observations)
         newplanet = Planet{Octofitter.orbittype(planet)}(planet.priors, planet.derived, newplanet_obs..., name=planet.name)
         return newplanet
     end
 
-    newstar_obs = []
+    newsys_obs = exclude_all ? AbstractLikelihood[] : filter(_isprior, system.observations)
 
     # Generate new system with observations
-    newsystem = System(system.priors, system.derived, newstar_obs..., newplanets..., name=system.name)
+    newsystem = System(system.priors, system.derived, newsys_obs..., newplanets..., name=system.name)
 
     return newsystem
 end

src/likelihoods/observable.jl

@@ -116,11 +116,13 @@ function Octofitter.ln_like(like::ObsPriorAstromONeil2019{<:PlanetRelAstromLikel
         ) + 2*(-2+e^2+e*cos(E)) *sin(E))
     end
 
-    jac *= cbrt(P) / sqrt(1-eccentricity(orbit)^2);
+    sqrt_eccen = sqrt(1-eccentricity(orbit)^2)
+    jac *= cbrt(P) / sqrt_eccen 
 
-    ln_prior += - 2log(jac)
+    ln_prior += 2log(jac)
 
     return ln_prior
 end
 
+
 # TODO: Add a RadialVelocity correction version in OctofitterRadialVelocity