Convergence of the sampler

[Peku995]

Hi,

I have a question about the behaviour of the sampler. I am currently doing simulations on scattering data and let the sampler start at the point estimation of a least-square simulation. Interestingly, the behaviour is as follows: In the attached picture you can see the sample log-function and as you can see, the sampler runs away from the optimum and instead settles at unfavourable values and does not reach the optimum anymore.
Is there an explanation for this behaviour? Maybe I misunderstood the concept, but I thought that the sampler should sample values with highter probability more often?

Interestingly, the adaption measure becomes better with increasing steps.

[shahmoradi]

Hi @Peku995, From the bounces in the shape of the adaptation measure, I can tell that you have likely set a range for targetAcceptanceRate.
That initial drop from the highest likelihood is entirely natural and indicates your problem likely has more than a few dimensions. The function value never returns to its original maximum at the starting point because it is extremely unlikely to find it again. The bulk of probability is around that peak or possibly on a mountain ridge. The fact that the likelihood function value remains stable is a good sign, indicating that your objective function implementation is numerically correct. Is this objective function constructed to fit a dataset? How many dimensions does data have?

Sometimes, the objective function is degenerate, like a mountain ridge (instead of a mountain peak). There could be a minor peak at which a deterministic optimizer would settle. But when you start exploring the space via a random sampler, the sampler is free to take jumps of any size out of the local peak and, therefore, readily gets out of the local minimum and falls on the ridge.

[Peku995]

Thank you for the answer and the explanations!
Yes, the objective function was programmed to fit measurement data, it has 7 dimensions.
The strange thing is that the majority of the data points are not clustered around the maximum likelihood estimator of the least square procedure but run relatively far away from it. I have inserted an image. The MLE is actually 1, but you can see that the frequency is more like 1.2. If I use the 1.2 to fit my measurement data, then you can also see that the result does not fit the measurement data.

I really don't understand that behaviour.

[shahmoradi]

also, do parameters converge to particular values throughout the simulation?

[shahmoradi]

Is the above figure a histogram of likelihoods sampled?
If so, consider the following function with a minimal peak,

Even though the function has a clear mode, a sampler will rarely visit (and sample) the function mode merely because it is probabilistically insignificant; that is, there is little area underneath the mode. Instead, the bulk of probability (the area under the function) is around the function mode. Therefore, the sampler will spend most of its time around the mode and not on it. If you were to draw the histogram of the sampled likelihoods (function values) for this objective function, you'd see the mode barely visited or perhaps never visited. But that is natural and correct. What matters is the average of a population (average likelihood here), not the best in a population (the function mode here). Now, as soon as you enter dimensions more than 3, it is almost impossible ever to sample the mode of a function because the odds of visiting a mode are practically zero or near zero.
I hope this helps!
p.s. What I find more disturbing than the rare mode visits is that the sampler can explore higher function values than the initial mode the optimizer found to be around f(x) = 1, unless the MLE you report in your comment is just an approximate os is incorrect.
I recommend starting your optimizer with the most likely parameter set you get from the sampling (which should be near MLE 1.4). This MLE value and corresponding parameters are reported toward the end of the output *_report.txt file. Then see if the optimizer gets back to the original mode it found earlier or if it converges to the same mode as the sampler finds.
Deterministic optimizers can only find local extrema. This may also be the case here.

[Peku995]

Thank you for the helpful hints!
However, I think I have expressed myself incorrectly. The histrogram shows the value that the sampler has determined for a variable over a chain with a length of 50000. The least-square method returns 1 as the optimal value for this variable, yet the sampler does not return most values near 1, but at 1.15. The MCMC also returns about 1 as the optimal value, which can be seen in the scatterplot below. What just puzzles me is that most of the values found are still at 1.15, so the values are well away from 1. And as can be seen in the histrogram above, one would clearly estimate the most probable value at 1.15 and definitely not at 1.0 because this looks not just like a broad distribution but rather a distinct peak.

[shahmoradi]

@Peku995 The objective function you are dealing with reminds me of the K2 mountains, a higher-peak mountain on the shoulder of another lower-profile mountain.

Flip the x and y axes of the K2 mountain picture, and we get the objective function along the axis you have plotted above and the histogram in your earlier post.
Looking at the adaptation measure you posted above, there is also a possibility that the sampler has not converged yet. I'd suggest switching off the targetAcceptanceRate and adjusting the scaleFactor up or down if the average efficiency reported on the screen is above 20% or too much below 10%. You may have to do this rescaling for a few short simulations. Ideally, you would trust the chain to be Markovian when the adaptation measure falls consistently below 0.001.
Also, posting a trace plot of the variable above, like the logFunc trace plot in your original question above, can help determine whether there is a lack of convergence. If there is a lack of convergence, the output sample in *_sample.txt would contain too few points. But you say you have 5000 final refined samples, which means the algorithm has likely converged. If there is no sign of a lack of convergence, then the results are expected, and there is a possibility that there are more complexities in other dimensions that are not visualized and not evident in the scatter plot above. If there is a lack of convergence, normally the sampler prints a warning message about the output sample size on screen and in the *_report.txt file.

[Peku995]

Thanks again for your comment! I have now tried what you suggested and not specified the targetAcceptanceRate but have set the scaleFactor instead. Now I also achieve a lower value for the AdaptionMeasure, see image.

However, the same problem still appears and even more clearly: The sampler runs away from the optimal solution and never samples in this range again, see picture of the log function.

I have to say that this is the chain and not the sample, however, the sample data looks similar and I got approx. 900 samples from a chain size of 50000.
I had a problem like this once before with other measurement data, but I was able to solve it by choosing a different starting point. That is not the case here. Do you have maybe any other ideas what I could try?

[shahmoradi]

Hi @Peku995, This logfunc trace plot looks natural to me, but of course, it implies that there is a tiny region of parameter space where the model provides ~ exp(3) times better fit to data. But the fact that the sampler quickly gets out of it indicates the isolated mode is not probabilistically significant. That's the whole point of sampling and uncertainty quantification, that is, to get out of the mode and learn the average behavior of the model. It initially stays in the mode for a long time because the sampler has not mixed well yet; that is, the initial visits around the mode are not i.i.d. samples from the likelihood, and the adaptation measure trace plot you posted corroborates this interpretation.
I suggest using the best-fit MLE value in the model reports in the *_report.txt file and reporting the 90% confidence band associated with each parameter mode. If you want to ensure the sampler's convergence (which seems achieved in your traceplot), you can start it from different locations in the parameter space and see if they all converge to the sample neighborhood in the parameter space.

In you subsequent sampling, you can set startPointVec to the value obtained from the previous simulation in the report file named stats.chain.refined.avgStd and set the proposalStartCovMat to the corresponding value
stats.chain.refined.covmat in the report file from the previous simulation. This way, the sampler won't waste time warming up starting from ground zero every time.
You can also run multiple chains in parallel to see if they all converge to the same neighborhood, a strong indicator of convergence.
If you still doubt the model's accuracy, you can share your model and data with me directly ([email protected]) to inspect this further. I am fond of examining pathological cases in sampling problems.