[tl;dr here is a collection of resources for SDEMEMs]
Mixed-effects models (MEM) are hierarchical models suited for “population inference”, where instead of fitting data from a single experiment, we are interested in learning characteristics common to 
In MEM, individual experiments (e.g. “subjects”) are modeled by introducing some set of parameters 
So we have mixed-effects, as some parameters vary randomly between subjects (
Example:
where 
When the 
then we have obtained a stochastic differential equation mixed-effects model (SDEMEM). Some papers consider SDEMEMs written as (2) (no measurement error) others consider it as the system (1)-(2).
This is a powerful class of models, since we can simultaneously consider three sources of variability: (i) variation between subjects: by estimating the parameter 
Now, inference for SDEMEMs is not trivial at all, essentially because inference for SDEs based on discrete observations can be quite tricky. However some literature is available, and the motivation for writing this post is to make reasearchers aware of a collection of resources for SDEMEMs I have set up.
Please contact me to help with the list, if you find some reference is missing.
. Suppose that 
. Finally, suppose that a non-negative and “unbiased” random estimator of 
. Then an intuitive thing to do is to use 
into a 
becomes (see 
is accepted with probability 
. Although it is not evident yet, what we have just written is an example of a pseudo-marginal method, but I will clarify this point below. Therefore in practice applying the method is relatively simple.
instead of the true and unavailable acceptance ratio 
(where 
) — has 
the generated chain is sampled from 
the number of MCMC iterations.
for every 
denotes expectation. However which probability measure do we consider for this expectation? Answering this question shows why the topic of this post is particularly relevant for sampling from intractable distributions or, in Bayesian inference, sampling from posterior distributions having intractable likelihoods.
we are not really interested in (we call it a “disturbance parameter” or an “auxiliary variable”). Let’s 
be their joint distribution. Assume the domain of the marginalized 
iid samples 
from the distribution of 
at every ![{\pi}(x)=\int {\pi}(x|\xi)p(\xi)d\xi = \mathbb{E}_\xi [{\pi}(x|\xi)]](https://s0.wp.com/latex.php?latex=%7B%5Cpi%7D%28x%29%3D%5Cint+%7B%5Cpi%7D%28x%7C%5Cxi%29p%28%5Cxi%29d%5Cxi+%3D+%5Cmathbb%7BE%7D_%5Cxi+%5B%7B%5Cpi%7D%28x%7C%5Cxi%29%5D&bg=ffffff&fg=333333&s=0&c=20201002)
and the notation 
emphasizes this fact.
we have shown that ![\mathbb{E}_\xi [{\pi}(x|\xi^n)] = \pi(x)](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D_%5Cxi+%5B%7B%5Cpi%7D%28x%7C%5Cxi%5En%29%5D+%3D+%5Cpi%28x%29&bg=ffffff&fg=333333&s=0&c=20201002)
and therefore by denoting with 
and![\mathbb{E}_\xi [\hat{\pi}(x)]=\mathbb{E}_\xi [\frac{1}{N}\sum_{n=1}^N {\pi}(x|\xi^n)]=\frac{1}{N}\sum_{n=1}^N \mathbb{E}_\xi [{\pi}(x|\xi^n)]=\pi(x)](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D_%5Cxi+%5B%5Chat%7B%5Cpi%7D%28x%29%5D%3D%5Cmathbb%7BE%7D_%5Cxi+%5B%5Cfrac%7B1%7D%7BN%7D%5Csum_%7Bn%3D1%7D%5EN+%7B%5Cpi%7D%28x%7C%5Cxi%5En%29%5D%3D%5Cfrac%7B1%7D%7BN%7D%5Csum_%7Bn%3D1%7D%5EN+%5Cmathbb%7BE%7D_%5Cxi+%5B%7B%5Cpi%7D%28x%7C%5Cxi%5En%29%5D%3D%5Cpi%28x%29&bg=ffffff&fg=333333&s=0&c=20201002)
. Also, the estimate is clearly non-negative.
for a parameter 
, given data 
. From Bayes’ theorem we have 
, with 
the prior of 
, 
the likelihood function and 
the “evidence” (if you are interested in state-space models you might want to check an 
. We show that in this case it is also possible to obtain 
from 
, an “auxiliary” variable we do not care for. Also in this case we consider sampling from an “augmented” space, which is the space 
over which the pair 
is defined.
even if we are ultimately interested in 
.
.
bacause of the assumed a-priori independence between 
.
at the denominator: this follows from the unbiasedness assumption as we have 
and the latter is by definition the evidence 
and then we plug the evidence into 
.
not 
.
. This move is generated from a proposal kernel 
. We write the proposed value 
or rejected with probability 
where
and write (assuming an arbitrary initial state 
)
compute 
accept 
, set 
otherwise store draws[r]:= 
end
means “copy the value at the right-hand-side into the left-hand-side”. After a few (sometimes many) iterations the chain reaches stationarity, and its stationary distribution happens to have 
as probability density (or probability mass function). In the end we will explore the content of 
for rejecting). So our decision is a random variable taking two possible values, hence is a Bernoulli random variable with “parameter” 
from a uniform distribution on the interval 
, that is 
which can be obtained in R via 
, and (ii) if 
accept 
then accept 
hence whenever 
is larger than 1, then certainly 
.
accept 
and 
loop it is fairly likely you may have acted the wrong way). In the implementation given in the introduction to this post (the one using the 
loop) you may notice that when a proposal is rejected I store the last accepted value as ![\texttt{draws[r]}:= x_{r-1}](https://s0.wp.com/latex.php?latex=%5Ctexttt%7Bdraws%5Br%5D%7D%3A%3D+x_%7Br-1%7D&bg=ffffff&fg=333333&s=0&c=20201002)
.
, so that we do not need to recompute it at the next iteration (since it is reused at the denominator of 
such that 
for some “prior” 
and “likelihood function” 
. In fact, it is often the case that the pointwise evaluation of the likelihood is far from being instantaneous, e.g. because the data has thousands of values, or the evaluation of the likelihood itself requires some numerical approximation.
your software will likely return zero (this is an underflow).
and you’ll see it returns zero .
accept 
without doing your calculations analytically first. That is, I see beginners trying to use the log-domain suggestion from point 4 to evaluate 
. This will not help because the log will be applied to a function 
to compute 
which is the safe way to evaluate the (normalised) standard Gaussian log-density. For other languages you may need to code the log-density yourself.
. This brings us to the next point.
but can get rid of the normalization constant 
.
then 
as part of the proportionality constant. When 
is not an integer but is a positive real number the evaluation of 
for every 
, that is a Gaussian density with mean 
. Then we have
denotes is the determinant of a matrix and 
denotes transposition. Clearly this 
is symmetric because of the quadratic form at the exponent.
is chosen by the experimenter and 
is sampled from a standard Gaussian. Therefore we have that 
. Clearly 
since we are assuming that the support of the targeted 
instead of 
. Then we evaluate 
, that is we exponentiate the proposal 
(which is Gaussian) and then by exponentiating we effectively sampled 
(which is log-normal).
we have (recall in this example 
.
. Hence the correct sampler does the following: simulate a 
is an invertible function. Say that for a generic 
we seek for the distribution of 
, and using the transformation rule linked above we have that 
. Notice here 
is the “Jacobian of the transformation”.
in the acceptance probability is the ratio of Jacobians.
(not-a-number) or that you do not get, say, 
(yes these things can happen if you code 
will evaluate to a logical 
statement and the proposed 
at the starting value 
(or 
) and we just put this value at the denominator of the ratio in the definition of 
for all iterations and the sampler will never move away from 
the importance weight for particle 
we have
could be plugged into a Bayesian procedure for sampling from the 
(more on this in future posts).
is problematic when particle 
produces a very small weight 
the datum 
might be a realization from a tail of the distribution with density 
. While this is legitimate, numerically it can create issues for the stability of our calculations. Say that 
, that is a Gaussian distribution with mean 
and standard deviation 0.5. Let’s see what happens if I anticipate some topic treated later on.
, then the weight becomes
then most of the particles 
generated from 
should take values not too far from 30. And here is a problem: from the figure we have 
and the computer implementation for the density function of 
might underflow and return zero for many particles (I will describe some tricks mitigating this problem in a future post). Clearly, all those particles having zero weight doom the values of the descendant particles, because each weight depends on the previous one. A similar scenario could happen even without outliers, if we simulate from a model not appropriate for the data. Clearly, as the time horizon 
.
as the probability associated to 
. If this can help, imagine the particles to be balls contained in an urn: some balls are large (large 
. Basically, we empty the urn then fill it up again with copies of the balls having the recorded indices. Say that we have extracted index 
five times, we put in the urn five copies of the ball with index 
(the resampling has destroyed the information on “how” we reached time 
).
) is illustrated in a 
, that is we propagate forward the state of the system by simulating particles 
. And how do we perform this step? We take advantage of the important particles we have just resampled, using the importance density to compute the move 
. For illustration see the picture below.
depicted with a curve. The particles weight is computed and you can see that some particles have a larger size (large weight), others are smaller and some have “disappeared”, which means that their weight is negligible. At this point the particles are resampled: notice at the resampling stage the largest particle in the figure happens to be resampled three times while others fewer times. Once resampling is performed, all the resampled particles get the same size because all weights are reset to 
as described above. Now it is time to propagate the resampled particles forward to time 
.
means that these are “blind” to data, since this importance density is unconditional to 
. Hence the propagation step does not take advantage of any information carried by the data. In some scenario this produces an inefficient sampling that may result, again, in particle degeneracy. A popular alternative is the 
(initialize) 
and assign 
, for all 
.
.
.
, for all 
and normalise weights 
.
and if 
go to step 2 otherwise go to 7.
.
is the unconditional density of an arbitrary initial state 
. The two formulations are completely equivalent as they only differ by a constant which is irrelevant for the purpose of assigning a weight to a particle. Also, since weights are going to be normalised, the 
. In step 7 I have explicitly reported the likelihood approximation, even though for parameter inference the product of the normalizing constants 
can be disregarded.
. If the degeneracy at 
) then resampling is performed, otherwise no resampling is performed at time 
we want to approximate the likelihood function
-dimensional integral. We can approach the problem from different angles, all interconnected.
as
with respect to the law 
. This means writing 
. The interpretation of an integral as a probabilistic expectation is at the core of 
, 
.
. Notice the conditioning on the sampled 
and the approximation improves for increasing 
.
regardless the dimension of 
gets evaluated, we want those particles that are “important” for the evaluation of the integral, and we wish most of the 
to end up in a region of high probability for 
is an arbitrary density function called “importance density” that is non-zero whenever 
is non-zero (the support of 
appearing in the penultimate equality. Therefore we should choose a 
easy to simulate from.
and 
. The first result is from the Markov property defined in the 
from the “history of 
prior to time 
independently, 
.
.
.
with 
. Clearly, generating at each time 
spanning the entire interval ![[0,t]](https://s0.wp.com/latex.php?latex=%5B0%2Ct%5D&bg=ffffff&fg=333333&s=0&c=20201002)
is not really computationally appealing, and it is not even clear how to construct the importance density. We need some strategy enabling a sort of sequential mechanism. Moreover, even if we are able to simulate the particles, are we able to evaluate all the densities appearing in 
. This is addressed in next post.
-dimensional integral bearable. Here I am going to show how to obtain the weights update above.
and 
:
depends on measurements up to time 
depends on measurements up to 
. This is because we have freedom in designing 
does not have to result in 
, though this choice is allowed. For simplicity, in the following I remove subscripts for 
.
, instead we can perform computations 
as we do not need to carry the constant term 
resulting from the denominator of the Bayes theorem. This term is independent of 
and the “evidence” 
,where the latter is really independent of 
, where 
denotes the state of the system at a time 
. The notation 
is latent (unobservable), and we can only get to know something about the system via some noisy measurements. I denote the available measurements (data) with 
and use 
to denote the process producing the actual observations at 
. Each 
) but it does not need to have the same dimension of the corresponding 
are affected with some measurement noise. For example we could have 
, with 
some random noise.
can be either continuous– or discrete–valued stochastic processes. However in the following I assume both processes to be defined on continuous spaces.
to denote a sequence 
. Therefore, data can be written 
to denote realizations of the process at times 
in place of 
, for 
. That is, “given the present state, the past is independent of the future”, so if time 
is the “present”, then 
. Also in this case, for simplicity we assume implicit the conditioning on 
. Specifically for our inference goals, we are interested in transitions between states corresponding to contiguous (integer) observational times, that is 
.
. That is, 
.
and 
and evidently it results in 
. If equality did hold 
would be independent of 
.]
. By convention we can set this arbitrary time to be 
. Therefore
is the likelihood function of 
.
.
(
) compute the sample means of the marginals, quantiles etc. This way we perform uncertainty quantification for all components of 
as a product of unconditional densities: instead we have
.
to 
are unknown.
Umberto Picchini's research blog 