Proposals

[zenna]

It's tricky to get the interface to proposals right.

Let's consider the requirements for different inference methods

• MH: sample proposal, and compute proposal density, compute p(x) and ratio
• InvolutiveMCMC - sample auxiliary proposal and define involution, compute p(x, u)
• RJMCMC - define move distribution and auxiliary distribution and reversible transform, compute p(x, u)

Of course, the user can a simply specify a proposal manually, which is just a function from ω to ω' , under some constraints:

• ω' should be valid in the sense that X(ω') is well defined.
  • But what is X?
    • I
  • And what if ω has components that aren't necessary for X in the sense that there is some restriction of ω such that X(ω) is well defined.
    • In finite dimensional models, these extra components won't have an effect. Or at least their effect will be averaged out
    • In trans dimensional models its a little less clear. Need to think about it

We need to be able to compute p(ω)

• If all the exogenous variables are accounted for, then we can simply just compute those
• There may be cases where we don't want to go back to the exogenous variables, e.g.:

One idea to specify proposals is to define prehooks

In the simplest case we might decide that:

1. logpdf(\omega) is just the logpdf of the sum of the exogenous variables
2. Proposals only happen by changing exogenous variables

Consequently, there is no need to resimulate the model after making a proposal, or even the need to consider the model at all

The problems with this are that:

1. If a dynamical model, changing the value of an exogenous variable can remove or add new variables, we need some mechanism (e.g., simulating the model) to find out what these new variables are so that the can be included in the new logpdf
2. If the model is conditioned, then changing one exogenous variable my determine the value of another, which would affect the logpdf, we need resimulation to find this out.

So from this perspective we can think of resimulation as taking a modified ω' and

1. completing it, so that x(ω') is well defined in that
   a. it is complete, and not missing components,
   b. it is free from inconsistency
2. and potentially that it is minimal as described above

[zenna]

For example, suppose we have the following normal-normal model:

μ ~ Normal(0, 1)
x ~ Normal(μ, 1)
μ | x = 5

In OmegaCore, you would express this as:

μ = 1 ~ Normal(0, 1)
x = 2 ~ Normal(μ, 1)
x_ = 0.123
μc = μ |ᶜ x ==ₚ x_

There is some amount of syntactic sugar going on here. The functions constructed in explicit form would be:

using Distributions: Normal, quantile, cdf

"std-normal inverse cdf"
q(x) = quantile(Normal(0, 1), x)

# Let's assume ω : Int -> [9, 1]
stdunif(id, ω) = ω[id]

n(id, ω, μ_, σ_) = q(stdunif(id, ω)) *  σ_ + + μ_

μ(ω) = n(1, ω, 0, 1)
x(ω) = n(2, ω, μ(ω), 1)
μc(ω) = x(ω) == x_ ? μ(ω) : error()

I'll consider different ways

Conventional way

In the conventional way, we:

• Make proposals over μ
• We compute the pdf using the chain rule p(μ = μ_, x = 5) = p(μ = μ_)p(x = 5 | μ = μ_)
• The proposal probability q(μ' = μnew_ | μ = μold_) is whatever it is
• Don't consider internals of μ or x

Old Omega Way

• Make proposal over stdunifs
• Compute logpdf over values of stdunifs (which is unnecessary is all uniform
• Compute soft predicate as quasilikelihood
• proposals are independent on exogenous

New Omega Way 1

• Make proposal on stdunif1 or stdunif2
• Compute pdf as p(u1 = u1_, u2 = u2)
• u2 is function of u1, in particular:
  x = 0.123 = n(2, w, mu(w), 1) where mu(w) = n(1, w, 0, 1) = q(stdunif(id, ω)) = q(u1), substituting:

x_ = n(2, w, q(u1), 1) = q(u2) + q(u1)

Hence q(u2) = x_ - q(u1) and u2 = invq(x - q(u1))

• To get the value of u2 to compute the pdf we'd need to do this, but if the base space were uniform it would be unnecessary
• We would have to compute the proposal density though, using the change of variables

Third Option

This option is a hybrid of the two:

• We find q(u2) as a function of u1
  • Then the logpdf would be p(u1 = u1_, qu2 = qu2)
  • in other words p(x = 5, u1 = u1) = p(x = 5 | u1 = u1)p(u1 = u1)