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

def bisect(f, x0, x1, e=1e-12):
    """
    Bisection method for root finding
    """
    x = (x0 + x1) / 2.0
    n = 0
    while (x1 - x0) / 2.0 > e:
        n += 1
        if f(x) == 0:
            return x
        elif f(x0) * f(x) < 0:
            x1 = x
        else:
            x0 = x
        x = (x0 + x1) / 2.0
    print("B\ti=%d\tx=%f\tf(x)=%f" % (n, x, f(x)))
    return x

def main():
    # linear regression with gradient descent
    test_gd()
    # function optimization with gradient descent
    manual_gd(f_)
    # root finding methods
    # bisection method
    bisect(f, 0, 2)
    # secant method
    secant(f, 1, 2)
    # newton method
    newton(f, f_, 1)
    # scipy
    print("SB\tx=%f" % optimize.bisect(f, 0, 2))
    print("SS\tx=%f" % optimize.newton(f, 1))
    print("SN\tx=%f" % optimize.newton(f, 1, f_))

def BinomialErrors_old(nobs, Nsamp, alpha=0.16):
    """
    DO NOT USE THIS! Use BinomialErrors instead!

    One sided confidence interval for a binomial test.

    If after Nsamp trials we obtain nobs
    trials that resulted in success, find c such that

    P(nobs/Nsamp < mle; theta = c) = alpha

    where theta is the success probability for each trial.

    Code stolen shamelessly from stackoverflow:
    http://stackoverflow.com/questions/13059011/is-there-any-python-function-library-for-calculate-binomial-confidence-intervals
    """
    warnings.warn('Use BinomialErrors instead!', DeprecationWarning)
    from scipy.stats import binom

    p0 = float(nobs) / float(Nsamp)
    upper_errfcn = lambda c: binom.cdf(nobs, Nsamp, c) - alpha
    lower_errfcn = lambda c: binom.cdf(nobs, Nsamp, c) - (1.0 - alpha)
    return p0, bisect(lower_errfcn, 0, 1), bisect(upper_errfcn, 0, 1)

def _get_wave_number(self, k, eps_inf, mu_inf, space):
        """Calculate the wave number for auxiliary fdtd using Newton's method.

        Keyword arguments:
        k -- normalized wave vector
        eps_inf -- permittivity which fills the auxiliary fdtd
        mu_inf -- permeability which fills the auxiliary fdtd
        space -- Cartesian instance

        """
        ds = np.array(space.dr)
        dt = space.dt
        v = 1 / sqrt(eps_inf * mu_inf)
        omega = 2 * pi * self.src_time.freq
        wave_number = omega / v
        wave_vector = wave_number * np.array(k)

        zeta = bisect(self._3d_dispersion_relation, 0, 2,
                      (v, omega, ds, dt, wave_vector))

        return zeta * wave_number

def _setup(self):
        super(MinValueEntropySearch, self)._setup()

        # Apply Gumbel sampling
        m = self.models[0]
        valid = self.feasible_data_index()

        # Work with feasible data
        X = self.data[0][valid, :]
        N = np.shape(X)[0]
        Xrand = RandomDesign(self.gridsize, self._domain).generate()
        fmean, fvar = m.predict_f(np.vstack((X, Xrand)))
        idx = np.argmin(fmean[:N])
        right = fmean[idx].flatten()# + 2*np.sqrt(fvar[idx]).flatten()
        left = right
        probf = lambda x: np.exp(np.sum(norm.logcdf(-(x - fmean) / np.sqrt(fvar)), axis=0))

        i = 0
        while probf(left) < 0.75:
            left = 2. ** i * np.min(fmean - 5. * np.sqrt(fvar)) + (1. - 2. ** i) * right
            i += 1

        # Binary search for 3 percentiles
        q1, med, q2 = map(lambda val: bisect(lambda x: probf(x) - val, left, right, maxiter=10000, xtol=0.01),
                          [0.25, 0.5, 0.75])
        beta = (q1 - q2) / (np.log(np.log(4. / 3.)) - np.log(np.log(4.)))
        alpha = med + beta * np.log(np.log(2.))

        # obtain samples from y*
        mins = -np.log(-np.log(np.random.rand(self.num_samples).astype(np_float_type))) * beta + alpha
        self.samples.set_data(mins)

def find_steady_state(self, a, b, method='brentq', **kwargs):
        """
        Compute the equilibrium value of capital stock (per unit
        effective labor).

        Parameters
        ----------
        a : float
            One end of the bracketing interval [a,b].
        b : float
            The other end of the bracketing interval [a,b]
        method : str (default=`brentq`)
            Method to use when computing the steady state. Supported
            methods are `bisect`, `brenth`, `brentq`, `ridder`. See
            `scipy.optimize` for more details (including references).
        kwargs : optional
            Additional keyword arguments. Keyword arguments are method
            specific see `scipy.optimize` for details.

        Returns
        -------
        x0 : float
            Zero of `f` between `a` and `b`.
        r : RootResults (present if ``full_output = True``)
            Object containing information about the convergence. In
            particular, ``r.converged`` is True if the routine
            converged.

        """
        if method == 'bisect':
            result = optimize.bisect(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'brenth':
            result = optimize.brenth(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'brentq':
            result = optimize.brentq(self.evaluate_k_dot, a, b, **kwargs)
        elif method == 'ridder':
            result = optimize.ridder(self.evaluate_k_dot, a, b, **kwargs)
        else:
            mesg = ("Method must be one of : 'bisect', 'brenth', 'brentq', " +
                    "or 'ridder'.")
            raise ValueError(mesg)

        return result