Monte Carlo errors estimation routine

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.

import numpy as np
import scipy.optimize
from scipy.optimize import curve_fit
from scipy.stats import kde

def MC_analysis_a():
x = (0.001, 0.002, 0.004, 0.008, 0.016, 0.032, 0.064, 0.128)
y_signal_a = (1.0, 0.83, 0.71, 0.62, 0.55, 0.49, 0.45, 0.42)
x = np.array(x, dtype = float)
y_signal_a = np.array(y_signal_a, dtype = float)
e_signal_a = (0.06, 0.05, 0.04, 0.05, 0.04, 0.02, 0.04, 0.03)

def func(x, a, b):
return a * np.exp(-b * x)

initial_guess = [1.0, 1.0]
fitting_parameters, covariance_matrix = curve_fit(func, x, y_signal_a, initial_guess)
print(round(fitting_parameters, 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 = 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(MC_pars[:,1]))
print(np.std(MC_pars[:,1]))

The output I get is apparently satisfactory and plausible.

• Welcome to Code Review! In order to get the best possible outcome from the review please change your code so that is is runnable and complete (How to create a MWE?, minimal is not so important here). This allows users to verify that it is actually on-topic as well as to get a better feeling on what you actually want to accomplish. Mar 28 '19 at 17:27
• Does the code need to be generalised in any particular way? Are you wanting to apply it to different datasets, or to fit a different family of curves? Do you want to store the output for further analysis? Mar 28 '19 at 17:29
• As @Alex says. Don't omit any of your imports, especially for external references like numpy. Mar 28 '19 at 18:22
• Hello! Thanks everybody for the feedback, I'll edit the post incorporating actual data to fit. @RussHyde I don't really aim at generalisation as of now, I don't think I'll use it for anything else than exponential decay actually... Mar 28 '19 at 18:59

Before starting with feedback on the actual code, I would like to share some thoughts about the style.

Chose a style, be consistent

Python comes with an official style guide called PEP8 which summarizes a whole lot of style advice and best practices generally to be used while coding Python. If you decide to stick with PEP8 or not is up to you, but if you chose a style, e.g. for variable names, stick to it. Your code has a mix of camelCase (e.g. iTrial), snake_case (e.g. initial_guess, actually most of the code is like this and also the PEP8 recommendation) and snake_swallowed_CASE_case (e.g. MC_analysis_a).

Another recommendation of PEP8 is to omit whitespace around = when used as keyword arguments to functions/methods. Incorporating the first one the line

yTrial = y_signal_a + np.random.normal(loc = y_signal_a, scale = e_signal_a, size = np.size(y_signal_a))

would become

y_trial = y_signal_a + np.random.normal(loc=y_signal_a, scale=e_signal_a, size=np.size(y_signal_a))

The same goes for np.array(...) and friends earlier in the code.

As a final note, do your future self a favor and document what you're doing in that function in a short little description right at the beginning. You might think you will remember your thoughts when looking at it in three months, but trust me, you won't. It is as simpe as:

def mc_analysis_a():
"""Running Monte Carlo analysis on some 1D sample data

This function uses scipy.optimize.curve_fit to estimate the parameters of
a simple x -> y timeseries.
"""

This type of documentation has the nice feature that it will be picked up by Python's built-in help(...) function.

The code itself

At the moment, there are some unused imports in your code (namely import scipy.optimize and from scipy.stats import kde, but that might be a symptom from bringing the code here on Code Review.

The next thing that caught my attention was that you do not use the estimated covariance matrix returned by curve_fit. If you really have no use for it, replace covariance_matrix by _ (underscore). This a convention used to tell yourself and others that you don't care about that (part of the) return value.

Next up comes total_iterations = 5000 as hardcoded value. If you ever try to increase or decrease that value you would have to edit your function before rerunning the experiment. A more convenient approach from my point of view would be to use it as input parameter with an default value of 5000. The same applys to initial_guess and friends. (Sidenote: if you were to use initial_guess' current value as default value and ever modified it in the function, funny things can happen. See this link for more information on this topic, since this would get to involved for this review.)

Now to the core loop of the function. While generating the noisy y values, np.size(y_signal_a) should be equal to y_signal_a.size. In your case, one could also use y_signal.shape in that spot which would then guarantee that the noise also fits the shape of the array if it was multidimensional.

Using try: ... catch: ... without specifying an exception type is considered a bad practice since it will catch all kinds of exceptions including keyboard interupts (Ctrl+C). So if you ever decide to stop the test because 5000 iterations where to much, well, you would have to kill the Python process itself. The documentation of curve_fit lists all the exceptions that might be raised while using this function. These are ValueError for bad input data or input options and Runtime Error in case the least-squares minimization failed. So these would likely be your candidates to catch, with the following construct:

try:
...
catch (ValueError, RuntimeError):
continue # or handle the exception however you want

This would still stop the execution if something unexpected has happened and, as you can see, works for multiple exceptions as well.

On the same loop you're using the NumPy array MC_pars to collect the results of those iterations step by step. As you might have noticed NumPy arrays are not really comfortable to use in an expanding way since their main purpose is to allow fast vectorized operations on "fixed" size numerical data (though their use is not limited to that). At this point using a simple Python list might be a better option. This would allow you to .append(...) the results of every run very conveniently. And, as you seem to know, it's quite easy to convert Python iterables (lists, tuples) into NumPy arrays.

For reference you can find the modified version of the total script below.

import numpy as np
from scipy.optimize import curve_fit

def mc_analysis_a(total_iterations=5000):
"""Running Monte Carlo analysis on some 1D sample data

This function uses scipy.optimize.curve_fit to estimate the parameters of
a simple x -> y timeseries.
"""
x = (0.001, 0.002, 0.004, 0.008, 0.016, 0.032, 0.064, 0.128)
y_signal_a = (1.0, 0.83, 0.71, 0.62, 0.55, 0.49, 0.45, 0.42)
x = np.array(x, dtype = float)
y_signal_a = np.array(y_signal_a, dtype = float)
e_signal_a = (0.06, 0.05, 0.04, 0.05, 0.04, 0.02, 0.04, 0.03)

def func(x, a, b):
return a * np.exp(-b * x)

initial_guess = [1.0, 1.0]
fitting_parameters, _ = curve_fit(func, x, y_signal_a, initial_guess)
print(round(fitting_parameters, 2))

# ---> PRODUCING PARAMETERS ESTIMATES
mc_pars = []
for _ in range(total_iterations):   # i was replaced by _ here, because it was not used
x_trial = x
y_trial = y_signal_a + np.random.normal(loc=y_signal_a, scale=e_signal_a, size=y_signal_a.shape)
try:
iteration_identifiers, _ = curve_fit(func, x_trial, y_trial, initial_guess)
except (ValueError, RuntimeError):
continue

mc_pars.append(iteration_identifiers)

# ---> SLICING THE ARRAY
mc_pars = np.array(mc_pars)
print(mc_pars.shape)
print(mc_pars[:, 1].std())

if __name__ == "__main__":
mc_analysis_a(5000):
• Thanks @Alex for the very complete review, I was indeed unaware of the PEP8 convention. To be honest, I am pretty consistent and always use snake_case type of names, but since this code was partially adapted from a source, then some features just stayed like that. With respect to the "useless" imports, this MC simulation is part of a much larger code, approx. 1000 lines, which I thought would have been inappropriate posting. Most importantly: the estimated covariance matrix from curve_fit is needed because I am actually interested in seeing how much MC does improve the fit. Thanks a lot again! Mar 28 '19 at 21:10
• Good to hear that, I hope you'll still get some value ouf of the review. (I tend to feel uneasy when thinking about 1000+ line scripts, especially if the parts are not clearly separated. But that's a topic for another time.) Mar 28 '19 at 21:27