Showing the error of lmfit-model center in plot

I'm working with fitting of data with the python lmfit package. I am looking at several peaks and fitting them with a pseudo-Voigt + constant model. I then extract the center position of the fitted peak and compile about a hundred of these peak positions. The hundred positions are then plotted and fitted with a sine-wave model, also with lmfit.
My issue, is that I am not able to find the information I need to show the uncertainty of the peak positions extracted from the pseudo-Voigt fitting. In many similar questions, I encounter the eval_uncertainty function also included in the lmfit package, but I am not able to use this to extract the confidence intervals of the center position, only uncertainties in the amplitude and sigma value of the peak.

The ideal result would look like the following image:
The ideal case I'm trying to achieve

And my current graph looks like this:
Current graph

Obviously I need to switch from a scatter plot to an errorbar plot, but I do not have any error-intervals to work with at the moment.

The following is the most simple workable "mini"-case I could construct, while still including my original fit-models. My original data is around 100 datapoints instead of four, and the sine-wave in this smaller case is therefore not "smooth", but this should not affect the method of introduction of errorbars.

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import PseudoVoigtModel, ConstantModel, SineModel
def Voigt_app(x,y,a,b):
psvo_mod = PseudoVoigtModel(prefix='psvo_')
pars = psvo_mod.guess(y,x=x)
con_mod = ConstantModel(prefix='con_')
pars = pars.update(con_mod.make_params(c = (y[0]+y[1])/2))
mod = psvo_mod + con_mod
out = mod.fit(y, pars, x=x)
details = out.fit_report(min_correl=0.25)
return out, details
def Sin_app(x,y):
sin_mod = SineModel(prefix='sin_')
pars = sin_mod.guess(y,x=x)
con_mod = ConstantModel(prefix='con_')
pars = pars.update(con_mod.make_params(c = np.mean(y)))
mod = sin_mod + con_mod
out = mod.fit(y,pars,x=x)
details = out.fit_report(min_correl=0.25)
return out, details
data = np.asarray(
[[14407.0 , 26989.0 , 31026.0 , 31194.0 , 31302.0 , 91996.0 , 112250.0 , 112204.0 , 113431.0 , 75097.0 , 55942.0 , 55826.0 , 56181.0 , 28266.0 , 14687.0],
[12842.0 , 13275.0 , 13581.0 , 13943.0 , 14914.0 , 20463.0 , 21132.0 , 21457.0 , 23308.0 , 32017.0 , 30808.0 , 30927.0 , 30337.0 , 21025.0 , 15216.0],
[15770.0 , 17677.0 , 19008.0 , 20911.0 , 25958.0 , 27984.0 , 28164.0 , 31161.0 , 33517.0 , 29430.0 , 28673.0 , 32725.0 , 28495.0 , 24527.0 , 25173.0],
[16299.0 , 20067.0 , 25102.0 , 34968.0 , 37757.0 , 44871.0 , 51347.0 , 60154.0 , 54985.0 , 53383.0 , 45776.0 , 40836.0 , 30816.0 , 27922.0 , 26066.0]])
a,b = 1932, 1947
centers = []
colors = ['red','black','green','blue']
fig, ax = plt.subplots()
for i in range(4):
ax.scatter(np.arange(15),data[i,:],label=f'data {i}',color=colors[i])
out, details = Voigt_app(x,data[i,:],a,b)
ax.plot(out.best_fit,label=f'fit {i}',color=colors[i])
centers.append(out.values['psvo_center'])
ax.set_title('Data points and pseudo-Voigt fits')
ax.legend()
fig, ax = plt.subplots()
ax.scatter(np.arange(4),centers,color='red',label='data')
out, details = Sin_app(np.arange(4),np.asarray(centers))
ax.plot(np.arange(4),out.best_fit,color='black',label='fit')
ax.set_title('Peak positions and sine-wave fit')
ax.legend()

From this I would like to switch out the scatter plot in the second graph with an errorbar plot with the confidence intervals of the pseudo-Voigt peak positions ("out.values['psvo_center']").

Here's how you could modify your Voigt_app function to do that:

def Voigt_app(x,y,a,b):
psvo_mod = PseudoVoigtModel(prefix='psvo_')
pars = psvo_mod.guess(y,x=x)
con_mod = ConstantModel(prefix='con_')
pars = pars.update(con_mod.make_params(c = (y[0]+y[1])/2))
mod = psvo_mod + con_mod
out = mod.fit(y, pars, x=x)
# Calculate the confidence interval for the fitted parameters
ci = lmfit.conf_interval(out, out)
print(lmfit.printfuncs.report_ci(ci))
# Extract the uncertainty for 'psvo_center'
center_err = (ci['psvo_center'][2][1] - ci['psvo_center'][1][1])/2
details = out.fit_report(min_correl=0.25)
return out, details, center_err

You can then use these uncertainties when plotting the error bars:

center_errs = []
for i in range(4):
out, details, center_err = Voigt_app(x, data[i,:], a, b)
# ...
center_errs.append(center_err)
# ...
ax.errorbar(np.arange(4), centers, yerr=center_errs, fmt='o', color='red', label='data')

If I understand your question correctly, you want to extract not only the best-fit value for the 'psvo_center' parameter for each of your fits, but also the standard error for 'psvo_center', and then use those as the data for fitting to a sinusoidal oscillation.

If that is correct, I think that what you want to do is to save the stderr attribute, by changing you code to

centers = []
stderrs = []  # < new
colors = ['red','black','green','blue']
fig, ax = plt.subplots()
for i in range(4):
ax.scatter(np.arange(15),data[i,:],label=f'data {i}',color=colors[i])
out, details = Voigt_app(np.arange(15), data[i,:], a, b)
ax.plot(out.best_fit,label=f'fit {i}',color=colors[i])
centers.append(out.params['psvo_center'].value)       # < new 
stderrs.append(out.params['psvo_center'].stderr)      # < new

You can then use that to make and error-bar plot with

ax.errorbar(np.arange(4), centers, stderrs, color='red',label='data')

You might also want to use those stderrs to weight the fit to the sinusoidal function, perhaps modifying your Sin_App to be

def Sin_app(x, y, dely):       # < new
sin_mod = SineModel(prefix='sin_')
pars = sin_mod.guess(y,x=x)
con_mod = ConstantModel(prefix='con_')
pars = pars.update(con_mod.make_params(c = np.mean(y)))
mod = sin_mod + con_mod
out = mod.fit(y, pars, x=x, weights=1./dely)    # < new
details = out.fit_report(min_correl=0.25)
return out, details

and calling that with

out, details = Sin_app(np.arange(4), np.asarray(centers), np.asarray(stderrs))

for the 4 example datasets you give, the variation in the standard-error isn't very big, so it does not make a big difference. But as you add more datasets, it might become more important. You might want to add a check that stderr is valid (like, not None and not 0) before using this "in production".

For "one more thing", you might also want to calculate the predicted uncertainty in the sinusoidal function and plot that range too. Your example has 4 variables for the sinusoidal fit and only 4 values, which results for this case in a near-perfect fit with absurdly small uncertainties on the fitted sine function. So to illustrate the approach, I fix the con_c parameter in the fit below. If fitting the results of more than 4 peak profiles, that should not be necessary. But, for illustration,
changing your Sin_app to be

def Sin_app(x, y, dely):
sin_mod = SineModel(prefix='sin_')
pars = sin_mod.guess(y,x=x)
con_mod = ConstantModel(prefix='con_')
pars = pars.update(con_mod.make_params(c = np.mean(y)))
pars['con_c'].vary = False  # needed with only 4 (x,y) values
mod = sin_mod + con_mod
out = mod.fit(y, pars, x=x, weights=1./dely)
details = out.fit_report(min_correl=0.25)
delta_fit = out.eval_uncertainty(sigma=1)   # < new
return out, details, delta_fit

and then plotting the results of using that with

out, details, delta_fit = Sin_app(np.arange(4), np.asarray(centers), np.asarray(stderrs))
ax.plot(np.arange(4),out.best_fit, 'o', color='black',label='fit')
ax.fill_between(np.arange(4), 
out.best_fit-delta_fit, out.best_fit+delta_fit,
color="#D7D7D7")
ax.set_title('Peak positions and sine-wave fit')
ax.legend()

This leads to the final plot of:

showing the "data" for the peak centers with error bars, the fitted values, and the 1-sigma range of the fit.