Top 3 alternatives 20 empirical moment conditions
976644
def calcMomentsMatlabEmpirical(params):
    """"""
    alpha = params[0]
    a = params[1:5]
    b = params[5:]
    p1 = alpha*a+(1-alpha)*b
    p21 = alpha*a[0]*a[1:]/(1-a[0])+(1-alpha)*b[0]*b[1:]/(1-b[0])
    p22 = alpha*a[1]*np.hstack((a[0],a[2:]))/(1-a[1])+(1-alpha)*b[1]*np.hstack((b[0],b[2:]))/(1-b[1])
    p23 = alpha*a[2]*np.hstack((a[:2],a[3]))/(1-a[2])+(1-alpha)*b[2]*np.hstack((b[:2],b[3]))/(1-b[2])
    p24 = alpha*a[3]*a[:3]/(1-a[3])+(1-alpha)*b[3]*b[:3]/(1-b[3])
    p3 = np.array([
        alpha*a[0]*a[2]*a[3]/(1-a[2])/(a[0]+a[1])+(1-alpha)*b[0]*b[2]*b[3]/(1-b[2])/(b[0]+b[1]),
        alpha*a[0]*a[1]*a[3]/(1-a[3])/(a[1]+a[2])+(1-alpha)*b[0]*b[1]*b[3]/(1-b[3])/(b[1]+b[2]),
        alpha*a[0]*a[1]*a[2]/(1-a[0])/(a[3]+a[2])+(1-alpha)*b[0]*b[1]*b[2]/(1-b[0])/(b[3]+b[2]),
        alpha*a[2]*a[1]*a[3]/(1-a[1])/(a[0]+a[3])+(1-alpha)*b[2]*b[1]*b[3]/(1-b[1])/(b[0]+b[3])
        ])
    return np.concatenate((p1,p21,p22,p23,p24,p3))
168700
def createExperimentArgs():
  """"""
  experimentArguments = []
  # for n in [300, 500, 700, 900, 1100, 1300, 1500, 1700, 1900, 2100, 2300,
  #           2500, 2700, 2900, 3100, 3300, 3500, 3700, 3900]:
  for n in [1500, 1700, 1900, 2100]:
    for a in [128]:
      # Some parameter combinations are just not worth running!
      if ( a==64 and n<=1500 ) or ( a==128 and n<= 1900 ) or ( a==256 ):
        experimentArguments.append(
          ("./sdr_calculations2", "results_errorbars/temp_"+str(n)+"_"+str(a)+".csv",
           "200000", str(n), str(a), "0"),
        )
  return experimentArguments
862727
def moment(expr, order, central=False):
    
    if not isinstance(order, six.integer_types):
        raise ValueError('Only integer-ordered moments are supported.')
    if order < 0:
        raise ValueError('Only non-negative orders are supported.')
    output_type = _stats_type(expr)
    return _reduction(expr, Moment, output_type, _order=order, _center=central)
901307
def d_distributionparams(self):
        """"""
        if not self._d_distributionparams:
            dx = 1.0e-3
            def f_alpha(alpha):
                return DiscreteGamma(alpha, self.beta_lambda, self.ncats)
            def f_beta(beta):
                return DiscreteGamma(self.alpha_lambda, beta, self.ncats)
            assert set(self.distributionparams) == {'alpha_lambda', 'beta_lambda'}
            for (param, f) in [('alpha_lambda', f_alpha), ('beta_lambda', f_beta)]:
                pvalue = getattr(self, param)
                dparam = scipy.misc.derivative(f, pvalue, dx, n=1, order=5)
                assert dparam.shape == (self.ncats,)
                for stepchange in [0.5, 2]: # make sure robust to step size
                    dparam2 = scipy.misc.derivative(f, pvalue, stepchange * dx,
                            n=1, order=5)
                    assert scipy.allclose(dparam, dparam2, atol=1e-5, rtol=1e-4), (
                            "Numerical derivative of {0} at {1} "
                            "differs for step {2} and {3}: {4} and {5}"
                            ", respectively.").format(param, pvalue,
                            dx, dx * stepchange, dparam, dparam2)
                self._d_distributionparams[param] = dparam
        return self._d_distributionparams
394547
private static function hsl2rgb($h, $s, $l, &$r, &$g, &$b)
    {
        if ($s == 0) {
            $r = round($l * 255);
            $g = $r;
            $b = $r;

            return;
        }

        $temp1;
        if ($l < 0.5) {
            $temp1 = $l * (1.0 + $s);
        } else {
            $temp1 = $l + $s - $l * $s;
        }

        $temp2 = 2 * $l - $temp1;

        $h = $h / 360;

        $tempR = $h + (1 / 3);
        $tempG = $h;
        $tempB = $h - (1 / 3);

        if ($tempR < 0) {
            $tempR++;
        } elseif ($tempR > 1) {
            $tempR--;
        }

        if ($tempG < 0) {
            $tempG++;
        } elseif ($tempG > 1) {
            $tempG--;
        }

        if ($tempB < 0) {
            $tempB++;
        } elseif ($tempB > 1) {
            $tempB--;
        }

        if (6 * $tempR < 1) {
            $r = $temp2 + ($temp1 - $temp2) * 6 * $tempR;
        } elseif (2 * $tempR < 1) {
            $r = $temp1;
        } elseif (3 * $tempR < 2) {
            $r = $temp2 + ($temp1 - $temp2) * ((2 / 3) - $tempR) * 6;
        } else {
            $r = $temp2;
        }

        if (6 * $tempG < 1) {
            $g = $temp2 + ($temp1 - $temp2) * 6 * $tempG;
        } elseif (2 * $tempG < 1) {
            $g = $temp1;
        } elseif (3 * $tempG < 2) {
            $g = $temp2 + ($temp1 - $temp2) * ((2 / 3) - $tempG) * 6;
        } else {
            $g = $temp2;
        }

        if (6 * $tempB < 1) {
            $b = $temp2 + ($temp1 - $temp2) * 6 * $tempB;
        } elseif (2 * $tempB < 1) {
            $b = $temp1;
        } elseif (3 * $tempB < 2) {
            $b = $temp2 + ($temp1 - $temp2) * ((2 / 3) - $tempB) * 6;
        } else {
            $b = $temp2;
        }

        $r = round($r * 255);
        $g = round($g * 255);
        $b = round($b * 255);
    }
986940
def _parameter_constraints(self, theta_E, gamma, q, phi_G, s_scale):
        
        if theta_E < 0:
            theta_E = 0
        if s_scale < 0.00000001:
            s_scale = 0.00000001
        if gamma < 1.2:
            gamma = 1.2
            theta_E = 0
        if gamma > 2.9:
            gamma = 2.9
            theta_E = 0
        if q < 0.01:
            q = 0.01
            theta_E = 0
        if q > 1:
            q = 1.
            theta_E = 0
        return theta_E, gamma, q, phi_G, s_scale
270140
def main():
    
    dataset = 'RD'

    observables = ['log(<\sigma v>)', '\Omega_{\chi}h^2', 'log(\sigma_p^{SI})']
    var  = ['log(m_{\chi})']
    var += ['log(C_1)', 'log(C_2)', 'log(C_3)', 'log(C_4)', 'log(C_5)', 'log(C_6)']
    var += observables

    plot_vs_mass(dataset, observables, 'mass_vs_observables.png')
    plot_oneD(dataset, var, 'oneD.png')
    pairplot(dataset, var, 'pairplot.png')
84413
def moments(self):
        """"""
        timepoints = self.ma.delays
        response = self.response
        moment1 = statstools.calc_mean_time(timepoints, response)
        moment2 = statstools.calc_mean_time_deviation(
            timepoints, response, moment1)
        return numpy.array([moment1, moment2])
1251521
def findbestparams_iso(xsamples):
    
    p = 0.5*np.sum(xsamples.T[3:]**2,axis=0)
    r = np.sum(xsamples.T[:3]**2,axis=0)
    return np.abs(leastsq(deltaH_iso,np.array([10.,10.]), Dfun = None , col_deriv=1,args=(p,r,))[0])
1269257
def _compute_deviations(self, sample, n):
        

        try:
            coeffs = fit_first_and_second_harmonics(self.sample.values[0],
                                                    self.sample.values[2])
            coeffs = coeffs[0]
            model = first_and_second_harmonic_function(self.sample.values[0],
                                                       coeffs)
            residual = self.sample.values[2] - model

            c = fit_upper_harmonic(residual, sample.values[2], n)
            covariance = c[1]
            ce = np.diagonal(covariance)
            c = c[0]

            a = c[1] / self.sma / sample.gradient
            b = c[2] / self.sma / sample.gradient

            # this comes from the old code. Likely it was based on
            # empirical experience with the STSDAS task, so we leave
            # it here without too much thought.
            gre = self.grad_r_error if self.grad_r_error is not None else 0.64

            a_err = abs(a) * np.sqrt((ce[1] / c[1])**2 + gre**2)
            b_err = abs(b) * np.sqrt((ce[2] / c[2])**2 + gre**2)

        except Exception:    # we want to catch everything
            a = b = a_err = b_err = None

        return a, b, a_err, b_err
229964
def rbf(ra, coeff, mag):
    
    a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = gc(coeff, mag)
    term1 = a1ca + a1cb * mag + a1cc * np.log(ra + a1cd*np.exp(a1ce*mag))
    term2 = a2ca + a2cb * mag
    term3 = a2cd*np.exp(a2ce*mag)
    return np.exp((term1 - term2) / a2cc) - term3
1245739
def params_of_mean(value=array([-.005, 1.]), tau=.1, rate=4.):
    

    def logp(value, tau, rate):
        if value[1] > 0 and value[1] + value[0] * 110 > 0:
            return normal_like(value[0], 0., tau) + \
                exponential_like(value[1], rate)
        else:
            return -Inf

    def random(tau, rate):
        val = zeros(2)
        val[0] = rnormal(0., tau)
        val[1] = rexponential(rate)
        while val[1] < 0 or val[1] + val[0] * 110 <= 0:
            val[0] = rnormal(0., tau)
            val[1] = rexponential(rate)
        return val
143923
def _generate_moment_matrix(self, n_vars, block_index, processed_entries,
                                monomialsA, monomialsB, ppt=False):
        
        row_offset = 0
        if block_index > 0:
            for block_size in self.block_struct[0:block_index]:
                row_offset += block_size ** 2
        N = len(monomialsA)*len(monomialsB)
        func = partial(assemble_monomial_and_do_substitutions,
                       monomialsA=monomialsA, monomialsB=monomialsB, ppt=ppt,
                       substitutions=self.substitutions,
                       pure_substitution_rules=self.pure_substitution_rules)
        if self._parallel:
            pool = Pool()
            # This is just a guess and can be optimized
            chunksize = int(max(int(np.sqrt(len(monomialsA) * len(monomialsB) *
                                        len(monomialsA) / 2) / cpu_count()),
                                1))
        for rowA in range(len(monomialsA)):
            if self._parallel:
                iter_ = pool.map(func, [(rowA, columnA, rowB, columnB)
                                        for rowB in range(len(monomialsB))
                                        for columnA in range(rowA,
                                                             len(monomialsA))
                                        for columnB in range((rowA == columnA)*rowB,
                                                             len(monomialsB))],
                                 chunksize)
                print_criterion = processed_entries + len(iter_)
            else:
                iter_ = imap(func, [(rowA, columnA, rowB, columnB)
                                    for columnA in range(rowA, len(monomialsA))
                                    for rowB in range(len(monomialsB))
                                    for columnB in range((rowA == columnA)*rowB,
                                                         len(monomialsB))])
            for columnA, rowB, columnB, monomial in iter_:
                processed_entries += 1
                n_vars = self._push_monomial(monomial, n_vars,
                                             row_offset, rowA,
                                             columnA, N, rowB,
                                             columnB, len(monomialsB),
                                             prevent_substitutions=True)
                if self.verbose > 0 and (not self._parallel or
                                         processed_entries == self.n_vars or
                                         processed_entries == print_criterion):
                    percentage = processed_entries / self.n_vars
                    time_used = time.time()-self._time0
                    eta = (1.0 / percentage) * time_used - time_used
                    hours = int(eta/3600)
                    minutes = int((eta-3600*hours)/60)
                    seconds = eta-3600*hours-minutes*60

                    msg = ""
                    if self.verbose > 1 and self._parallel:
                        msg = ", working on block {:0} with {:0} processes with a chunksize of {:0d}"\
                              .format(block_index, cpu_count(),
                                      chunksize)
                    msg = "{:0} (done: {:.2%}, ETA {:02d}:{:02d}:{:03.1f}"\
                          .format(n_vars, percentage, hours, minutes, seconds) + \
                          msg
                    msg = "\r\x1b[KCurrent number of SDP variables: " + msg + ")"
                    sys.stdout.write(msg)
                    sys.stdout.flush()

        if self._parallel:
            pool.close()
            pool.join()
        if self.verbose > 0:
            sys.stdout.write("\r")
        return n_vars, block_index + 1, processed_entries
843012
function generateVariationsOverConditions(conditionList) {
      const curCondition = conditionList[0];

      if (!curCondition) {
        return [];
      }

      // get combinations from the n - 1th dimension
      const nextVariations = generateVariationsOverConditions(
        conditionList.slice(1)
      );
      const variations = [];

      if (!nextVariations.length) {
        for (const curConditionValue of conditionalEnv[curCondition]) {
          const variationConditions = {};
          variationConditions[curCondition] = curConditionValue;
          variations.push(variationConditions);
        }
      }

      // multiply the combinations of the n - 1 dimension by the cominations of this nth dimension
      for (const nextVariation of nextVariations) {
        for (const curConditionValue of conditionalEnv[curCondition]) {
          const variationConditions = Object.assign({}, nextVariation);
          variationConditions[curCondition] = curConditionValue;
          variations.push(variationConditions);
        }
      }

      return variations;
    }
398734
def moments(I_xy_input, x, y):
    
    I_xy = copy.deepcopy(I_xy_input)
    background = np.minimum(0, np.min(I_xy))
    I_xy -= background
    x_ = np.sum(I_xy * x)
    y_ = np.sum(I_xy * y)
    r = (np.max(x) - np.min(x)) / 3.
    mask = mask_util.mask_sphere(x, y, center_x=x_, center_y=y_, r=r)
    Q_xx = np.sum(I_xy * mask * (x - x_) ** 2)
    Q_xy = np.sum(I_xy * mask * (x - x_) * (y - y_))
    Q_yy = np.sum(I_xy * mask * (y - y_) ** 2)
    return Q_xx, Q_xy, Q_yy, background / np.mean(I_xy)
459173
def _get_a2(bbar, dbar, slip_moment, mmax):
        
        return ((dbar - bbar) / bbar) * (slip_moment / _scale_moment(mmax))
1074973
def moments(data,cntr):
    
    total = data.sum()
    #X, Y = np.indices(data.shape)
    #x = (X*data).sum()/total
    #y = (Y*data).sum()/total
    x,y = cntr
    xi = int(x)
    yi = int(y)
    if xi < 0 or xi >= data.shape[1] or yi < 0 or yi >= data.shape[0]:
        raise ValueError
    col = data[:, xi]
    width_x = np.sqrt(abs(((np.arange(col.size)-y)**2*col).sum()/col.sum()))
    row = data[yi, :]
    width_y = np.sqrt(abs(((np.arange(row.size)-x)**2*row).sum()/row.sum()))
    height = data.max()
    return height, x, y, width_x, width_y