Python scipy.optimize 模块，least_squares() 实例源码

def least_squares_jacobian(x, u, y):
    """Jacobian for least squares. Used by scipy.optimize.least_squares."""
    J = np.empty((u.size, x.size))
    den = u ** 2 + x[2] * u + x[3]
    num = u ** 2 + x[1] * u
    J[:, 0] = num / den
    J[:, 1] = x[0] * u / den
    J[:, 2] = -x[0] * num * u / den ** 2
    J[:, 3] = -x[0] * num / den ** 2
    return J

# import md5
#
# def remove_duplicate_junctions(igs_list):
#     """Remove igs which have same junction."""
#     igs, juncs = [], []
#     md5_list = []
#     for ig in igs_list:
#         junc = extra.junction_re(ig.junction)
#         md = md5.new(junc).digest()
#         if md not in md5_list:
#             igs.append(ig)
#             juncs.append(junc)
#             md5_list.append(md)
#     return igs, juncs

def estimate(y, yl, x, poly='beta'):
    """
    Fit MIDAS model

    Args:
       y (Series): Low-frequency data
       yl (DataFrame): Lags of low-frequency data
       x (DataFrame): High-frequency regressors

    Returns:
        scipy.optimize.OptimizeResult
    """

    weight_method = polynomial_weights(poly)

    xw, w = weight_method.x_weighted(x, weight_method.init_params())

    # First we do OLS to get initial parameters
    c = np.linalg.lstsq(np.concatenate([np.ones((len(xw), 1)), xw.reshape((len(xw), 1)), yl], axis=1), y)[0]

    f = lambda v: ssr(v, x.values, y.values, yl.values, weight_method)
    jac = lambda v: jacobian(v, x.values, y.values, yl.values, weight_method)

    opt_res = least_squares(f,
                            np.concatenate([c[0:2], weight_method.init_params(), c[2:]]),
                            jac,
                            xtol=1e-9,
                            ftol=1e-9,
                            max_nfev=5000,
                            verbose=0)

    return opt_res

def doleast_squares(chivecfn, pinit, chi2args):
    from scipy.optimize import least_squares
    res = least_squares(chivecfn, pinit, method='lm',
                        x_scale='jac', args=chi2args)#, verbose=1)
    return res

# ISN'T WELL DOCUMENTED BY PROSPECTOR AND MAYBE THIS IS THE KEYYYYYYYYYYY

def doleast_squares(chivecfn, pinit, chi2args):
    from scipy.optimize import least_squares
    res = least_squares(chivecfn, pinit, method='lm',
                        x_scale='jac', args=chi2args)#, verbose=1)
    return res

# ISN'T WELL DOCUMENTED BY PROSPECTOR AND MAYBE THIS IS THE KEYYYYYYYYYYY

def solve():
    starting_params = PatientState.schnider_params()

    params_array = convert_params_structure_to_vector(starting_params)

    result = least_squares(get_residuals_for_all_patients, params_array, bounds=[-1000, 1000], ftol=1e-001, xtol=1e-001, gtol=1e-001, x_scale=1.0, max_nfev=1000, verbose=0)
    return result.x

def corner_plot(self,filter_params=(),num_samples=100000,*args,**kwargs):
        """
        Create a "corner plot" that shows distributions of values for each
        parameter, as well as cross-correlations between parameters.

        Parameters
        ----------
        filter_params : list-like
            strings used to search parameter names.  if the string matches, 
            the parameter is *excluded* from the plot.
        num_samples : int
            how many samples to generate

        Wraps the Fitter.plot_corner method.  Least squares does not generate 
        samples but corner.corner samples.  Use the Jacobian spit out by 
        least_squares to generate a whole bunch of fake samples.  These are 
        then deleted. 

        Approximate the covariance matrix as $(2*J^{T} \dot J)^{-1}$, then perform
        cholesky factorization on the covariance matrix.  This can then be
        multiplied by random normal samples to create distributions that come
        from this covariance matrix. 

        See:       
        https://stackoverflow.com/questions/40187517/getting-covariance-matrix-of-fitted-parameters-from-scipy-optimize-least-squares
        https://stats.stackexchange.com/questions/120179/generating-data-with-a-given-sample-covariance-matrix
        """

        J = self._fit_result.jac

        cov = np.linalg.inv(2*np.dot(J.T,J))
        chol_cov = np.linalg.cholesky(cov).T

        self._samples = np.dot(np.random.normal(size=(num_samples,chol_cov.shape[0])),chol_cov)
        self._samples = self._samples + self.estimate       

        fig = Fitter.corner_plot(self,filter_params,*args,**kwargs)

        del self._samples

        return fig

def __init__(self, model, tmin=None, tmax=None, noise=True, freq='D',
                 weights=None, **kwargs):
        BaseSolver.__init__(self)

        # Update the kwargs going to the solver
        self.default_kwargs = dict(ftol=1e-3)
        kwargs = self.update_kwargs(kwargs)

        parameters = model.parameters.initial.values

        # Set the boundaries
        pmin = np.where(model.parameters.pmin.isnull(), -np.inf,
                        model.parameters.pmin)
        pmax = np.where(model.parameters.pmax.isnull(), np.inf,
                        model.parameters.pmax)
        bounds = (pmin, pmax)

        if False in model.parameters.vary.values.astype('bool'):
            logger.warning("Fixing parameters is not supported with this"
                           "solver. Please use LmfitSolve or apply small"
                           "boundaries as a solution.")

        self.fit = least_squares(self.objfunction, x0=parameters,
                                 bounds=bounds,
                                 args=(tmin, tmax, noise, model, freq,
                                       weights), **kwargs)
        self.optimal_params = self.fit.x
        self.report = None

def least_squares_mdl(x, u):
    """Model for least squares. Used by scipy.optimize.least_squares."""
    return x[0] * (u ** 2 + x[1] * u) / (u ** 2 + x[2] * u + x[3])