Ytdyklly

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[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 = 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.

Many thanks in advance to any contributor!

edit: see comments.
edit2: typo.

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[1], 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):