search

florianhartig / BayesianTools

General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics
https://cran.r-project.org/web/packages/BayesianTools/index.html
119 stars 29 forks source link

Unwritable Outputs For Some Random Seeds #235

Open hosseinkavianih opened 2 years ago

hosseinkavianih commented 2 years ago

Hi,

I am trying to run the DREAM(ZS) functions for a number of parameters and trying to repeat the experiment for different random seeds. When I was processing the outputs, I came to realize that for the same setup, some outputs are not being generated for different random seeds. I was wondering if you could assist me on this matter. I appreciate your help. Hossein.

woodwards commented 2 years ago

Are you using multiple threads? I run BT using the parallel package and have to use clusterSetRNGStream() to control the random seeds.

Or maybe your model crashes for certain parameter values and returns NA/NaN etc?

florianhartig commented 2 years ago

Hi @hosseinkavianih - could you provide a reproducible example, or at least a more detailed description about how this is happening?

hosseinkavianih commented 2 years ago

Thank you for your response @florianhartig. I put the script in my Bayesian Repository (https://github.com/hosseinkavianih/Bayesian). The script is TestError.R and you can find Hyfunction.R function which is being called in the script and PR.txt, E.txt, Qsyntrue.txt, and DREAMsamples2.txt as the input files. If you run the script as it is with set.seed(1) on line 72, you'll get an output file, but if you run it with set.seed(3), you'll get this error: "Error in AdaptpCR(CR, delta, lCR, settings, Npop) : object 'pCR' not found". Please let me know if you have any questions.

florianhartig commented 2 years ago

Hi Hossein, thumbs up for the reproducible example, I can reproduce the error - I will have to investigate what's going wrong, but given that >10 evaluations of the algorithm before work fine without problem, it seems likely that some kind of numerical problem is occurring that is not properly handled.

A comment independent of that: why do you set nrChains = 17? This doesn't seem to make sense. Maybe this is a misunderstanding about what this parameter controls - these are not the internal DE chains, with nrChains you are running 17 completely independent replicate runs of the entire algorithm, which is rarely sensible. The only reason to replicate is to get good convergence diagnostics, and 3-4 chains is enough for that, the rest is wasted. If you want to change the number of internal chains, you should change the starting values of the algorithm.

p.s.: if you keep message = T, you can monitor the progress of the algorithms that are evaluated. It seems the error occurs in the 15th DREAM run that is calculated.

hosseinkavianih commented 2 years ago

Hi Florian. Thank you so much for taking a look at this. Also, sorry for my late reply. For this analysis, we’re trying to see the sensitivity of the algorithm’s convergence to different hyper parameters, including the number of chains, and have sampled 2-20. We found some examples in the literature that went up to 20 chains so that was the reason that we came up with this range to study the tradeoffs.

hosseinkavianih commented 2 years ago

Hi,

I went through the source code for the difference between internal and external chains and I have a few questions. How would you define internal chains? Is there an argument within the DREAMZs function that controls that? I though DEpairs is the argument that defines number of internal chains. Before this, I didn't came across to "internal" and "external" terms in the literature, so I am not quite sure about the differences. I am trying to come up for a reasonable range for my experiments, and in some literature, number of chains even goes up to 100 and I think this would be the internal chain that you mentioned here. I'd appreciate your feedback.

florianhartig commented 2 years ago

Hi Hossein,

as I said, I believe that what you refer to when you say 20 chains work well is the internal chains of the DE algorithms (the DREAMzs is a population MCMC, so it works via running several MCMC chains in parallel and exchanging information, this is what I refer to by "internal chains"). The number of internal chains can be changed by changing the dimension of the starting conditions.

The argument nrChains changes the number of INDEPENDENT DREAMzs Algorithms that are run. So, currently, you are running 17 completely independent DREAMzs calibration, probably with 3 internal chains each.

Regarding the literature: I'm not sure that 20 internal chains is a good value / necessary for the zs version of the DREAM algorithm. When you read about 20, it may be that people referred to the standard DREAM algorithm. Given the way the memory works in the zs, the number of chains should be of lesser importance. Given the way the algorithms can be parallelized, however, it may make sense to choose the number of internal chains equal or n-1 the number of cores on your cluster node, and parallelise the algorithm (this only makes sense for models that evaluate with > 20ms or so).