Recently I need to use Bayesian method to estimate the parameters of a chemical kinetic model which can be described by ordinary differential equations. I found Turing.jl or DiffeqBayes.jl might be the correct way to do this job.

The kinetic model has four variables, however, they are not experimentally determined at the same time. For example, at time point ti, variable x1 and x3 have values but variable x2 and x4 only have missing data. In this case, how can I do Bayesian estimations of ODE parameters?

This should work. To fit data, the loss function should be defined as the sum of weighted squared errors. The weighted is the reciprocal of the variance of the measurement.

Loss function and maximum likelihood observables are effectively the same thing. Just do like:

for i in 1:length(predicted)
if i in values_for_dataset1
data1[i] ~ MvNormal(predicted[1:2, i], σ^2 * I)
else
data2[i] ~ MvNormal(predicted[3:4, i], σ^2 * I)
end
end

which would split a 4 observable system to make 1:2 to data1 and 3:4 to data2. Of course there are details to fill in there but that should be a sufficient hint.