correlation coefficient
        input:
        s: simulated
        o: observed
        output:
        correlation: correlation coefficient
973720
def correlation(s, o):
    
    # s,o = filter_nan(s,o)
    if s.size == 0:
        corr = np.NaN
    else:
        corr = np.corrcoef(o, s)[0, 1]
    return corr
1172483
def analyzeParameters(expName, suite):
  
  print("\n================",expName,"=====================")
  try:
    expParams = suite.get_params(expName)
    pprint.pprint(expParams)

    for p in ["boost_strength", "k", "learning_rate", "weight_sparsity",
              "k_inference_factor", "boost_strength_factor",
              "c1_out_channels", "c1_k", "learning_rate_factor",
              "batches_in_epoch",
              ]:
      if p in expParams and type(expParams[p]) == list:
        print("\n",p)
        for v1 in expParams[p]:
          # Retrieve the last totalCorrect from each experiment
          # Print them sorted from best to worst
          values, params = suite.get_values_fix_params(
            expName, 0, "testerror", "last", **{p:v1})
          v = np.array(values)
          try:
            print("Average/min/max for", p, v1, "=", v.mean(), v.min(), v.max())
            # sortedIndices = v.argsort()
            # for i in sortedIndices[::-1]:
            #   print(v[i],params[i]["name"])
          except:
            print("Can't compute stats for",p)

  except:
    print("Couldn't load experiment",expName)
91461
def paired(x, y):
    
    samples = ''
    while samples not in ['i', 'r', 'I', 'R', 'c', 'C']:
        print('\nIndependent or related samples, or correlation (i,r,c): ',)
        samples = raw_input()

    if samples in ['i', 'I', 'r', 'R']:
        print('\nComparing variances ...',)
        # USE O'BRIEN'S TEST FOR HOMOGENEITY OF VARIANCE, Maxwell & delaney, p.112
        r = obrientransform(x, y)
        f, p = F_oneway(pstat.colex(r, 0), pstat.colex(r, 1))
        if p < 0.05:
            vartype = 'unequal, p=' + str(round(p, 4))
        else:
            vartype = 'equal'
        print(vartype)
        if samples in ['i', 'I']:
            if vartype[0] == 'e':
                t, p = ttest_ind(x, y, 0)
                print('\nIndependent samples t-test:  ', round(t, 4), round(p, 4))
            else:
                if len(x) > 20 or len(y) > 20:
                    z, p = ranksums(x, y)
                    print('\nRank Sums test (NONparametric, n>20):  ', round(z, 4), round(p, 4))
                else:
                    u, p = mannwhitneyu(x, y)
                    print('\nMann-Whitney U-test (NONparametric, ns<20):  ', round(u, 4), round(p, 4))
        else:  # RELATED SAMPLES
            if vartype[0] == 'e':
                t, p = ttest_rel(x, y, 0)
                print('\nRelated samples t-test:  ', round(t, 4), round(p, 4))
            else:
                t, p = ranksums(x, y)
                print('\nWilcoxon T-test (NONparametric):  ', round(t, 4), round(p, 4))
    else:  # CORRELATION ANALYSIS
        corrtype = ''
        while corrtype not in ['c', 'C', 'r', 'R', 'd', 'D']:
            print('\nIs the data Continuous, Ranked, or Dichotomous (c,r,d): ',)
            corrtype = raw_input()
        if corrtype in ['c', 'C']:
            m, b, r, p, see = linregress(x, y)
            print('\nLinear regression for continuous variables ...')
            lol = [['Slope', 'Intercept', 'r', 'Prob', 'SEestimate'], [round(m, 4), round(b, 4), round(r, 4), round(p, 4), round(see, 4)]]
            pstat.printcc(lol)
        elif corrtype in ['r', 'R']:
            r, p = spearmanr(x, y)
            print('\nCorrelation for ranked variables ...')
            print("Spearman's r: ", round(r, 4), round(p, 4))
        else: # DICHOTOMOUS
            r, p = pointbiserialr(x, y)
            print('\nAssuming x contains a dichotomous variable ...')
            print('Point Biserial r: ', round(r, 4), round(p, 4))
    print('\n\n')
    return None
976545
def coefficients(self):
        """"""
        out = numpy.array([self.A[key] for key in self.keys])
        out = numpy.rollaxis(out, -1)
        return out
383521
def obfn_dfd(self):
        
        """
        XF = sl.rfftn(self.obfn_fvar(), mp_Nv, mp_axisN)
        DX = np.moveaxis(sl.irfftn(sl.inner(mp_Df, XF, mp_axisM),
                                   mp_Nv, mp_axisN), mp_axisM,
                         self.cri.axisM)
        return np.sum((self.W*(DX-self.S))**2)/2.0
497411
def _u_distance_correlation_sqr_naive(x, y, exponent=1):
    """"""
    return _distance_sqr_stats_naive_generic(
        x, y,
        matrix_centered=_u_distance_matrix,
        product=u_product,
        exponent=exponent).correlation_xy
576138
def resample_zeroinflation_variables(self):
        
        # TODO: move this to cython?
        T, N, C, D, b = self.T, self.D_emission, self.C, self.D, self.emission_distn.b
        indptr = [0]
        indices = []
        vals = []
        offset = 0
        X = np.hstack((self.gaussian_states, self.inputs))
        for t in range(T):
            # Evaluate probability of data
            y_t = np.zeros(N)
            ns_t = self.data.indices[self.data.indptr[t]:self.data.indptr[t+1]]
            y_t[ns_t] = self.data.data[self.data.indptr[t]:self.data.indptr[t+1]]
            ll = self.emission_distn._elementwise_log_likelihood((X[t], y_t))
            ll = ll.ravel()

            # Evaluate the probability that each emission was "exposed",
            # i.e. p(z_tn = 1 | y_tn, x_tn)
            log_p_exposed = np.log(self.rho) + ll
            log_p_exposed -= np.log(np.exp(log_p_exposed) + (1-self.rho) * (y_t == 0))

            # Sample zero inflation mask
            z_t = np.random.rand(N) < np.exp(log_p_exposed)

            # Construct the sparse matrix
            t_inds = np.where(z_t)[0]
            indices.append(t_inds)
            vals.append(y_t[t_inds])
            offset += t_inds.size
            indptr.append(offset)

        # Construct a sparse matrix
        vals = np.concatenate(vals)
        indices = np.concatenate(indices)
        indptr = np.array(indptr)
        self.masked_data = csr_matrix((vals, indices, indptr), shape=(T, N))
1272647
def infer(args):
    

    expected_output_files = ["inferred.pkl"]
    if args.plotting:
        expected_output_files.extend([each.format(args.plot_format) \
            for each in "chain.{}", "corner.{}", "acceptance-fractions.{}",
            "autocorrelation.{}"])

    # Optimise them first.
    model, data, metadata, optimised_theta = optimise(args,
        expected_output_files=expected_output_files, __return_result=True)

    # Get the inference parameters from the model configuration.
    kwargs = model._configuration.get("infer", {})
    [kwargs.pop(k, None) \
        for k in ("debug", "full_output", "initial_proposal", "data")]

    try:
        theta, chains, lnprobability, acceptance_fractions, sampler, info_dict \
            = model.infer(data, initial_proposal=optimised_theta, 
                full_output=True, debug=args.debug, 
                __keep_convolution_functions=True,
                __show_progress_bar=True, **kwargs)

    except:
        logger.exception("Failed to infer model parameters")
        raise

    metadata["inferred"] = {
        "theta": theta,
        "chi_sq": info_dict["chi_sq"],
        "dof": info_dict["dof"],
        "r_chi_sq": info_dict["chi_sq"]/info_dict["dof"]
    }

    logger.info("Inferred parameters are:")
    _announce_theta(theta)
    
    # Write the results to file.
    _write_output(_prefix(args, "inferred.pkl"), metadata)

    # Write the chains, etc to disk.
    if args.save_chain_files:
        filename = _prefix(args, "chains.pkl")
        with open(filename, "wb") as f:
            pickle.dump(
                (chains, lnprobability, acceptance_fractions, info_dict), f, -1)
        logger.info("Saved chains to {}".format(filename))

    # Make plots.
    if args.plotting:
        burn = info_dict["burn"]

        # Any truth values to plot?
        truths = model._configuration.get("truths", None)
        if truths:
            truths = [truths.get(p, np.nan) for p in model.parameters]

        # Labels?
        labels = model._configuration.get("labels", {})
        labels = [labels.get(p, p) for p in info_dict["parameters"]]
        
        # Acceptance fractions.
        fig = sick.plot.acceptance_fractions(acceptance_fractions,
            burn=burn)
        _ = _prefix(args, "acceptance-fractions.{}".format(args.plot_format))
        fig.savefig(_)
        logger.info("Saved acceptance fractions figure to {}".format(_))

        # Autocorrelation.
        fig = sick.plot.normalised_autocorrelation_function(chains, burn=burn)
        _ = _prefix(args, "auto-correlation.{}".format(args.plot_format))
        fig.savefig(_)
        logger.info("Saved auto-correlation figure to {}".format(_))

        # Chains.
        fig = sick.plot.chains(chains, labels=labels, burn=burn,
            truths=truths)
        _ = _prefix(args, "chains.{}".format(args.plot_format))
        fig.savefig(_)
        logger.info("Saved chains figure to {}".format(_))

        # Corner plots (astrophysical + all).
        N = len(model.grid_points.dtype.names)
        fig = sick.plot.corner(chains[:, burn:, :N].reshape(-1, N),
            labels=labels, truths=truths[:N] if truths else None)
        _ = _prefix(args, "corner.{}".format(args.plot_format))
        fig.savefig(_)
        logger.info("Saved corner plot (astrophysical parameters) to {}"\
            .format(_))

        if len(model.parameters) > N:
            fig = sick.plot.corner(chains[:, burn:, :].reshape(-1, len(theta)),
                labels=labels, truths=truths)
            _ = _prefix(args, "corner-all.{}".format(args.plot_format))
            fig.savefig(_)
            logger.info("Saved corner plot (all parameters) to {}".format(_))

        # Projections.
        # Note here we need to scale the chains back to redshift so the data
        # are generated properly.
        fig = sick.plot.projection(data, model, 
            chains=chains[:, burn:, :]/info_dict["scales"],
            parameters=theta.keys())
        _ = _prefix(args, "projection.{}".format(args.plot_format))
        fig.savefig(_)
        logger.info("Saved projection plot to {}".format(_))

    model._destroy_convolution_functions()

    return None
79668
def column_correlations(self, X):
        """"""
        utils.validation.check_is_fitted(self, 's_')

        X_global = self._build_X_global(X)
        row_pc = self._row_coordinates_from_global(X_global)

        return pd.DataFrame({
            component: {
                feature: row_pc[component].corr(X_global[feature].to_dense())
                for feature in X_global.columns
            }
            for component in row_pc.columns
        })
1082084
def correlation_model(prediction, fm):
    
    (_, r_x) = calc_resize_factor(prediction, fm.image_size)
    fdm = compute_fdm(fm, scale_factor = r_x)
    return np.corrcoef(fdm.flatten(), prediction.flatten())[0,1]
772271
def cross_correlate(observed, template_dispersion, template_fluxes,
    rebin="template", wavelength_range=None, continuum_degree=-1,
    apodize=0.10, rescale=False, z_limits=None, full_output=False):
    

    template_fluxes = np.atleast_2d(template_fluxes)
    template_dispersion = np.array(template_dispersion)
    if template_dispersion.shape[0] != template_fluxes.shape[1]:
        raise ValueError("template dispersion must have size (N_pixels,) and "\
            "template fluxes must have size (N_models, N_pixels)")

    try:
        continuum_degree = int(continuum_degree)
    except (TypeError, ValueError):
        raise TypeError("continuum order must be an integer-like object")

    if not (1 > apodize >= 0):
        raise ValueError("apodize fraction must be between 0 and 1")

    if rebin.lower() == "template":
        # Put the template fluxes onto the observed dispersion map.
        dispersion = observed.disp
        template_flux = template_fluxes \
            * resample(template_dispersion, observed.disp)
        observed_flux = observed.flux.copy()
        observed_ivar = observed.ivariance.copy()

    elif rebin.lower() == "observed":
        # Put the observed fluxes onto the template dispersion map.
        dispersion = template_dispersion
        template_flux = template_fluxes
        mat = resample(observed.disp, template_dispersion)
        observed_flux = observed.flux.copy() * mat
        observed_ivar = observed.ivariance.copy() * mat
        
    else:
        raise ValueError("rebin must be either `template` or `observed`")

    if wavelength_range is not None:
        if not isinstance(wavelength_range, (tuple, list, np.ndarray)) \
        or len(wavelength_range) != 2:
            raise TypeError("wavelength range must be a two length tuple")

        indices = np.clip(
            template_dispersion.searchsorted(wavelength_range) + [0, 1],
            0, template_dispersion.size)
        template_dispersion = template_dispersion[indices[0]:indices[1]]
        template_flux = template_flux[indices[0]:indices[1]]
        observed_flux = observed_flux[indices[0]:indices[1]]

    N_templates = template_fluxes.shape[0]

    # Clip out the non-finite edges, if they exist.
    _ = np.where(np.isfinite(template_flux[0] * observed_flux))[0]
    l_idx, u_idx = _[0], _[-1]

    dispersion = dispersion[l_idx:u_idx]
    observed_flux = observed_flux[l_idx:u_idx]
    observed_ivar = observed_ivar[l_idx:u_idx]
    template_flux = template_flux[:, l_idx:u_idx]

    # Ensure an even number of points.
    N = u_idx - l_idx
    N = N - 1 if N % 2 > 0 else N
    dispersion = dispersion[:N]
    observed_flux = observed_flux[:N]
    observed_ivar = observed_ivar[:N]
    template_flux = template_flux[:, :N]

    finite = np.isfinite(observed_flux * observed_ivar) * (observed_flux > 1e-3)
    
    # Continuum.
    if continuum_degree >= 0:
        coefficients = np.polyfit(dispersion[finite], observed_flux[finite],
            continuum_degree)
        observed_flux /= np.polyval(coefficients, dispersion)

    # Interpolate over non-finite pixels.
    observed_flux[~finite] = np.interp(dispersion[~finite], dispersion[finite],
        observed_flux[finite], left=1, right=1)

    # Scale the flux level to that the template intensities
    if rescale:
        observed_flux = (observed_flux * template_flux.ptp()) \
            + template_flux.min()

    # Apodize edges.
    edge_buffer = apodize * (dispersion[-1] - dispersion[0])
    low_w_indices = np.nonzero(dispersion < dispersion[0] + edge_buffer)[0]
    high_w_indices = np.nonzero(dispersion > dispersion[-1] - edge_buffer)[0]

    apod_curve = np.ones(N, dtype='d')
    apod_curve[low_w_indices] = (1.0 + np.cos(np.pi*(
        1.0 - (dispersion[low_w_indices] - dispersion[0])/edge_buffer)))/2.
    apod_curve[high_w_indices] = (1.0 + np.cos(np.pi*(
        1.0 - (dispersion[-1] - dispersion[high_w_indices])/edge_buffer)))/2.

    apod_observed_flux = observed_flux * apod_curve
    apod_template_flux = template_flux * apod_curve

    fft_observed_flux = np.fft.fft(apod_observed_flux)
    fft_template_flux = np.fft.fft(apod_template_flux)
    template_flux_corr = (fft_observed_flux * fft_template_flux.conjugate())
    template_flux_corr /= np.sqrt(
        np.inner(apod_observed_flux, apod_observed_flux))

    z_array = np.array(dispersion.copy())/dispersion[N/2] - 1.0

    z = np.ones(N_templates) * np.nan
    z_err = np.ones(N_templates) * np.nan
    R = np.ones(N_templates) * np.nan
    for i in range(N_templates):

        denominator = np.sqrt(
            np.inner(apod_template_flux[i, :], apod_template_flux[i, :]))
        flux_correlation = template_flux_corr[i, :] / denominator 
        correlation = np.fft.ifft(flux_correlation).real

        # Reflect about zero
        ccf = np.zeros(N)
        ccf[:N/2] = correlation[N/2:]
        ccf[N/2:] = correlation[:N/2]

        # Get height and redshift of best peak
        h = ccf.max()

        # Scale the CCF
        #ccf -= ccf.min()
        #ccf *= (h/ccf.max())
        
        R[i] = h
        z[i] = z_array[ccf.argmax()]
        try:
            z_err[i] = (np.ptp(z_array[np.where(ccf >= 0.5*h)])/2.35482)**2
        except ValueError:
            continue

    # Limits:
    if z_limits is not None:
        R[~((z_limits[1] >= z) * (z >= z_limits[0]))] = np.nan

    # Re-measure the velocity at the best peak.
    index = np.nanargmax(R)

    denominator = np.sqrt(
        np.inner(apod_template_flux[index, :], apod_template_flux[index, :]))

    flux_correlation = template_flux_corr[index, :] / denominator
    correlation = np.fft.ifft(flux_correlation).real

    # FReflect about zero.
    ccf = np.zeros(N)
    ccf[:N/2] = correlation[N/2:]
    ccf[N/2:] = correlation[:N/2]

    h = ccf.max()
    ccf -= ccf.min()
    ccf *= h/ccf.max()

    # Fit +/- 5 pixels
    idx = np.argmax(ccf) - 3, np.argmax(ccf) + 3

    coeffs = np.polyfit(z_array[idx[0]:idx[1]], ccf[idx[0]:idx[1]], 2)
    x_i = np.linspace(z_array[idx[0]], z_array[idx[1]], 1000)
    y_i = np.polyval(coeffs, x_i)

    # Update the value
    z[index] = x_i[y_i.argmax()]

    """
    # Fit the profile peak.
    fwhm = np.ptp(z_array[np.where(ccf >= 0.5)])
    p0 = np.array([z[index], fwhm/2.355, 1.0]) # mu, sigma, peak.
    use = (z_array > (p0[0] - 3 * p0[1])) * ((p0[0] + 3*p0[1]) > z_array)
    x, y = z_array[use], ccf[use]

    # Fit profile:
    f = lambda x, mu, sigma, peak: peak * np.exp(-(mu - x)**2/(2*sigma**2)) 
    
    import matplotlib.pyplot as plt
    plt.plot(z_array, ccf, c='k')
    plt.plot(x, y, c='b')
    plt.plot(x, f(x, *p0), c='r')

    p1 = curve_fit(f, x, y, p0=p0)
    plt.plot(x, f(x, *p1[0]), c='m')

    plt.gca().set_ylim(0,1)
    raise a

    p1 = curve_fit(f, x, y, p0=p0)
    """

    return (z * c, z_err * c, R)
576185
def coefficient_of_variation(data, data_mean=None, data_stdev=None):
    
    data_mean = data_mean or mean(data)
    data_stdev = data_stdev or stdevp(data, data_mean)
    if data_mean == 0:
        return float("inf") if data_stdev != 0 else 0
    else:
        return data_stdev / data_mean
452467
def readout(self):
        """"""

        elec = self.simulate_poisson_variate()

        elec_pre = self.saturate(elec)

        elec_f = self.pre_readout(elec_pre)

        adu_r = self.base_readout(elec_f)

        adu_p = self.post_readout(adu_r)

        self.clean_up()

        return adu_p
1237663
private void calculateSCL(double[] x) {
    //System.out.println("Checking at: "+x[0]+" "+x[1]+" "+x[2]);
    value = 0.0;
    Arrays.fill(derivative, 0.0);
    double[] sums = new double[numClasses];
    double[] probs = new double[numClasses];
    double[] counts = new double[numClasses];
    Arrays.fill(counts, 0.0);
    for (int d = 0; d < data.length; d++) {
      //       if (d == testMin) {
      //         d = testMax - 1;
      //         continue;
      //       }
      int[] features = data[d];
      // activation
      Arrays.fill(sums, 0.0);
      for (int c = 0; c < numClasses; c++) {
        for (int f = 0; f < features.length; f++) {
          int i = indexOf(features[f], c);
          sums[c] += x[i];
        }
      }
      // expectation (slower routine replaced by fast way)
      // double total = Double.NEGATIVE_INFINITY;
      // for (int c=0; c<numClasses; c++) {
      //   total = SloppyMath.logAdd(total, sums[c]);
      // }
      double total = ArrayMath.logSum(sums);
      int ld = labels[d];
      for (int c = 0; c < numClasses; c++) {
        probs[c] = Math.exp(sums[c] - total);
        for (int f = 0; f < features.length; f++) {
          int i = indexOf(features[f], c);
          derivative[i] += probs[ld] * probs[c];
        }
      }
      // observed
      for (int f = 0; f < features.length; f++) {
        int i = indexOf(features[f], labels[d]);
        derivative[i] -= probs[ld];
      }
      value -= probs[ld];
    }
    // priors
    if (true) {
      for (int i = 0; i < x.length; i++) {
        double k = 1.0;
        double w = x[i];
        value += k * w * w / 2.0;
        derivative[i] += k * w;
      }
    }
  }
1321761
def correlation(a, b):
    ""
    if isinstance(a, list):
        a = np.array(a)
    if isinstance(b, list):
        b = np.array(b)
    a = a.reshape(1, -1)
    b = b.reshape(1, -1)
    return cdist(a, b, 'correlation')
821859
def simulate_measurement(self, index: int) -> bool:
        
        args = self._shard_num_args({'index': index})
        prob_one = np.sum(self._pool.map(_one_prob_per_shard, args))
        result = bool(np.random.random() <= prob_one)

        args = self._shard_num_args({
            'index': index,
            'result': result,
            'prob_one': prob_one
        })
        self._pool.map(_collapse_state, args)
        return result
653447
public int compare(Performance o1, Performance o2) {
    int	result;
    double	p1;
    double	p2;

    if (o1.getEvaluation() != o2.getEvaluation())
      throw new IllegalArgumentException("Comparing different types of performances!");

    p1 = o1.getPerformance();
    p2 = o2.getPerformance();

    if (p1 < p2)
      result = -1;
    else if (p1 > p2)
      result = 1;
    else
      result = o1.getValues().compareTo(o2.getValues());

    // only correlation coefficient/accuracy/kappa obey to this order, for the
    // errors (and the combination of all three), the smaller the number the
    // better -> hence invert them
    if (getMetrics().invert(getEvaluation()))
      result = -result;

    return result;
  }