Let $f$ be a function such that:

$$f~:~(x,~\theta)\in\mathbb{R}^{3}\times\mathbb{R}^{12} \rightarrow f(x,~\theta)\in\mathbb{R}^3$$

My observations $y$ are noisy values taken by the function $f(\cdot ,~\theta)$ for known values of $x$. I would like to estimate the conditional distribution of $\theta \in \mathbb{R}^{12}$ given the observations $y$. I assume the prior distribution of $\theta$ is uniform, and the likelihood of the observations given $\theta$ is Gaussian.

I would like to have some insights about the estimation of $\theta$ with MCMC Metropolis-Hastings simulations. I use Gaussian proposals. To test the algorithm, I have created synthetic noiseless $y$ given a known $\theta$, and I try to use MCMC to find point-estimates really close to the true $\theta$. To make things harder, I start far from the true values, but even if I started "close" (in $\mathbb{R}^{12}$) to the solution, the convergence is not fast or obvious. I guess this is due to the quite complex non linear expression of $f$.

1) First of all, since $y\in\mathbb{R}^3$, would it be relevant to define separate Gaussian likelihood functions, one for each component of $y$, with potentially different error variances $\sigma_i^2$. Or is it better to just define one error variance $\sigma^2$ and jointly consider all the components of $y$ ?

2) Second, if I want to provide an estimate of $\sigma^2$ to the algorithm, I could just compute the empirical variance, which is roughly $\frac{1}{n-p}$ times the sum of square differences between $y$ and $f(x,~\theta_{init})$. However, I have noticed, by looking at the each simulated component of $\theta$ independently, that the chain does not mix well unless I give a much higher value to $\sigma^2$ (like omitting the factor $\frac{1}{n-p}$ where $n$ can be of the order of $10^5$ and $p=12$). Would this be due to the high dimensionality of the estimation problem? As I understand it, the higher $\sigma^2$, the higher the probability to accept a new move even if the move "looks bad", and thus the higher the acceptance rate. So, to put my second question in other words, should I be patient and drastically increase the number of simulations, or should I just force the chain to sample more space more quickly by increasing $\sigma^2$? I would love to get a result in $10^5$ simulations, but I am afraid I will need a lot more because of the dimensionality and non linearity.

  • 1
    $\begingroup$ What do your acceptance rates look like? $\endgroup$ – jerad Feb 26 '13 at 22:14
  • $\begingroup$ If I initialize $\sigma^2$ with the empirical variance of the error with $\theta=\theta_{init}$, the acceptance rate is about 1%. $\endgroup$ – Wok Feb 26 '13 at 23:05

An acceptance rate of $1\%$ is extremely low. You have to aim for an acceptance rate of around $15\%-30\%$. In the normal case the optimal rate is the famous $0.234$. However, this is a tricky diagnostic tool since this "optimal rate" can also be achieved despite a terrible sampling (for example when you sample well from an entry and terribly from another one, or some other combinations). You should also use other diagnostic tools, such as the log-posterior, autocorrelation plots, among others to check the performance of your MCMC.

In addition, it is recommended to use (if possible) an adaptive method as a second opinion. (If you are using R, Python, or C, take a look at these packages)

  • $\begingroup$ Thanks. I have kept the empirical variance for $\sigma^2$, but I have managed to quicken the convergence by using another covariance matrix for the proposal. $\endgroup$ – Wok Feb 27 '13 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.