def nlsq_fit(x, y, dy, func, params_init, verbose=False, **kwargs):
    """Perform a non-linear least squares fit

    Inputs:
        x: one-dimensional numpy array of the independent variable
        y: one-dimensional numpy array of the dependent variable
        dy: absolute error (square root of the variance) of the dependent
            variable. Either a one-dimensional numpy array or None. In the array
            case, if any of its elements is NaN, the whole array is treated as
            NaN (= no weighting)
        func: a callable with the signature
            func(x,par1,par2,par3,...)
        params_init: list or tuple of the first estimates of the
            parameters par1, par2, par3 etc. to be fitted
        `verbose`: if various messages useful for debugging should be printed on
            stdout.

        other optional keyword arguments will be passed to leastsq().

    Outputs: p, dp, statdict where
        p: list of fitted values of par1, par2 etc.
        dp: list of estimated errors
        statdict: dictionary of various statistical parameters:
            'DoF': Degrees of freedom
            'Chi2': Chi-squared
            'Chi2_reduced': Reduced Chi-squared
            'R2': Coefficient of determination
            'num_func_eval': number of function evaluations during fit.
            'func_value': the function evaluated in the best fitting parameters
            'message': status message from leastsq()
            'error_flag': integer status flag from leastsq() ('ier')
            'Covariance': covariance matrix (variances in the diagonal)
            'Correlation_coeffs': Pearson's correlation coefficients (usually
                denoted by 'r') in a matrix. The diagonal is unity.

    Notes:
        for the actual fitting, scipy.optimize.leastsq() is used.
    """
    if verbose:
        t0 = time.monotonic()
        print("nlsq_fit starting.")
    else:
        t0 = 0
    func_orig = func
    params_init_orig = params_init
    func, params_init = hide_fixedparams(func_orig, params_init_orig)
    if (dy is None) or (dy == np.nan).sum() > 0 or (dy <= 0).sum() > 0:
        if verbose:
            print("nlsq_fit: no weighting")
        dy = None

    def objectivefunc(params, x, y, dy):
        """The target function for leastsq()."""
        if dy is None:
            return (func(x, *(params.tolist())) - y)
        else:
            return (func(x, *(params.tolist())) - y) / dy
    # do the fitting
    if verbose:
        print("nlsq_fit: now doing the fitting...")
        t1 = time.monotonic()
    else:
        t1 = 0
    par, cov, infodict, mesg, ier = leastsq(objectivefunc,
                                            np.array(params_init),
                                            (x, y, dy), full_output=True,
                                            **kwargs)
    if verbose:
        print("nlsq_fit: fitting done in %.2f seconds." % (time.monotonic() - t1))
        print("nlsq_fit: status from scipy.optimize.leastsq(): %d (%s)" % (ier, mesg))
        print("nlsq_fit: extracting statistics.")
    # test if the covariance was singular (cov is None)
    if cov is None:
        cov = np.ones((len(par), len(par))) * np.nan  # set it to a NaN matrix
    # calculate the Pearson's R^2 parameter (coefficient of determination)
    if dy is None:
        sserr = np.sum(((func(x, *(par.tolist())) - y)) ** 2)
        sstot = np.sum((y - np.mean(y)) ** 2)
    else:
        sserr = np.sum(((func(x, *(par.tolist())) - y) / dy) ** 2)
        sstot = np.sum((y - np.mean(y)) ** 2 / dy ** 2)
    r2 = 1 - sserr / sstot
    # assemble the statistics dictionary
    statdict = {'DoF' : len(x) - len(par),  # degrees of freedom
                'Chi2' : (infodict['fvec'] ** 2).sum(),
                'R2' : r2,
                'num_func_eval' : infodict['nfev'],
                'func_value' : func(x, *(par.tolist())),
                'message' : mesg,
                'error_flag' : ier,
               }
    statdict['Chi2_reduced'] = statdict['Chi2'] / statdict['DoF']
    statdict['Covariance'] = cov * statdict['Chi2_reduced']
    par, statdict['Covariance'] = resubstitute_fixedparams(par, params_init_orig, statdict['Covariance'])
    # calculate the estimated errors of the fit parameters
    dpar = np.sqrt(statdict['Covariance'].diagonal())
    # Pearson's correlation coefficients (usually 'r') in a matrix.
    statdict['Correlation_coeffs'] = statdict['Covariance'] / np.outer(dpar,
                                                                       dpar)
    if verbose:
        print("nlsq_fit: returning with results.")
        print("nlsq_fit: total time: %.2f sec." % (time.monotonic() - t0))
    return par, dpar, statdict