Python scipy.optimize.fmin() Examples

def get_min(self, x0=None, **kwds):
        """Minimize fit function by `scipy.optimize.fmin()`.

        Parameters
        ----------
        x0 : 1d array, optional
            Initial guess. If not given then `points[i,...]` at the min of
            `values` is used.
        **kwds : keywords to fmin()

        Returns
        -------
        1d array (ndim,)
        """
        _kwds = dict(disp=0, xtol=1e-12, ftol=1e-8, maxfun=1e4, maxiter=1e4)
        _kwds.update(kwds)
        if x0 is None:
            idx = self.values.argmin()
            x0 = self.points[idx,...]
        xopt = optimize.fmin(self, x0, **_kwds)
        return xopt


# Need to inherit first Fit1D such that Fit1D.get_min() is used instead of
# PolyFit.get_min(). 

def gammamomentcond2(distfn, params, mom2, quantile=None):
    '''estimate distribution parameters based method of moments (mean,
    variance) for distributions with 1 shape parameter and fixed loc=0.

    Returns
    -------
    difference : array
        difference between theoretical and empirical moments

    Notes
    -----
    first test version, quantile argument not used

    The only difference to previous function is return type.

    '''
    alpha, scale = params
    mom2s = distfn.stats(alpha, 0.,scale)
    return np.array(mom2)-mom2s



######### fsolve doesn't move in small samples, fmin not very accurate 

def test_neldermead_initial_simplex_bad(self):
        # Check it fails with a bad simplices
        bad_simplices = []

        simplex = np.zeros((3, 2))
        simplex[...] = self.startparams[:2]
        for j in range(2):
            simplex[j+1,j] += 0.1
        bad_simplices.append(simplex)

        simplex = np.zeros((3, 3))
        bad_simplices.append(simplex)

        for simplex in bad_simplices:
            if self.use_wrapper:
                opts = {'maxiter': self.maxiter, 'disp': False,
                        'return_all': False, 'initial_simplex': simplex}
                assert_raises(ValueError,
                              optimize.minimize, self.func, self.startparams, args=(),
                              method='Nelder-mead', options=opts)
            else:
                assert_raises(ValueError, optimize.fmin, self.func, self.startparams,
                              args=(), maxiter=self.maxiter,
                              full_output=True, disp=False, retall=False,
                              initial_simplex=simplex) 

def _compute_reg_bw(self, bw):
        if not isinstance(bw, string_types):
            self._bw_method = "user-specified"
            return np.asarray(bw)
        else:
            # The user specified a bandwidth selection method e.g. 'cv_ls'
            self._bw_method = bw
            res = self.bw_func[bw]
            X = np.std(self.exog, axis=0)
            h0 = 1.06 * X * \
                 self.nobs ** (- 1. / (4 + np.size(self.exog, axis=1)))

            func = self.est[self.reg_type]
            bw_estimated = optimize.fmin(res, x0=h0, args=(func, ),
                                         maxiter=1e3, maxfun=1e3, disp=0)
            return bw_estimated 

def _cv_ls(self):
        r"""
        Returns the cross-validation least squares bandwidth parameter(s).

        Notes
        -----
        For more details see pp. 16, 27 in Ref. [1] (see module docstring).

        Returns the value of the bandwidth that maximizes the integrated mean
        square error between the estimated and actual distribution.  The
        integrated mean square error (IMSE) is given by:

        .. math:: \int\left[\hat{f}(x)-f(x)\right]^{2}dx

        This is the general formula for the IMSE.  The IMSE differs for
        conditional (``KDEMultivariateConditional``) and unconditional
        (``KDEMultivariate``) kernel density estimation.
        """
        h0 = self._normal_reference()
        bw = optimize.fmin(self.imse, x0=h0, maxiter=1e3, maxfun=1e3, disp=0,
                           xtol=1e-3)
        bw = self._set_bw_bounds(bw)  # bound bw if necessary
        return bw 

def get_min(self, x0=None, **kwds):
        """Return [x,y] where z(x,y) = min(z) by minimizing z(x,y) w/
        scipy.optimize.fmin().

        Parameters
        ----------
        x0 : sequence, length (2,), optional
            Initial guess. If None then use the data grid point with the
            smallest `z` value.

        Returns
        -------
        [xmin, ymin]: 1d array (2,)
        """
        _kwds = dict(disp=0, xtol=1e-12, ftol=1e-8, maxfun=1e4, maxiter=1e4)
        _kwds.update(kwds)
        if x0 is None:
            idx0 = self.values.argmin()
            x0 = [self.xx[idx0], self.yy[idx0]]
        xopt = optimize.fmin(self, x0, **_kwds)
        return xopt 

def test_minimize_l_bfgs_b_maxfun_interruption(self):
        # gh-6162
        f = optimize.rosen
        g = optimize.rosen_der
        values = []
        x0 = np.ones(7) * 1000

        def objfun(x):
            value = f(x)
            values.append(value)
            return value

        # Look for an interesting test case.
        # Request a maxfun that stops at a particularly bad function
        # evaluation somewhere between 100 and 300 evaluations.
        low, medium, high = 30, 100, 300
        optimize.fmin_l_bfgs_b(objfun, x0, fprime=g, maxfun=high)
        v, k = max((y, i) for i, y in enumerate(values[medium:]))
        maxfun = medium + k
        # If the minimization strategy is reasonable,
        # the minimize() result should not be worse than the best
        # of the first 30 function evaluations.
        target = min(values[:low])
        xmin, fmin, d = optimize.fmin_l_bfgs_b(f, x0, fprime=g, maxfun=maxfun)
        assert_array_less(fmin, target) 

def validate_models(self, n_users, n_debug=None):
        df_train, df_val = self.load_train_val(n_users, n_debug=n_debug)

        preds_mat = np.vstack([model.fit_and_predict(df_train, df_val, validate=True) for model in self.models]).T

        def opt_coefs(coefs):
            preds = preds_mat.dot(coefs)
            df_val["preds"] = preds
            mrr = mrr_fast(df_val, "preds")
            print(mrr, coefs)
            return -mrr

        best_coefs = fmin(opt_coefs, [model.weight for model in self.models])
        best_coefs = fmin_powell(opt_coefs, best_coefs)

        preds = preds_mat.dot(best_coefs)
        df_val["click_proba"] = preds
        print("MRR {:4f}".format(mrr_fast(df_val, "click_proba")))
        print("Best coefs: ", best_coefs) 

def compute_reference_frame(self,epipole,F):
		'''
		Method to compute the reference frame of the reconstruction (i.e. plane at infinity in an affine or metric space).
		Args: 
			epipole: the epipole
			F: the fundamental matrix
		Returns:
			p: the reference plane
			h: the homography [e]xF

		'''
		H=self._eg_utils.compute_homography(epipole,F)	#compute the homography [e]xF 
		# get the reference plane 
		p = np.sum(np.divide(np.eye(3)-H, np.transpose(np.asarray([epipole, epipole, epipole]))),axis=0)/3
		# adjust reference plane to make the first two projection matrices as equal as possible
		p=fmin(self.init_plane,np.append(p,1),xtol=1e-25,ftol=1e-25,args=(H.real,epipole.real));
		p=p[0:3]
		return p, H 

def Get_RpRs(d, **kwargs):
    '''
    Returns the value of the planet radius over the stellar radius
    for a given depth :py:obj:`d`, given
    the :py:class:`everest.pysyzygy` transit :py:obj:`kwargs`.

    '''
    if ps is None:
            raise Exception("Unable to import `pysyzygy`.")

    def Depth(RpRs, **kwargs):
        return 1 - ps.Transit(RpRs=RpRs, **kwargs)([kwargs.get('t0', 0.)])

    def DiffSq(r):
        return 1.e10 * (d - Depth(r, **kwargs)) ** 2

    return fmin(DiffSq, [np.sqrt(d)], disp=False) 

def Get_rhos(dur, **kwargs):
    '''
    Returns the value of the stellar density for a given transit
    duration :py:obj:`dur`, given
    the :py:class:`everest.pysyzygy` transit :py:obj:`kwargs`.

    '''
    if ps is None:
            raise Exception("Unable to import `pysyzygy`.")

    assert dur >= 0.01 and dur <= 0.5, "Invalid value for the duration."

    def Dur(rhos, **kwargs):
        t0 = kwargs.get('t0', 0.)
        time = np.linspace(t0 - 0.5, t0 + 0.5, 1000)
        try:
            t = time[np.where(ps.Transit(rhos=rhos, **kwargs)(time) < 1)]
        except:
            return 0.
        return t[-1] - t[0]

    def DiffSq(rhos):
        return (dur - Dur(rhos, **kwargs)) ** 2

    return fmin(DiffSq, [0.2], disp=False) 

def gammamomentcond2(distfn, params, mom2, quantile=None):
    '''estimate distribution parameters based method of moments (mean,
    variance) for distributions with 1 shape parameter and fixed loc=0.

    Returns
    -------
    difference : array
        difference between theoretical and empirical moments

    Notes
    -----
    first test version, quantile argument not used

    The only difference to previous function is return type.

    '''
    alpha, scale = params
    mom2s = distfn.stats(alpha, 0.,scale)
    return np.array(mom2)-mom2s



######### fsolve doesn't move in small samples, fmin not very accurate 

def fit_minic(self):
        #this doesn't make sense, since number of parameters stays unchanged
        #need leave-one-out, gcv; or some penalization for weak priors
        #added extra penalization for lambd
        def get_bic(lambd):
            #return self.fit(lambd).bic #+lambd #+ 1./lambd  #added 1/lambd for checking
            #return self.fit(lambd).gcv()
            #return self.fit(lambd).cv()
            return self.fit(lambd).aicc()

        from scipy import optimize
        lambd = optimize.fmin(get_bic, 1.)
        return lambd

#TODO:
#I need the hatmatrix in the model if I want to do iterative fitting, e.g. GCV
#move to model or use it from a results instance inside the model,
#    each call to fit returns results instance 

def _compute_reg_bw(self, bw):
        if not isinstance(bw, string_types):
            self._bw_method = "user-specified"
            return np.asarray(bw)
        else:
            # The user specified a bandwidth selection method e.g. 'cv_ls'
            self._bw_method = bw
            res = self.bw_func[bw]
            X = np.std(self.exog, axis=0)
            h0 = 1.06 * X * \
                 self.nobs ** (- 1. / (4 + np.size(self.exog, axis=1)))

            func = self.est[self.reg_type]
            bw_estimated = optimize.fmin(res, x0=h0, args=(func, ),
                                         maxiter=1e3, maxfun=1e3, disp=0)
            return bw_estimated 

def isf(STAT, alpha, df, resels, n, Q=None, version='spm12'):
	'''
	Inverse survival function
	'''
	if isinstance(alpha, (int,float)):
		alpha       = [alpha]
	zstar = []
	for aaa in alpha:
		z0    = _approx_threshold(STAT, aaa, df, resels, n)
		fn    = lambda x : (rft(1, 0, STAT, x[0], df, resels, n, Q, False, version)[0] - aaa)**2
		zzz   = optimize.fmin(fn, z0, xtol=1e-9, disp=0)[0]
		zstar.append(zzz)
	return np.asarray(zstar)




################################
# Convienence classes
################################ 

def expected_loss_for_range(estimate_range, true_s, function='absolute', u=1,o=1,u_f=1,o_f=1, r=1):
    """Function to attain expected loss and Bayes action based on an absolute-error or
           squared error-loss function for a defined range of estimates.

                Args:
                    estimate_range (np.array): Range of value estimates.
                    true_s (np.array): Array of possible true value occurrences (from a probability distribution)
                    u (int or float, optional): Underestimation re-weighting factor.
                    o (int or float, optional): Overestimation re-weighting factor.
                    u_f (int or float, optional): Fatal underestimation re-weighting factor.
                    o_f (int or float, optional): Fatal overestimation re-weighting factor.
                    r (int, float or np.array, optional): Risk-affinity re-weighting factor.
                Returns:
                    [0]: Expected loss for the range of estimates.
                    [1]: Bayes action (estimate with minimal expected loss)
                    [2]: Expected loss of Bayes action.
        """
    expected_loss_s = lambda estimate_s, r: expected_loss_for_estimate(estimate_s, true_s, function, u, o, u_f, o_f, r)
    loss_e = [expected_loss_s(e, r) for e in estimate_range]
    bayes_action = sop.fmin(expected_loss_s, -40, args=(r,), disp=False)
    bayes_action_loss_e = expected_loss_s(bayes_action, r)
    return loss_e, bayes_action, bayes_action_loss_e 

def _find_threshold(errors, z_range):
    """Find the ideal threshold.

    The ideal threshold is the one that minimizes the z_cost function. Scipy.fmin is used
    to find the minimum, using the values from z_range as starting points.

    Args:
        errors (ndarray):
            Array of errors.
        z_range (list):
            List of two values denoting the range out of which the start points for the
            scipy.fmin function are chosen.

    Returns:
        float:
            Calculated threshold value.
    """
    mean = errors.mean()
    std = errors.std()

    min_z, max_z = z_range
    best_z = min_z
    best_cost = np.inf
    for z in range(min_z, max_z):
        best = fmin(z_cost, z, args=(errors, mean, std), full_output=True, disp=False)
        z, cost = best[0:2]
        if cost < best_cost:
            best_z = z[0]

    return mean + best_z * std 

def _cv_ml(self):
        """
        Returns the cross validation maximum likelihood bandwidth parameter.

        Notes
        -----
        For more details see p.16, 18, 27 in Ref. [1] (see module docstring).

        Returns the bandwidth estimate that maximizes the leave-out-out
        likelihood.  The leave-one-out log likelihood function is:

        .. math:: \ln L=\sum_{i=1}^{n}\ln f_{-i}(X_{i})

        The leave-one-out kernel estimator of :math:`f_{-i}` is:

        .. math:: f_{-i}(X_{i})=\frac{1}{(n-1)h}
                        \sum_{j=1,j\neq i}K_{h}(X_{i},X_{j})

        where :math:`K_{h}` represents the Generalized product kernel
        estimator:

        .. math:: K_{h}(X_{i},X_{j})=\prod_{s=1}^
                        {q}h_{s}^{-1}k\left(\frac{X_{is}-X_{js}}{h_{s}}\right)
        """
        # the initial value for the optimization is the normal_reference
        h0 = self._normal_reference()
        bw = optimize.fmin(self.loo_likelihood, x0=h0, args=(np.log, ),
                           maxiter=1e3, maxfun=1e3, disp=0, xtol=1e-3)
        bw = self._set_bw_bounds(bw)  # bound bw if necessary
        return bw 

def loss_for_range(estimate_range, true_s, function='absolute', u=1,o=1,u_f=1,o_f=1, r=1):
    """Function to attain loss and Bayes action based on an absolute-error or
               squared error-loss function for a defined range of estimates
               given one single true value.

                    Args:
                        estimate_range (np.array): Range of value estimates.
                        true_value (int or float): Array of possible true value occurrences (from a probability distribution)
                        u (int or float, optional): Underestimation re-weighting factor.
                        o (int or float, optional): Overestimation re-weighting factor.
                        u_f (int or float, optional): Fatal underestimation re-weighting factor.
                        o_f (int or float, optional): Fatal overestimation re-weighting factor.
                        r (int, float or np.array, optional): Risk-affinity re-weighting factor.
                    Returns:
                        [0]: Loss for the range of estimates.
                        [1]: Bayes action (estimate with minimal expected loss)
                        [2]: Expected loss of Bayes action.

                        Note: Since the is only one possible true value,
                        the Bayes action will always equal this value.
            """
    incurred_loss = lambda estimate_s, r: loss_for_estimate(estimate_s, true_s, function, u, o, u_f, o_f, r)
    loss_i = [incurred_loss(e, r) for e in estimate_range]
    bayes_action = sop.fmin(incurred_loss, -40, args=(r,), disp=False)
    bayes_action_loss = incurred_loss(bayes_action, r)
    return loss_i, bayes_action, bayes_action_loss 

def fitquantilesgmm(distfn, x, start=None, pquant=None, frozen=None):
    if pquant is None:
        pquant = np.array([0.01, 0.05,0.1,0.4,0.6,0.9,0.95,0.99])
    if start is None:
        if hasattr(distfn, '_fitstart'):
            start = distfn._fitstart(x)
        else:
            start = [1]*distfn.numargs + [0.,1.]
    #TODO: vectorize this:
    xqs = [stats.scoreatpercentile(x, p) for p in pquant*100]
    mom2s = None
    parest = optimize.fmin(lambda params:np.sum(
        momentcondquant(distfn, params, mom2s,(pquant,xqs), shape=None)**2), start)
    return parest 

def beta_opt(self, method='mxie', apply=False, **kwargs):
        if method == 'mxie':
            beta_orig_x, beta_orig_y = self.betax, self.betay
            beta_orig = np.mean([beta_orig_x, beta_orig_y])
            
            fel_copy = deepcopy(self)
            def f(x, method=method):
                fel_copy.betax = fel_copy.betay = x
                fel_copy.eval(method=method)
                return fel_copy.lg3
            
            err_dict = np.geterr()
            np.seterr(all='ignore')
            beta_opt = fmin(f, beta_orig, disp=0, **kwargs)
            np.seterr(**err_dict)
            
        elif method == 'ssy_opt':
            beta_opt = self.beta_opt_calc
        
        else:
            _logger.error('method should be in ["mxie", "ssy_opt"]')
            raise ValueError('method should be in ["mxie", "ssy_opt"]')
            
        if apply:
            self.betax = beta_opt
            self.betay = beta_opt
            self.eval(method)
        else:
            return beta_opt[0] 

def fit(self, start_params=None, maxiter=5000, method='fmin', tol=1e-08):
        if start_params is None:
            start_params = np.concatenate((0.05*np.ones(self.nar + self.nma), [1]))
        mlefit = super(Arma, self).fit(start_params=start_params,
                maxiter=maxiter, method=method, tol=tol)
        return mlefit 

def fit(self, start_params=None, maxiter=5000, method='fmin', tol=1e-08):
        '''estimate model by minimizing negative loglikelihood

        does this need to be overwritten ?
        '''
        if start_params is None and hasattr(self, '_start_params'):
            start_params = self._start_params
        #start_params = np.concatenate((0.05*np.ones(self.nar + self.nma), [1]))
        mlefit = super(TSMLEModel, self).fit(start_params=start_params,
                maxiter=maxiter, method=method, tol=tol)
        return mlefit 

def _fit_nm(f, score, start_params, fargs, kwargs, disp=True,
                maxiter=100, callback=None, retall=False,
                full_output=True, hess=None):
    xtol = kwargs.setdefault('xtol', 0.0001)
    ftol = kwargs.setdefault('ftol', 0.0001)
    maxfun = kwargs.setdefault('maxfun', None)
    retvals = optimize.fmin(f, start_params, args=fargs, xtol=xtol,
                            ftol=ftol, maxiter=maxiter, maxfun=maxfun,
                            full_output=full_output, disp=disp, retall=retall,
                            callback=callback)
    if full_output:
        if not retall:
            xopt, fopt, niter, fcalls, warnflag = retvals
        else:
            xopt, fopt, niter, fcalls, warnflag, allvecs = retvals
        converged = not warnflag
        retvals = {'fopt': fopt, 'iterations': niter,
                   'fcalls': fcalls, 'warnflag': warnflag,
                   'converged': converged}
        if retall:
            retvals.update({'allvecs': allvecs})
    else:
        xopt = retvals
        retvals = None

    return xopt, retvals 

def fmin_pos(func, x0, *args, **kwargs):
  '''fmin with positivity constraint''' 
  def pos_func(x, *blargs):
    return func(np.exp(x), *blargs)

  out = fmin(pos_func, np.log(x0), *args, **kwargs)
  out = np.exp(out)
  return out 

def _cv_ml(self):
    bw = self._normal_reference() # use normal_reference as initialization for bandwidth
    x0 = np.concatenate([bw, np.array([self.epsilon])], axis=0)

    def loologli(x):
      assert x.shape[0] == self.ndim_y + 1
      bw = x[:self.ndim_y]
      eps = float(x[self.ndim_y])
      return - self.loo_likelihood(bw, epsilon=eps)

    x_opt = optimize.fmin(loologli, x0=x0, maxiter=_MAX_ITER_CV_ML_OPTIMIZER, disp=0)
    bw_opt, eps_opt = x_opt[:self.ndim_y], x_opt[self.ndim_y]
    return bw_opt, eps_opt 

def _fit_nm(f, score, start_params, fargs, kwargs, disp=True,
                maxiter=100, callback=None, retall=False,
                full_output=True, hess=None):
    xtol = kwargs.setdefault('xtol', 0.0001)
    ftol = kwargs.setdefault('ftol', 0.0001)
    maxfun = kwargs.setdefault('maxfun', None)
    retvals = optimize.fmin(f, start_params, args=fargs, xtol=xtol,
                            ftol=ftol, maxiter=maxiter, maxfun=maxfun,
                            full_output=full_output, disp=disp, retall=retall,
                            callback=callback)
    if full_output:
        if not retall:
            xopt, fopt, niter, fcalls, warnflag = retvals
        else:
            xopt, fopt, niter, fcalls, warnflag, allvecs = retvals
        converged = not warnflag
        retvals = {'fopt': fopt, 'iterations': niter,
                   'fcalls': fcalls, 'warnflag': warnflag,
                   'converged': converged}
        if retall:
            retvals.update({'allvecs': allvecs})
    else:
        xopt = retvals
        retvals = None

    return xopt, retvals 

def fit(self, correct_gear, gears, engine_speeds_out, times, velocities,
            accelerations, motive_powers, velocity_speed_ratios, stop_velocity):
        from ..mechanical import calculate_gear_box_speeds_in
        self.clear()
        self.velocity_speed_ratios = velocity_speed_ratios
        self.update(_identify_gear_shifting_velocity_limits(
            gears, velocities, stop_velocity
        ))
        if dfl.functions.CMV.ENABLE_OPT_LOOP:
            from co2mpas.utils import mae
            gear_id, velocity_limits = zip(*sorted(self.items())[1:])
            max_gear, _inf = gear_id[-1], float('inf')
            update, predict = self.update, self.predict

            def _update_gvs(vel_limits):
                self[0] = (0, vel_limits[0])
                self[max_gear] = (vel_limits[-1], _inf)
                update(dict(zip(gear_id, _grouper(vel_limits[1:-1], 2))))

            def _error_fun(vel_limits):
                _update_gvs(vel_limits)

                g_pre = predict(
                    times, velocities, accelerations, motive_powers,
                    correct_gear=correct_gear
                )

                speed_pred = calculate_gear_box_speeds_in(
                    g_pre, velocities, velocity_speed_ratios, stop_velocity
                )
                return np.float32(mae(speed_pred, engine_speeds_out))

            x0 = [self[0][1]].__add__(
                list(itertools.chain(*velocity_limits))[:-1]
            )
            from scipy.optimize import fmin
            _update_gvs(fmin(_error_fun, x0, disp=False))

        return self 

def modelNFR(self, boundaries = (35,115)):
        """Model NFR distribution with gamma distribution"""
        b = np.where(self.fragmentsizes.get(self.lower,boundaries[1]) == max(self.fragmentsizes.get(self.lower,boundaries[1])))[0][0] + self.lower
        boundaries = (min(boundaries[0],b), boundaries[1])
        x = np.arange(boundaries[0],boundaries[1])        
        y = self.fragmentsizes.get(boundaries[0],boundaries[1]) 
        def gamma_fit(X,o,p):
            k = p[0]
            theta = p[1]
            a = p[2]
            x_mod = X-o
            res = np.zeros(len(x_mod))
            if k>=1:
                nz = x_mod >= 0
            else:
                nz = x_mod > 0
            res[nz] = a * x_mod[nz]**(k-1) * np.exp(-x_mod[nz]/theta) / (theta **k * gamma(k))
            return res 
        res_score = np.ones(boundaries[0]+1)*np.float('inf')
        res_param = [0 for i in range(boundaries[0]+1)]
        pranges = ((0.01,10),(0.01,150),(0.01,1))
        for i in range(15,boundaries[0]+1):
            f = lambda p: np.sum((gamma_fit(x,i,p) - y)**2)
            tmpres = optimize.brute(f, pranges,  full_output=True,
                              finish=optimize.fmin)
            res_score[i] = tmpres[1]
            res_param[i] = tmpres[0]
        whichres = np.argmin(res_score)
        res = res_param[whichres]
        self.nfr_fit0 = FragmentSizes(self.lower,self.upper, vals = gamma_fit(np.arange(self.lower,self.upper),whichres,res_param[whichres]))
        nfr = np.concatenate((self.fragmentsizes.get(self.lower,boundaries[1]), self.nfr_fit0.get(boundaries[1],self.upper))) 
        nfr[nfr==0] = min(nfr[nfr!=0])*0.01
        self.nfr_fit = FragmentSizes(self.lower,self.upper, vals = nfr)
        nuc = np.concatenate((np.zeros(boundaries[1]-self.lower),
                            self.fragmentsizes.get(boundaries[1],self.upper) -
                            self.nfr_fit.get(boundaries[1],self.upper)))
        nuc[nuc<=0]=min(min(nfr)*0.1,min(nuc[nuc>0])*0.001)
        self.nuc_fit = FragmentSizes(self.lower, self.upper, vals = nuc) 

def fit_mps(dist, data, x0=None):
    '''Estimate distribution parameters with Maximum Product-of-Spacings

    Parameters
    ----------
    params : array_like, tuple ?
        parameters of the distribution funciton
    xsorted : array_like
        data that is already sorted
    dist : instance of a distribution class
        only cdf method is used

    Returns
    -------
    x : ndarray
        estimates for the parameters of the distribution given the data,
        including loc and scale


    '''
    xsorted = np.sort(data)
    if x0 is None:
        x0 = getstartparams(dist, xsorted)
    args = (xsorted, dist)
    print(x0)
    #print(args)
    return optimize.fmin(logmps, x0, args=args)