Helper function for calculation of split-half cross-validation model
        error and signal reliability.
430819
def _xval_model_error(self, fit_spectra, reject_outliers, fit_lb, fit_ub,
                           fitter, func):
        
        set1 = fit_spectra[::2]
        set2 = fit_spectra[1::2]
        errs = []
        signal_select = [] 
        models = []
        signals = []
        for this_set in [set1, set2]:
            choose_transients, model, signal, params, this_idx =\
                self._fit_helper(this_set, reject_outliers,
                                 fit_lb, fit_ub, fitter)
            models.append(np.nanmean(model[choose_transients], 0))
            signals.append(np.nanmean(signal[choose_transients], 0))

            signal_select.append(signal[choose_transients])
                
        #Cross-validation error estimation:
        model_err = np.mean([ut.rmse(models[0], signals[1]),
                              ut.rmse(models[1], signals[0])])
        # Also for the signal:
        signal_err = ut.rmse(np.nanmean(signal_select[0], 0),
                             np.nanmean(signal_select[1], 0))
        # Based on the errors, choose a function. Also report errors:
        return model_err, signal_err
1340767
def np_lst_sq_xval(vecMdl, aryFuncChnk, aryIdxTrn, aryIdxTst):
    

    varNumXval = aryIdxTrn.shape[-1]

    varNumVoxChnk = aryFuncChnk.shape[-1]

    # pre-allocate ary to collect cross-validation
    # error for every xval fold
    aryResXval = np.empty((varNumVoxChnk,
                           varNumXval),
                          dtype=np.float32)

    # loop over cross-validation folds
    for idxXval in range(varNumXval):
        # Get pRF time course models for trn and tst:
        vecMdlTrn = vecMdl[aryIdxTrn[:, idxXval], :]
        vecMdlTst = vecMdl[aryIdxTst[:, idxXval], :]
        # Get functional data for trn and tst:
        aryFuncChnkTrn = aryFuncChnk[
            aryIdxTrn[:, idxXval], :]
        aryFuncChnkTst = aryFuncChnk[
            aryIdxTst[:, idxXval], :]

        # Numpy linalg.lstsq is used to calculate the
        # parameter estimates of the current model:
        vecTmpPe = np.linalg.lstsq(vecMdlTrn,
                                   aryFuncChnkTrn,
                                   rcond=-1)[0]

        # calculate model prediction time course
        aryMdlPrdTc = np.dot(vecMdlTst, vecTmpPe)

        # calculate residual sum of squares between
        # test data and model prediction time course
        aryResXval[:, idxXval] = np.sum(
            (np.subtract(aryFuncChnkTst,
                         aryMdlPrdTc))**2, axis=0)

    return aryResXval
586554
def _xval_choose_funcs(self, fit_spectra, reject_outliers, fit_lb, fit_ub,
                           fitters=[ana.fit_gaussian,ana.fit_two_gaussian],
                           funcs = [ut.gaussian, ut.two_gaussian]):
        

        set1 = fit_spectra[::2]
        set2 = fit_spectra[1::2]

        errs = []
        signal_select = [] 
        # We can loop over functions and try each one out, checking the
        # error in each:
        for fitter in fitters:
            models = []
            signals = []
            for this_set in [set1, set2]:
                choose_transients, model, signal, params, this_idx =\
                    self._fit_helper(this_set, reject_outliers,
                                    fit_lb, fit_ub, fitter)
                models.append(np.nanmean(model[choose_transients], 0))
                signals.append(np.nanmean(signal[choose_transients], 0))

                signal_select.append(signal[choose_transients])
                
            #Cross-validate!
            errs.append(np.mean([ut.rmse(models[0], signals[1]),
                                 ut.rmse(models[1], signals[0])]))
        # We really only need to look at the first two:
        signal_err = ut.rmse(np.nanmean(signal_select[0], 0),
                             np.nanmean(signal_select[1], 0))
        # Based on the errors, choose a function. Also report errors:
        return (fitters[np.argmin(errs)], funcs[np.argmin(errs)], np.min(errs),
                signal_err)
91265
def normal_errors(axes, covariance_matrix, **kwargs):
    
    level = kwargs.pop('level',1)
    traditional_layout = kwargs.pop('traditional_layout',True)
    d = N.diagonal(covariance_matrix)
    ell = ellipse(**kwargs)

    if axes[2,2] < 0:
        axes *= -1

    # Not sure where this factor comes from but it
    # seems to make things work better
    c1 = 2
    axis_lengths = d[:2]
    f = N.linalg.norm(
        ell*axis_lengths,axis=1)

    e0 = -ell.T*d[2]*c1
    e = N.vstack((e0,f))

    _ = dot(e.T,axes).T

    if traditional_layout:
        lon,lat = stereonet_math.cart2sph(_[2],_[0],-_[1])
    else:
        lon,lat = stereonet_math.cart2sph(-_[1],_[0],_[2])
    return list(zip(lon,lat))
225582
function decision2DStumpTrain(data, labels, ix, options) {
    
    options = options || {};
    var numtries = options.numTries || 10;
    
    // choose a dimension at random and pick a best split
    var N = ix.length;
    
    var ri1= 0;
    var ri2= 1;
    if(data[0].length > 2) {
      // more than 2D data. Pick 2 random dimensions
      ri1= randi(0, data[0].length);
      ri2= randi(0, data[0].length);
      while(ri2 == ri1) ri2= randi(0, data[0].length); // must be distinct!
    }
    
    // evaluate class entropy of incoming data
    var H= entropy(labels, ix);
    var bestGain=0; 
    var bestw1, bestw2, bestthr;
    var dots= new Array(ix.length);
    for(var i=0;i<numtries;i++) {
        
        // pick random line parameters
        var alpha= randf(0, 2*Math.PI);
        var w1= Math.cos(alpha);
        var w2= Math.sin(alpha);
        
        // project data on this line and get the dot products
        for(var j=0;j<ix.length;j++) {
          dots[j]= w1*data[ix[j]][ri1] + w2*data[ix[j]][ri2];
        }
        
        // we are in a tricky situation because data dot product distribution
        // can be skewed. So we don't want to select just randomly between
        // min and max. But we also don't want to sort as that is too expensive
        // let's pick two random points and make the threshold be somewhere between them.
        // for skewed datasets, the selected points will with relatively high likelihood
        // be in the high-desnity regions, so the thresholds will make sense
        var ix1= ix[randi(0, N)];
        var ix2= ix[randi(0, N)];
        while(ix2==ix1) ix2= ix[randi(0, N)]; // enforce distinctness of ix2
        var a= Math.random();
        var dotthr= dots[ix1]*a + dots[ix2]*(1-a);
        
        // measure information gain we'd get from split with thr
        var l1=1, r1=1, lm1=1, rm1=1; //counts for Left and label 1, right and label 1, left and minus 1, right and minus 1
        for(var j=0;j<ix.length;j++) {
            if(dots[j] < dotthr) {
              if(labels[ix[j]]==1) l1++;
              else lm1++;
            } else {
              if(labels[ix[j]]==1) r1++;
              else rm1++;
            }
        }
        var t= l1+lm1; 
        l1=l1/t;
        lm1=lm1/t;
        t= r1+rm1;
        r1=r1/t;
        rm1= rm1/t;
        
        var LH= -l1*Math.log(l1) -lm1*Math.log(lm1); // left and right entropy
        var RH= -r1*Math.log(r1) -rm1*Math.log(rm1);
        
        var informationGain= H - LH - RH;
        //console.log("Considering split %f, entropy %f -> %f, %f. Gain %f", thr, H, LH, RH, informationGain);
        if(informationGain > bestGain || i === 0) {
            bestGain= informationGain;
            bestw1= w1;
            bestw2= w2;
            bestthr= dotthr;
        }
    }
    
    model= {};
    model.w1= bestw1;
    model.w2= bestw2;
    model.dotthr= bestthr;
    return model;
}
900222
def decision_function(self, X, **kwargs):
        ""
        if X is None:
            return self._hy_test
        X = self.convert_features(X)
        if len(X) < self.nvar:
            _ = 'Number of variables differ, trained with %s given %s' % (self.nvar, len(X))
            raise RuntimeError(_)
        hist = self._hist
        for node in hist:
            if node.height:
                node.eval(hist)
            else:
                node.eval(X)
        node.normalize()
        r = node.hy
        for i in hist[:-1]:
            i.hy = None
            i.hy_test = None
        gc.collect()
        return r
771556
def adaptive_model_average(X, penalization, method):
    
    n_trials = 100
    print("Adaptive ModelAverage with:")
    print("   estimator: QuicGraphicalLasso (default)")
    print("   n_trials: {}".format(n_trials))
    print("   penalization: {}".format(penalization))
    print("   adaptive-method: {}".format(method))

    # if penalization is random, first find a decent scalar lam_ to build
    # random perturbation matrix around. lam doesn't matter for fully-random.
    lam = 0.5
    if penalization == "random":
        cv_model = QuicGraphicalLassoCV(
            cv=2, n_refinements=6, n_jobs=1, init_method="cov", score_metric=metric
        )
        cv_model.fit(X)
        lam = cv_model.lam_
        print("   lam: {}".format(lam))

    model = AdaptiveGraphicalLasso(
        estimator=ModelAverage(
            n_trials=n_trials, penalization=penalization, lam=lam, n_jobs=1
        ),
        method=method,
    )
    model.fit(X)
    lam_norm_ = np.linalg.norm(model.estimator_.lam_)
    print("   ||lam_||_2: {}".format(lam_norm_))
    return model.estimator_.covariance_, model.estimator_.precision_, lam_norm_
993492
def roc_curve(input:Tensor, targ:Tensor):
    ""
    targ = (targ == 1)
    desc_score_indices = torch.flip(input.argsort(-1), [-1])
    input = input[desc_score_indices]
    targ = targ[desc_score_indices]
    d = input[1:] - input[:-1]
    distinct_value_indices = torch.nonzero(d).transpose(0,1)[0]
    threshold_idxs = torch.cat((distinct_value_indices, LongTensor([len(targ) - 1]).to(targ.device)))
    tps = torch.cumsum(targ * 1, dim=-1)[threshold_idxs]
    fps = (1 + threshold_idxs - tps)
    if tps[0] != 0 or fps[0] != 0:
        fps = torch.cat((LongTensor([0]), fps))
        tps = torch.cat((LongTensor([0]), tps))
    fpr, tpr = fps.float() / fps[-1], tps.float() / tps[-1]
    return fpr, tpr
954488
def _do_analysis_no_cross_validation(self):
        

        # first model is just the mean
        response_term = [Term([LookupFactor(self.y)])]
        model_terms = [Term([])] # empty term is the intercept
        all_model_terms_dict = {x:Term([LookupFactor(x)]) for x in self.list_of_x}
        # ...then add another term for each candidate
        #model_terms += [Term([LookupFactor(c)]) for c in candidates]
        model_desc = ModelDesc(response_term, model_terms)
        self._list_of_fits.append(fm.ols(model_desc, data=self.df).fit())
        # try to improve the model until no improvements can be found

        while all_model_terms_dict:
            # try each x and overwrite the best_fit if we find a better one
            # the first best_fit is the one from the previous round
            ref_fit = self._list_of_fits[-1]
            best_fit = self._list_of_fits[-1]
            best_bic = best_fit.bic
            for x, term in all_model_terms_dict.items():
                # make new_fit, compare with best found so far
                model_desc = ModelDesc(response_term, ref_fit.model.formula.rhs_termlist + [term])
                fit = fm.ols(model_desc, data=self.df).fit()
                if fit.bic < best_bic:
                    best_bic = fit.bic
                    best_fit = fit
                    best_x = x
            # Sometimes, the obtained fit may be better, but contains unsignificant parameters.
            # Correct the fit by removing the unsignificant parameters and estimate again
            best_fit = self._prune(best_fit, p_max=self.p_max)

            # if best_fit does not contain more variables than ref fit, exit
            if len(best_fit.model.formula.rhs_termlist) == len(ref_fit.model.formula.rhs_termlist):
                break
            else:
                self._list_of_fits.append(best_fit)
                all_model_terms_dict.pop(best_x)
        self._fit = self._list_of_fits[-1]
110284
private double computeOutputError(double output, double[] gradient) {

        double error = 0.0;
        double out = outputLayer.output[0];
        double g = output - out;

        error += (0.5*g * g);


        gradient[0] = g;

        return error;
    }
1357104
def crossvalidation_error(errors):
    """"""
    import numpy as np
    return np.mean(errors), np.std(errors)
823922
def scale_errors(cov_axes, confidence_level=0.95):
    
    dof = len(cov_axes)
    x2t = chi2.ppf(confidence_level,dof)
    return N.sqrt(x2t*cov_axes)
262307
def get_err(snr):
    

    quad_terms = np.array(
            [3.11e-3, 1.10e-5, 6.95e-6, 5.05e-6, 4.65e-6, 4.10e-6])
    lin_terms = np.array(
            [-0.869, -2.07e-3, -1.40e-3, -1.03e-3, -1.13e-3, -7.29e-4])
    consts = np.array([104, 0.200, 0.117, 0.114, 0.156, 0.0624])
    err = quad_terms[:,None] * snr**2 + lin_terms[:,None] * snr + consts[:,None]

    # find the minimum of the quadratic function
    min_snr = -lin_terms / (2*quad_terms)
    min_err = quad_terms * min_snr**2 + lin_terms * min_snr + consts
    mask = (snr[:,None] > min_snr).T
    for i in range(0,len(min_err)):
        err[i][mask[i]] = min_err[i]

    return err
786897
public static int computeBestPredictor(int[] data, int dataLen, double[] residualBitsPerSample) {
        int lastError0 = data[-1];
        int lastError1 = data[-1] - data[-2];
        int lastError2 = lastError1 - (data[-2] - data[-3]);
        int lastError3 = lastError2 - (data[-2] - 2 * data[-3] + data[-4]);
        int error, save;
        int totalError0 = 0, totalError1 = 0, totalError2 = 0, totalError3 = 0, totalError4 = 0;
        int i, order;
        
        for (i = 0; i < dataLen; i++) {
            error = data[i];
            totalError0 += Math.abs(error);
            save = error;
            error -= lastError0;
            totalError1 += Math.abs(error);
            lastError0 = save;
            save = error;
            error -= lastError1;
            totalError2 += Math.abs(error);
            lastError1 = save;
            save = error;
            error -= lastError2;
            totalError3 += Math.abs(error);
            lastError2 = save;
            save = error;
            error -= lastError3;
            totalError4 += Math.abs(error);
            lastError3 = save;
        }
        
        if (totalError0 < Math.min(Math.min(Math.min(totalError1, totalError2), totalError3), totalError4))
            order = 0;
        else if (totalError1 < Math.min(Math.min(totalError2, totalError3), totalError4))
            order = 1;
        else if (totalError2 < Math.min(totalError3, totalError4))
            order = 2;
        else if (totalError3 < totalError4)
            order = 3;
        else
            order = 4;
        
        // Estimate the expected number of bits per residual signal sample.
        // 'total_error*' is linearly related to the variance of the residual
        // signal, so we use it directly to compute E(|x|)
        residualBitsPerSample[0] = (double) ((totalError0 > 0) ? Math.log(M_LN2 * (double) totalError0 / (double) dataLen) / M_LN2 : 0.0);
        residualBitsPerSample[1] = (double) ((totalError1 > 0) ? Math.log(M_LN2 * (double) totalError1 / (double) dataLen) / M_LN2 : 0.0);
        residualBitsPerSample[2] = (double) ((totalError2 > 0) ? Math.log(M_LN2 * (double) totalError2 / (double) dataLen) / M_LN2 : 0.0);
        residualBitsPerSample[3] = (double) ((totalError3 > 0) ? Math.log(M_LN2 * (double) totalError3 / (double) dataLen) / M_LN2 : 0.0);
        residualBitsPerSample[4] = (double) ((totalError4 > 0) ? Math.log(M_LN2 * (double) totalError4 / (double) dataLen) / M_LN2 : 0.0);
        
        return order;
    }
1183373
def error_rate(predictions, labels):
  """"""
  return 100.0 - (
      100.0 *
      np.sum(np.argmax(predictions, 1) == np.argmax(labels, 1)) /
      predictions.shape[0])
1307370
public void crossValidateSetSigma(GeneralDataset<L, F> dataset,int kfold, final Scorer<L> scorer, LineSearcher minimizer) {
    System.err.println("##in Cross Validate, folds = " + kfold);
    System.err.println("##Scorer is " + scorer);

    featureIndex = dataset.featureIndex;
    labelIndex = dataset.labelIndex;

    final CrossValidator<L, F> crossValidator = new CrossValidator<L, F>(dataset,kfold);
    final Function<Triple<GeneralDataset<L, F>,GeneralDataset<L, F>,CrossValidator.SavedState>,Double> score =
      new Function<Triple<GeneralDataset<L, F>,GeneralDataset<L, F>,CrossValidator.SavedState>,Double> ()
      {
        public Double apply (Triple<GeneralDataset<L, F>,GeneralDataset<L, F>,CrossValidator.SavedState> fold) {
          GeneralDataset<L, F> trainSet = fold.first();
          GeneralDataset<L, F> devSet   = fold.second();

          double[] weights = (double[])fold.third().state;
          double[][] weights2D;

          weights2D = trainWeights(trainSet, weights,true); // must of course bypass sigma tuning here.

          fold.third().state = ArrayUtils.flatten(weights2D);

          LinearClassifier<L, F> classifier = new LinearClassifier<L, F>(weights2D, trainSet.featureIndex, trainSet.labelIndex);

          double score = scorer.score(classifier, devSet);
          //System.out.println("score: "+score);
          System.out.print(".");
          return score;
        }
      };

    Function<Double,Double> negativeScorer =
      new Function<Double,Double> ()
      {
        public Double apply(Double sigmaToTry) {
          //sigma = sigmaToTry;
          setSigma(sigmaToTry);
          Double averageScore = crossValidator.computeAverage(score);
          System.err.print("##sigma = "+getSigma()+" ");
          System.err.println("-> average Score: "+averageScore);
          return -averageScore;
        }
      };

    double bestSigma = minimizer.minimize(negativeScorer);
    System.err.println("##best sigma: " + bestSigma);
    setSigma(bestSigma);
  }
265211
def cross_validate(data=None, folds=5, repeat=1, metrics=None,
                   reporters=None, model_def=None, **kwargs):
    
    md_kwargs = {}
    if model_def is None:
        for arg in ModelDefinition.params:
            if arg in kwargs:
                md_kwargs[arg] = kwargs.pop(arg)
        model_def = ModelDefinition(**md_kwargs)
    if metrics is None:
        metrics = []
    if reporters is None:
        reporters = []
    metrics = [MetricReporter(metric) for metric in metrics]
    results = modeling.cross_validate(model_def, data, folds, repeat=repeat, **kwargs)
    for r in reporters + metrics:
        r.process_results(results)
    return CVResult(results, reporters, metrics)