I would value your opinion on the following piece of code. I am rather new to both Python and Monte Carlo analysis, so I was wondering whether the routine makes sense to more experienced and knowledgeable users.

`def MC_analysis_a(): x = spin_lock_durations y_signal_a = (a_norm1, a_norm2, a_norm3, a_norm4, a_norm5, a_norm6, a_norm7, a_norm8) x = np.array(x, dtype = float) y_signal_a = np.array(y_signal_a, dtype = float) def func(x, a, b): return a * np.exp(-b * x) initial_guess = [1.0, 1.0] fitting_parameters, covariance_matrix = optimize.curve_fit(func, x, y_signal_a, initial_guess) print(round(fitting_parameters[1], 2)) # ---> PRODUCING PARAMETERS ESTIMATES total_iterations = 5000 MC_pars = np.array([]) for iTrial in range(total_iterations): xTrial = x yTrial = y_signal_a + np.random.normal(loc = y_signal_a, scale = e_signal_a, size = np.size(y_signal_a)) try: iteration_identifiers, covariance_matrix = optimize.curve_fit(func, xTrial, yTrial, initial_guess) except: dumdum = 1 continue # ---> STACKING RESULTS if np.size(MC_pars) < 1: MC_pars = np.copy(iteration_identifiers) else: MC_pars = np.vstack((MC_pars, iteration_identifiers)) # ---> SLICING THE ARRAY print(np.shape(MC_pars)) # print(np.median(aFitpyars[:,1])) print(np.std(MC_pars[:,1])) `

The output I get is apparently satisfactory and plausible.

Many thanks in advance to any contributor!