reciprocal function to the power n to fit convergence data
516055
def reciprocal(x, a, b, n):
    
    if n < 1:
        n = 1
    elif n > 5:
        n = 5
    if isinstance(x, list):
        y_l = []
        for x_v in x:
            y_l.append(a + b / x_v ** n)
        y = np.array(y_l)
    else:
        y = a + b / x ** n
    return y
1314882
def fpbspl(t, n, k, x, l, h, hh):
    
    h[0] = 1.0
    for j in range(1, k + 1):
        hh[0:j] = h[0:j]
        h[0] = 0.0
        for i in range(j):
            li = l+i
            f = hh[i]/(t[li] - t[li - j])
            h[i] = h[i] + f*(t[li] - x)
            h[i + 1] = f*(x - t[li - j])
775158
def generator(self, t):
        """"""
        self.check_fit()

        return (1.0 / self.theta) * (np.power(t, -self.theta) - 1)
179612
def _mse_converged(self):
        

        mse = mean_squared_error(self.local_prior, self.local_posterior_,
                                 multioutput='uniform_average')
        if mse > self.threshold:
            return False, mse
        else:
            return True, mse
502498
def _separable_approx2(h, N=1):
    
    return np.cumsum([np.outer(fy, fx) for fy, fx in _separable_series2(h, N)], 0)
382100
def p0_reciprocal(xs, ys):
    
    a0 = ys[len(ys) - 1]
    b0 = ys[0]*xs[0] - a0*xs[0]
    return [a0, b0, 1]
1035235
def fit(self, X, y, n_iter=None):
        """"""
        self.n_iter = self.n_iter if n_iter is None else n_iter
        X = getattr(X, 'values', X).reshape(len(X), 1)
        X_1 = self.homogenize(X)
        for i in range(self.n_iter):
            for i in range(0, len(X), 10):  # minibatch learning for numerical stability
                batch = slice(i, min(i + 10, len(X)))
                Xbatch, ybatch = X[batch, :], y[batch]
                X_1_batch = X_1[batch, :]
                self.W += (self.alpha / len(X) ** 1.5) * (
                    self.delta(Xbatch, ybatch).reshape((len(Xbatch), 1)).T.dot(X_1_batch))
        return self
1204957
def converge(f, step, tol, max_h):
    
    g = f(0)
    dx = 10000
    h = step
    while (dx > tol):
        g2 = f(h)
        dx = abs(g - g2)
        g = g2
        h += step

        if h > max_h:
            raise Exception("Did not converge before {}".format(h))
    return g
741166
function bgrat(a, b, x) {
      var T, fTs, u, logH, logM, lnx2s, lnx2sn, ps, cs, j, ser;

      T = a + 0.5 * (b - 1);
      fTs = 1 / (4 * T * T);
      u = -T * Math.log(x);
      logH = Math.log(b) + log1p(gam1(b)) + b * Math.log(u) - u;
      logM = logH - (algdiv(b, a) + b * Math.log(T));
      cs = [1];   // cs[n] = 1 / (2n + 1)!
      ps = [1];   // ps[n] = p_n in (9.3)
      lnx2s = Math.log(x) * 0.5;   // (ln(x) / 2) ^ 2
      lnx2s *= lnx2s;
      lnx2sn = 1;   // lnx2s ^ n
      j = gratR(b, u, logH);
      ser =
      Math.log(
         series(function(n) {
            var mbmn; // holds m*b - n

            mbmn = -n;
            if (n > 0) {
               j = (b + 2 * n - 2) * (b + 2 * n - 2 + 1) * j +
                     (u + b + 2 * n - 2 + 1) * lnx2sn;
               j *= fTs;
               cs[n] = cs[n - 1] / (2 * n) / (2 * n + 1);
               ps[n] = (b - 1) * cs[n] + 1 / n * series(function(m) {
                  if (m === 0) { return 0; }
                  mbmn += b;
                  return mbmn * cs[m] * ps[n - m];
               }, n);
               lnx2sn *= lnx2s;
            }
            return ps[n] * j;
         })
      );
      return logM + ser;
   }
197406
def fit(self, x0=None, distribution='lognormal', n=None, **kwargs):
        '''
        '''
        dist = {'lognormal': PSDLognormal, 
                'GGS': PSDGatesGaudinSchuhman, 
                'RR': PSDRosinRammler}[distribution]
        
        if distribution == 'lognormal':
            if x0 is None:
                d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                s = 0.4
                x0 = [d_characteristic, s]
        elif distribution == 'GGS':
            if x0 is None:
                d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                m = 1.5
                x0 = [d_characteristic, m]
        elif distribution == 'RR':
            if x0 is None:
                x0 = [5E-6, 1e-2]
        from scipy.optimize import minimize
        return minimize(self._fit_obj_function, x0, args=(dist, n), **kwargs)
1348820
def extrapolate_reciprocal(xs, ys, n, noise):
    
    if len(xs) > 4 and noise:
        y1 = (ys[-3] + ys[-4]) / 2
        y2 = (ys[-1] + ys[-2]) / 2
        x1 = (xs[-3] + xs[-4]) / 2
        x2 = (xs[-1] + xs[-2]) / 2
        try:
            b = (y1 - y2) / (1/x1**n - 1/x2**n)
            a = y2 - b / x2**n
        except IndexError:
            print_and_raise_error(xs, ys, 'extrapolate_reciprocal')
    else:
        try:
            b = (ys[-2] - ys[-1]) / (1/(xs[-2])**n - 1/(xs[-1])**n)
            a = ys[-1] - b / (xs[-1])**n
        except IndexError:
            print_and_raise_error(xs, ys, 'extrapolate_reciprocal')
    return [a, b, n]
1025141
def _converged(self, X):
        """"""

        if len(self.responsibilities) < 2:
            return False

        diff = np.linalg.norm(self.responsibilities[-1] - self.responsibilities[-2])
        return diff <= self.tolerance
106260
def _mm(n_items, data, initial_params, alpha, max_iter, tol, mm_fun):
    
    if initial_params is None:
        params = np.zeros(n_items)
    else:
        params = initial_params
    converged = NormOfDifferenceTest(tol=tol, order=1)
    for _ in range(max_iter):
        nums, denoms = mm_fun(n_items, data, params)
        params = log_transform((nums + alpha) / (denoms + alpha))
        if converged(params):
            return params
    raise RuntimeError("Did not converge after {} iterations".format(max_iter))
556559
def from_function(self, f, domain=None, N=None):
        
        # rescale f to the unit domain
        domain = self.get_default_domain(domain)
        a,b = domain[0], domain[-1]
        map_ui_ab = lambda t: 0.5*(b-a)*t + 0.5*(a+b)
        args = {'f': lambda t: f(map_ui_ab(t))}
        if N is not None: # N is provided
            nextpow2 = int(np.log2(N))+1
            args['kmin'] = nextpow2
            args['kmax'] = nextpow2+1
            args['raise_no_convergence'] = False
        else:
            args['raise_no_convergence'] = True

        # Find out the right number of coefficients to keep
        coeffs = self.dichotomy(**args)

        return self.from_coeff(coeffs, domain)
461528
def _MStep(self, l, K, x, z_h1, pi_h1, p_h1, p_h, delta_i_j_s):
        
        for k in range(K):
            normconst = 0
            if l == 0: # only need to compute pi at first MM iteration
                pi_h1[k] = np.sum(z_h1.T[k]) / len(z_h1)
            for j in range(self.m):
                omega_k_j = EMMMixPLAggregator.omega(k, j, z_h1, x) # numerator
                denom_sum = 0
                for i in range(self.n):
                    sum1 = 0
                    for t in range(len(x[i])):
                        sum2 = 0
                        sum3 = 0
                        for s in range(t, self.m):
                            sum2 += p_h[k][EMMMixPLAggregator.c(x[i], s)]
                        for s in range(t, self.m + 1):
                            sum3 += delta_i_j_s[i][j][s]
                        sum1 += z_h1[i][k] * (sum2 ** -1) * sum3
                    denom_sum += sum1
                p_h1[k][j] = omega_k_j / denom_sum
                normconst += p_h1[k][j]
            for j in range(self.m):
                p_h1[k][j] /= normconst
711236
def nhill_i(self,x,threshold=0.1,power=2):
        
        x_pow = np.power(x,power)
        threshold_pow = np.power(threshold,power)
        return threshold_pow / (x_pow + threshold_pow) * (1 - x_pow)
742312
def proper(self,x):
        """"""
        x[x < 0] = self.max_fit
        x[np.isnan(x)] = self.max_fit
        x[np.isinf(x)] = self.max_fit
        return x