3.1.1 returns ValueError

[sarahshi]

Hi Patricio, how are you. I have developed a new package, PyIRoGlass, for fitting baselines to FTIR spectra with MC3. I have created a corresponding Google Colab for running this code on the cloud. The Google Colab version appears to have stopped working as I am now getting this error:

ValueError: Cannot populate an initial sample set of parameters, try updating the parameters initial guess to avoid sampling beyond the parameter boundaries or where the model returns non-finite values.

Do you know what changes in 3.1.1 could be returning this new error? Nothing in my code has changed with the initial set of parameters. The local version of the code also remains functional, as I am running 3.0.13 and have not updated to 3.1.1 yet. I have used the 'trf' option in the mc3.sample function and all the least squares best-fitting parameters are further within my range of initials. Here is what is returned from least squares:

Least-squares best-fitting parameters:
[1.5850219, -0.0199404429, 0.711353268, -0.229082749, 
-0.00998274823, 1428.13013, 30.0187476, 0.108494754,
 1517.87931, 34.8462181, 0.110622216, 0.673443334,
 0.0359489381, 0.0949624095, 0.000125732721, 1.25690699]

compared against what I provide as initials:

params = np.array([1.25, 2.00, 0.25, 0.01, 0.01, 1430, 25.0, 0.0100, 1510, 25.0, 0.0100, 0.10, 0.02, 0.01, 5e-4, 0.70])
pmin = np.array([0.00, -5.00, -1.00, -0.75, -0.75, 1415, 22.5, 0.0000, 1500, 22.5, 0.0000, 0.00, -0.50, -0.50, -5e-2, -1.00])
pmax = np.array([5.00, 8.00, 1.00, 0.75, 0.75, 1445, 40.0, 3.0000, 1535, 40.0, 3.0000, 3.00, 0.50, 0.50, 5e-2, 3.00])
pstep = np.array([0.30, 0.50, 0.20, 0.20, 0.20, 3.25, 2.25, 0.0005, 6.0, 2.25, 0.0005, 0.25, 0.75, 0.75, 0.002, 0.20])

Any insight you have into this would be immensely helpful!

[pcubillos]

Hi Sarah, Thanks for the heads up! I think I know what's going on.

One of the changes implemented between 3.0.13 and 3.1.1 was how the code draws samples before the MCMC to have a starting population. Before, the code was simply drawing random samples around your initial point (in this case the best-fit pars) and if any value was out of ranges, the value was clipped into the (pmin, pmax) boundaries. This worked OK, but it's sub-optimal (particularly if the starting point is close to the edge of the parameter space).

In the new version, when the code tries these initial draws it checks if any of the values is out of bounds, and if so, it throws that away and tries a new draw, repeating until we have a enough samples. The problem is that the code only makes a certain number of attempts at creating this starting sample, if not, it gives up and throws the error you are seeing.

There are two things that could help:

• (if possible and if it makes sense for your case) one option is to reduce the pstep value for some of the parameters (as these values are the width of the Gaussian used to draw this starting sample). For example, for a couple of parameters you have, the pstep is on the same order of magnitude of the [pmin-pmax] range, so it is expected that many draws will fall outside the wanted range.

• On my end, I will increase the threshold of how many attempts the code makes to initialize this sample (the reason I don't have this in a while-True loop is because the code could end up hanging indefinitely, which happened a few times for the worst case scenarios when I was editing the code). This should be simple, so, it'll be done by tomorrow.

I hope that helped and that is clear enough!

[sarahshi]

Hi Patricio,

Glad to hear about this update -- it sounds super useful. It definitely makes sense to drop my pstep value; this is a bit of an oversight on my end. I think I set these values up at some point and didn't modify them too much. I'll change my setup so stepping can be something like 5% of the absolute value of the total pmin-pmax range.

Best, Sarah

[sarahshi]

It looks like reducing the pstep to 5% of absolute value of the total pmin-pmax range is still causing this ValueError. I released a testPyPi package without mc3 version specifications, so it automatically installs 3.1.1.

[pcubillos]

Ahh, OK, I think I found the issue: when you make the mc3.sample() call in PyIRoGlass/__init__.py (around line 920), you missed the pstep argument. When I changed to something like this below, the run does not break:

    mc3_output = mc3.sample(
        data=data, uncert=uncert, func=func,
        params=params, indparams=indparams,
        pstep=pstep, pmin=pmin, pmax=pmax,
        # No prior?
        priorlow=priorlow, priorup=priorup,
        pnames=pnames, texnames=texnames, sampler='snooker',
        nsamples=1e6, nchains=9, ncpu=3, burnin=5000, thinning=5,
        leastsq='trf', grtest=True, grbreak=1.01, grnmin=0.5,
        plots=False, log=log, savefile=savefile,
    )

In this case, mc3 is guessing the stepsizes to a tenth of the input param value (so the mu's and sigma parameter were throwing most of the trials outside the boundaries)

[pcubillos]

Note also that you are setting some Gaussian priors (nonzero priorlow and priorup), but you are not setting the mean of the Gaussian (prior), so, the code is ignoring these priors.

[sarahshi]

Hi Patricio, thank you so much for your response. I have this working on my end now with 3.1.1, but I'm finding that the algorithm appears to be quite a bit slower. Is this because I'm setting a smaller pstep than the former automated value? I've tested this both with and without specifying priors (as I lacked before), to see if this could have impacted things.

I am finding that the posteriors appear to capture some baselines that appear to be like outliers. I've tested this with essentially no priors and all uniform distributions (as you saw my code had before), but with the reduced pstep of 5% of the range between np.abs(pmin-pmax). This is my old chunk of code to plot some of the sampled baselines in the posteriors.

                posteriorerror = np.load(path_beg+savefilepath+files+'.npz')
                samplingerror = posteriorerror['posterior'][:, 0:5]
                samplingerror = samplingerror[0: np.shape(posteriorerror['posterior'][:, :])[0] :int(np.shape(posteriorerror['posterior'][:, :])[0] / 200), :]
                lineerror = posteriorerror['posterior'][:, -2:None]
                lineerror = lineerror[0:np.shape(posteriorerror['posterior'][:, :])[0]:int(np.shape(posteriorerror['posterior'][:, :])[0] / 200), :]
                Baseline_Array = np.array(samplingerror * PCAmatrix[:, :].T)
                Baseline_Array_Plot = Baseline_Array

                ax4 = plt.subplot2grid((2, 3), (0, 2), rowspan = 2)
                for i in range(np.shape(Baseline_Array)[0]):
                    Linearray = Linear(Wavenumber, lineerror[i, 0], lineerror[i, 1])
                    Baseline_Array_Plot[i, :] = Baseline_Array[i, :] + Linearray
                    plt.plot(Wavenumber, Baseline_Array_Plot[i, :], 'dimgray', linewidth = 0.25)

This would return something like this attached figure. AC4_EUH102_030920_256s_15x20_a_old.pdf

With 3.1.1 and using the same MCMC function, I instead get this resulting figure. AC4_EUH102_030920_256s_15x20_a_new.pdf

I have altered the code to remove the burn-in with zmask, to see if that was causing this.

posteriorerror = np.load(path_beg+savefilepath+'AC4_EUH102_030920_256s_15x20_a.npz')
posterior = posteriorerror['posterior']
mask = posteriorerror['zmask']
masked_posterior = posterior[mask]
samplingerror = posterior[mask][:, 0:5]
samplingerror = samplingerror[0: np.shape(masked_posterior[:, :])[0] :int(np.shape(masked_posterior[:, :])[0] / 200), :]
lineerror = posteriorerror['posterior'][:, -2:None]
lineerror = lineerror[0:np.shape(masked_posterior[:, :])[0]:int(np.shape(masked_posterior[:, :])[0] / 200), :]
Baseline_Array = np.array(samplingerror * PCAmatrix[:, :].T)
Baseline_Array_Plot = Baseline_Array

plt.figure(figsize = (12, 12))
for i in range(np.shape(Baseline_Array)[0]):
    Linearray = pig.Linear(Wavenumber, lineerror[i, 0], lineerror[i, 1])
    Baseline_Array_Plot[i, :] = Baseline_Array[i, :] + Linearray
    plt.plot(Wavenumber, Baseline_Array_Plot[i, :], 'dimgray', linewidth = 0.25)

Could the pstep be what is causing this change in behavior, or would you expect that this is more related to 3.1.1 updates?

[pcubillos]

Hi Sarah, yes you have to remove the burn-in samples because you have to let the MCMC converge first. How many samples to burn is a bit of a manual process because it differs for every different problem. You can check this by looking at the trace plots (basically you want to see that the trace is not moving up/down over time, or that there are not too many outliers).

About the performance, what times are you seeing now? and before? I ran your Google Colab example, and it took about a minute per MCMC. This sounds reasonable to me. If the run takes more than a few minutes, one way to quantify this is to compare the MCMC vs a forward-model run of your model. Say, for example, run your model in a loop for 100 times, and calculate how much time it takes per evaluation. This should be about the same order of magnitude of an MCMC's runtime / (nsamples/ncpu)

[sarahshi]

Hi Patricio. I think I have resolved some of the slower behavior now by optimizing a few other things on my end.

Thank you for all of your hard work on the package. I am enjoying the new plotting with the stated parameters with the +/- values. I was wondering if it would be possible to add a fig.suptitle command into the mix for the histogram, pairwise, trace, and modelfit figures? I have done this in my local version of the old plotting code. In the case of spectral fitting where I'm working with hundreds of spectra, it would be helpful to have the ability to have the sample name as a title in the output figures.

I have also tested some of the pairwise plots -- is there a way to reduce the linewidth of the contourf lines? I am noticing that your example in the documentation here shows thinner lines, and this would be helpful for visualization. Not sure if this is purely a result of the size of this figure or the number of parameters.

AC4_EUH102_030920_256s_15x20_a_pairwise.pdf

[pcubillos]

Hi Sarah, OK, good to hear things are working better now. I will close this issue then.

• The problem with the linewidths is fixed when you update to 3.1.2.
• About the subtitles, unfortunately the code does not do that, but you can re-create the plots and add a subtitle manually with these lines after you call the MCMC routines:

    ifree = mc3_output['ifree']
    posterior, zchain, zmask = mc3.utils.burn(mc3_output)
    post = mc3.plots.Posterior(
        posterior, pnames=mc3_output['texnames'][ifree],
        bestp=mc3_output['bestp'][ifree],
    )
    savefile_root = os.path.splitext(savefile)[0]

    # Pairwise posteriors:
    fig = post.plot()
    plt.suptitle('PLOT1 SUBTITLE')
    savefile = f'{savefile_root}_pairwise_posterior.png'
    plt.savefig(savefile, dpi=300)

    # Histograms:
    post.plot_histogram()
    plt.suptitle('PLOT2 SUBTITLE')
    savefile = f'{savefile_root}_marginal_posterior.png'
    plt.savefig(savefile, dpi=300)