Fisher's method for combining independent p-values.
2332
def fishers_method(pvals):
    """"""
    pvals = np.asarray(pvals)
    degrees_of_freedom = 2 * pvals.size
    chisq_stat = np.sum(-2*np.log(pvals))
    fishers_pval = stats.chi2.sf(chisq_stat, degrees_of_freedom)
    return fishers_pval
271692
def _calc_norm_coeff(h, J, bias_range, quadratic_range, ignored_variables,
                     ignored_interactions):
    """"""

    if ignored_variables is None or ignored_interactions is None:
        raise ValueError('ignored interactions or variables cannot be None')

    def parse_range(r):
        if isinstance(r, Number):
            return -abs(r), abs(r)
        return r

    def min_and_max(iterable):
        if not iterable:
            return 0, 0
        return min(iterable), max(iterable)

    if quadratic_range is None:
        linear_range, quadratic_range = bias_range, bias_range
    else:
        linear_range = bias_range

    lin_range, quad_range = map(parse_range, (linear_range,
                                              quadratic_range))

    lin_min, lin_max = min_and_max([v for k, v in h.items()
                                    if k not in ignored_variables])
    quad_min, quad_max = min_and_max([v for k, v in J.items()
                                      if not check_isin(k,
                                                        ignored_interactions)])

    inv_scalar = max(lin_min / lin_range[0], lin_max / lin_range[1],
                     quad_min / quad_range[0], quad_max / quad_range[1])

    if inv_scalar != 0:
        return 1. / inv_scalar
    else:
        return 1.
548767
def mae(x_values, y_values, drop_missing = True):
    ''''''
    num_points = len(x_values)
    assert(num_points == len(y_values) and num_points > 0)
    return numpy.sum(numpy.apply_along_axis(numpy.abs, 0, numpy.subtract(x_values, y_values))) / float(num_points)
310191
def compute_skew_stats(intra, inter):
    """"""
    # Intra (within a group) stats
    intra_skew = {}
    for k, v in intra.items():
        skew = st.skew(v)
        try:
            skew_zstat, skew_pval = st.skewtest(v)
        except ValueError:  # if sample size too small
            skew_zstat, skew_pval = (0, 1)
        intra_skew[k] = (skew, skew_zstat, skew_pval)

    # Inter (between groups) stats
    inter_skew = {}
    for k, v in inter.items():
        # Inter skew stats
        skew_sep = st.skew(v.flatten())
        try:
            skew_sep_zstat, skew_sep_pval = st.skewtest(v.flatten())
        except ValueError:
            skew_sep_zstat, skew_sep_pval = (0, 1)
        inter_skew['-'.join(k)] = (skew_sep, skew_sep_zstat, skew_sep_pval)

        # Significance of difference between intra and inter distributions
        for intra_key in k:
            try:
                separation_zstat, separation_pval = mannwhitneyu(intra[intra_key],
                                                                 v.flatten(),
                                                                 alternative='less')
            except ValueError:  # All numbers are identical in mannwhitneyu
                separation_zstat, separation_pval = (0, 1)
            inter_skew['{}<{}'.format(intra_key, '-'.join(k))] = (separation_zstat, separation_pval)

    return intra_skew, inter_skew
860038
def get_fisher_scores_vs_background(self):
        '''
        
        '''
        df = self.get_term_and_background_counts()
        odds_ratio, p_values = self._get_fisher_scores_from_counts(
            df['corpus'], df['background'])
        df['Odds ratio'] = odds_ratio
        df['Bonferroni-corrected p-values'] = p_values * len(df)
        df.sort_values(by=['Bonferroni-corrected p-values', 'Odds ratio'],
                       ascending=[True, False])
        return df
494846
private void computeAxisVariance(List<P> points) {
		int numPoints = points.size();

		for( int i = 0; i < N; i++ ) {
			mean[i] = 0;
			var[i] = 0;
		}

		// compute the mean
		for( int i = 0; i < numPoints; i++ ) {
			P p = points.get(i);

			for( int j = 0; j < N; j++ ) {
				mean[j] += distance.valueAt(p,j);
			}
		}

		for( int i = 0; i < N; i++ ) {
			mean[i] /= numPoints;
		}

		// compute the variance * N
		for( int i = 0; i < numPoints; i++ ) {
			P p = points.get(i);

			for( int j = 0; j < N; j++ ) {
				double d = mean[j] - distance.valueAt(p,j);
				var[j] += d*d;
			}
		}
	}
1349474
def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True):
    

    # Contour over a range of length scales and signal/noise ratios.
    length_scales = np.linspace(0.1, 60., resolution)
    log_SNRs = np.linspace(-3., 4., resolution)

    try:import pods
    except ImportError:
        print('pods unavailable, see https://github.com/sods/ods for example datasets')
        return
    data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
    # data['Y'] = data['Y'][0::2, :]
    # data['X'] = data['X'][0::2, :]

    data['Y'] = data['Y'] - np.mean(data['Y'])

    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)
    if plot:
        pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
        ax = pb.gca()
        pb.xlabel('length scale')
        pb.ylabel('log_10 SNR')

        xlim = ax.get_xlim()
        ylim = ax.get_ylim()

    # Now run a few optimizations
    models = []
    optim_point_x = np.empty(2)
    optim_point_y = np.empty(2)
    np.random.seed(seed=seed)
    for i in range(0, model_restarts):
        # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
        kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))

        m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern)
        m.likelihood.variance = np.random.uniform(1e-3, 1)
        optim_point_x[0] = m.rbf.lengthscale
        optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        # optimize
        if optimize:
            m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)

        optim_point_x[1] = m.rbf.lengthscale
        optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        if plot:
            pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
        models.append(m)

    if plot:
        ax.set_xlim(xlim)
        ax.set_ylim(ylim)
    return m
1212739
def reject_sv(m, s, y):
    
    mp = m + 0.5 * s**2 * (-1. + y**2 * np.exp(-m))
    ntries = 0
    while True:
        ntries += 1
        x = stats.norm.rvs(loc=mp, scale=s)
        u = stats.uniform.rvs()
        if np.log(u) < - 0.5 * y**2 * (np.exp(-x) - np.exp(-m) * (1. + m - x)):
            break
        if ntries > 1000:
            print('1000 failed attempt, m,s,y=%f, %f, %f' % (m, s, y))
            break
    return x
1355879
def sin_fisher(y,k):
    
    return 1/np.sqrt(1-y**2) * (k/(np.sinh(k)) * y * (np.cosh(k*np.sqrt(1-y**2))))
444484
def Kdiag(self, X):
        """"""
        K_ff = np.zeros(X.shape[0])
        for i,x in enumerate(X):
            K_ff[i] = self.k_ff(x[0],x[0],self.lengthscale[0])
        return K_ff * self.variances[0]
371373
def _rm_anova(self, dv=None, within=None, subject=None, detailed=False,
              correction='auto', export_filename=None):
    """"""
    aov = rm_anova(data=self, dv=dv, within=within, subject=subject,
                   correction=correction, detailed=detailed,
                   export_filename=export_filename)
    return aov
1139004
def ts_stats_significance_bootstrap(ts, stats_ts, stats_func, B=1000, b=3):
    
    pvals = []
    for tp in np.arange(0, len(stats_ts)):
        pf = partial(stats_func, t=tp)
        bs = bootstrap_ts(ts, pf, B=B, b=b)
        ci = get_ci(bs, blockratio=b / len(stats_ts))
        pval = abs(get_pvalue(stats_ts[tp], ci))
        pvals.append(pval)
    return pvals
317284
def qualNorm_filter_genes(data, qualitative, pval_threshold=0.05, smoothing=1e-5, eps=1e-5):
    
    genes, cells = data.shape
    clusters = qualitative.shape[1]
    output = np.zeros((genes, clusters))
    missing_indices = []
    genes_included = []
    qual_indices = []
    thresholds = qualitative.min(1) + (qualitative.max(1) - qualitative.min(1))/2.0
    pvals = np.zeros(genes)
    for i in range(genes):
        if qualitative[i,:].max() == -1 and qualitative[i,:].min() == -1:
            missing_indices.append(i)
            continue
        qual_indices.append(i)
        threshold = thresholds[i]
        data_i = data[i,:]
        if sparse.issparse(data):
            data_i = data_i.toarray().flatten()
        assignments, means = poisson_cluster(data_i.reshape((1, cells)), 2)
        means = means.flatten()
        high_i = 1
        low_i = 0
        if means[0]>means[1]:
            high_i = 0
            low_i = 1
        # do a p-value test
        p_val = poisson_test(data_i[assignments==low_i], data_i[assignments==high_i], smoothing=smoothing)
        pvals[i] = p_val
        if p_val <= pval_threshold:
            genes_included.append(i)
        else:
            continue
        high_mean = np.median(data_i[assignments==high_i])
        low_mean = np.median(data_i[assignments==low_i]) + eps
        for k in range(clusters):
            if qualitative[i,k]>threshold:
                output[i,k] = high_mean
            else:
                output[i,k] = low_mean
    output = output[genes_included,:]
    pvals = pvals[genes_included]
    return output, pvals, genes_included
864897
private static DataTable2D bivariateMatrix(Dataframe dataSet, BivariateType type) {
        DataTable2D bivariateMatrix = new DataTable2D();
        
        //extract values of first variable
        Map<Object, TypeInference.DataType> columnTypes = dataSet.getXDataTypes();
        Object[] allVariables = columnTypes.keySet().toArray();
        int numberOfVariables = allVariables.length;
        
        TransposeDataList transposeDataList;
        for(int i=0;i<numberOfVariables;++i) {
            Object variable0 = allVariables[i];
            if(columnTypes.get(variable0)==TypeInference.DataType.CATEGORICAL) {
                continue;
            }
            
            transposeDataList = new TransposeDataList();
            
            //extract values of first variable
            transposeDataList.put(0, dataSet.getXColumn(variable0));
            
            for(int j=i;j<numberOfVariables;++j) {
                Object variable1 = allVariables[j];
                if(columnTypes.get(variable1)==TypeInference.DataType.CATEGORICAL) {
                    continue;
                }
            
                transposeDataList.put(1, dataSet.getXColumn(variable1));
                
                double value = 0.0;
                if(type==BivariateType.COVARIANCE) {
                    value = Descriptives.covariance(transposeDataList, true);
                }
                else if(type==BivariateType.PEARSONCORRELATION) {
                    if(variable0.equals(variable1)) {
                        value=1.0;
                    }
                    else {
                        value = PearsonCorrelation.calculateCorrelation(transposeDataList);
                    }
                }
                else if(type==BivariateType.SPEARMANCORRELATION) {
                    if(variable0.equals(variable1)) {
                        value=1.0;
                    }
                    else {
                        value = SpearmanCorrelation.calculateCorrelation(transposeDataList);
                    }
                }
                else if(type==BivariateType.KENDALLTAUCORRELATION) {
                    if(variable0.equals(variable1)) {
                        value=1.0;
                    }
                    else {
                        value = KendallTauCorrelation.calculateCorrelation(transposeDataList);
                    }
                }
                
                //bivariateMatrix.internalData.get(variable0).internalData.put(variable1, value);
                bivariateMatrix.put2d(variable0, variable1, value);
                
                if(!variable0.equals(variable1)) {
                    /*
                    if(!bivariateMatrix.internalData.containsKey(variable1)) {
                        bivariateMatrix.internalData.put(variable1, new AssociativeArray());
                    }
                    bivariateMatrix.internalData.get(variable1).internalData.put(variable0, value);
                    */
                    bivariateMatrix.put2d(variable1, variable0, value);
                }
            }
            //transposeDataList = null;
        }
        
        return bivariateMatrix;
    }
511391
def gpinv(p, k, sigma):
    """"""
    x = np.empty(p.shape)
    x.fill(np.nan)
    if sigma <= 0:
        return x
    ok = (p > 0) & (p < 1)
    if np.all(ok):
        if np.abs(k) < np.finfo(float).eps:
            np.negative(p, out=x)
            np.log1p(x, out=x)
            np.negative(x, out=x)
        else:
            np.negative(p, out=x)
            np.log1p(x, out=x)
            x *= -k
            np.expm1(x, out=x)
            x /= k
        x *= sigma
    else:
        if np.abs(k) < np.finfo(float).eps:
            # x[ok] = - np.log1p(-p[ok])
            temp = p[ok]
            np.negative(temp, out=temp)
            np.log1p(temp, out=temp)
            np.negative(temp, out=temp)
            x[ok] = temp
        else:
            # x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
            temp = p[ok]
            np.negative(temp, out=temp)
            np.log1p(temp, out=temp)
            temp *= -k
            np.expm1(temp, out=temp)
            temp /= k
            x[ok] = temp
        x *= sigma
        x[p == 0] = 0
        if k >= 0:
            x[p == 1] = np.inf
        else:
            x[p == 1] = -sigma / k
    return x
692432
def _exec_vector(self, a, bd, xy, xy_orig, mask,
                     n_withdrifts, spec_drift_grids):
        """"""

        npt = bd.shape[0]
        n = self.X_ADJUSTED.shape[0]
        zero_index = None
        zero_value = False

        a_inv = scipy.linalg.inv(a)

        if np.any(np.absolute(bd) <= self.eps):
            zero_value = True
            zero_index = np.where(np.absolute(bd) <= self.eps)

        if self.UNBIAS:
            b = np.zeros((npt, n_withdrifts+1, 1))
        else:
            b = np.zeros((npt, n_withdrifts, 1))
        b[:, :n, 0] = - self.variogram_function(self.variogram_model_parameters, bd)
        if zero_value:
            b[zero_index[0], zero_index[1], 0] = 0.0

        i = n
        if self.regional_linear_drift:
            b[:, i, 0] = xy[:, 0]
            i += 1
            b[:, i, 0] = xy[:, 1]
            i += 1
        if self.point_log_drift:
            for well_no in range(self.point_log_array.shape[0]):
                log_dist = np.log(np.sqrt((xy[:, 0] - self.point_log_array[well_no, 0])**2 +
                                          (xy[:, 1] - self.point_log_array[well_no, 1])**2))
                if np.any(np.isinf(log_dist)):
                    log_dist[np.isinf(log_dist)] = -100.0
                b[:, i, 0] = - self.point_log_array[well_no, 2] * log_dist
                i += 1
        if self.external_Z_drift:
            b[:, i, 0] = self._calculate_data_point_zscalars(xy_orig[:, 0], xy_orig[:, 1])
            i += 1
        if self.specified_drift:
            for spec_vals in spec_drift_grids:
                b[:, i, 0] = spec_vals.flatten()
                i += 1
        if self.functional_drift:
            for func in self.functional_drift_terms:
                b[:, i, 0] = func(xy[:, 0], xy[:, 1])
                i += 1
        if i != n_withdrifts:
            warnings.warn("Error in setting up kriging system. "
                          "Kriging may fail.", RuntimeWarning)
        if self.UNBIAS:
            b[:, n_withdrifts, 0] = 1.0

        if (~mask).any():
            mask_b = np.repeat(mask[:, np.newaxis, np.newaxis],
                               n_withdrifts+1, axis=1)
            b = np.ma.array(b, mask=mask_b)

        if self.UNBIAS:
            x = np.dot(a_inv, b.reshape((npt, n_withdrifts+1)).T).reshape((1, n_withdrifts+1, npt)).T
        else:
            x = np.dot(a_inv, b.reshape((npt, n_withdrifts)).T).reshape((1, n_withdrifts, npt)).T
        zvalues = np.sum(x[:, :n, 0] * self.Z, axis=1)
        sigmasq = np.sum(x[:, :, 0] * -b[:, :, 0], axis=1)

        return zvalues, sigmasq
398579
def fit(ts, fs=[], all_params=[], fit_vars=None,
        alg='leastsq', make_bounded=True):
    
    if fit_vars is None:
        fit_vars = [f._peakargs for f in fs]
    initc = [min(ts.values)]
    for f, peak_params, to_fit in zip(fs, all_params, fit_vars):
        if 'v' in to_fit:
            to_fit.remove('v')

        if make_bounded and hasattr(f, '_pbounds'):
            new_v = _to_unbnd_p({i: peak_params[i] for i in to_fit},
                                f._pbounds)
            initc += [new_v[i] for i in to_fit]
        else:
            initc += [peak_params[i] for i in to_fit]

    def errfunc_lsq(fit_params, t, y, all_params):
        # first value in fit_params is baseline
        # fit_y = np.ones(len(t)) * fit_params[0]
        fit_y = np.zeros(len(t))
        param_i = 1
        for f, peak_params, to_fit in zip(fs, all_params, fit_vars):
            for k in to_fit:
                peak_params[k] = fit_params[param_i]
                param_i += 1
            if make_bounded and hasattr(f, '_pbounds'):
                fit_y += f(t, **_to_bound_p(peak_params, f._pbounds))
            else:
                fit_y += f(t, **peak_params)
        return fit_y - y

    def errfunc(p, t, y, all_params):
        return np.sum(errfunc_lsq(p, t, y, all_params) ** 2)

    if alg == 'simplex':
        fit_p, _ = fmin(errfunc, initc, args=(ts.index, ts.values,
                                              peak_params))
#    elif alg == 'anneal':
#        fit_p, _ = anneal(errfunc, initc, args=(ts.index, ts.values,
#                                                peak_params))
    elif alg == 'lbfgsb':
        # TODO: use bounds param
        fitp, _ = fmin_l_bfgs_b(errfunc, fit_p,
                                args=(ts.index, ts.values, peak_params),
                                approx_grad=True)
    elif alg == 'leastsq':
        fit_p, _ = leastsq(errfunc_lsq, initc, args=(ts.index, ts.values,
                                                     all_params))
    # else:
    #     r = minimize(errfunc, initc, \
    #                  args=(ts.index, ts.values, all_params), \
    #                  jac=False, gtol=1e-2)
    #     #if not r['success']:
    #     #    print('Fail:' + str(f))
    #     #    print(r)
    #     #if np.nan in r['x']:  # not r['success']?
    #     #    fit_p = initc
    #     #else:
    #     #    fit_p = r['x']

    fit_pl = fit_p.tolist()
    v = fit_pl.pop(0)  # noqa
    fitted_params = []
    for f, to_fit in zip(fs, fit_vars):
        fit_p_dict = {v: fit_pl.pop(0) for v in to_fit}
        # fit_p_dict['v'] = v
        if make_bounded and hasattr(f, '_pbounds'):
            fitted_params.append(_to_bound_p(fit_p_dict, f._pbounds))
        else:
            fitted_params.append(fit_p_dict)

    # calculate r^2 of the fit
    ss_err = errfunc(fit_p, ts.index, ts.values, fitted_params)
    ss_tot = np.sum((ts.values - np.mean(ts.values)) ** 2)
    r2 = 1 - ss_err / ss_tot
    res = {'r^2': r2}

    return fitted_params, res