Draw parameters from a mean-field normal family
431502
def draw_normal(self):
        
        means, scale = self.get_means_and_scales()
        return np.random.normal(means,scale,size=[self.sims,means.shape[0]]).T
247081
def calc_average_parameters(parameter_layers):
        
        mean_layers = [numpy.mean(x) if x[0] else 0 for x in parameter_layers]
        overall_mean = numpy.mean([x for x in mean_layers if x])
        return mean_layers, overall_mean
1114578
def draw(self, **kwargs):
        
        # Compute the curves
        x = self.n_feature_subsets_
        means = self.cv_scores_.mean(axis=1)
        sigmas = self.cv_scores_.std(axis=1)


        # Plot one standard deviation above and below the mean
        self.ax.fill_between(x, means - sigmas, means+sigmas, alpha=0.25)

        # Plot the curve
        self.ax.plot(x, means, 'o-')

        # Plot the maximum number of features
        self.ax.axvline(
            self.n_features_, c='k', ls='--',
            label="n_features = {}\nscore = {:0.3f}".format(
                self.n_features_, self.cv_scores_.mean(axis=1).max()
            )
        )

        return self.ax
17899
def normal_log_q_initial(self,z):
                     
        means, scale = self.get_means_and_scales_from_q()
        return ss.norm.logpdf(z,loc=means,scale=scale)
571195
def fit_mle(self, data, fix_mean=False):
        

        if not fix_mean:
            sigma, _, scale = stats.lognorm.fit(data, floc=0)
            return np.log(scale), sigma

        else:
            mean = np.mean(data)

            # MLE fxn to be optmimized
            mle = lambda sigma, x, mean: -1 *\
                                    np.sum(self._pdf_w_mean(x, mean, sigma))

            sigma = optim.fmin(mle, np.array([np.std(np.log(data), ddof=1)]),
                                            args=(data, mean), disp=0)[0]

            return self.translate_args(mean, sigma)
947035
def _synDelayParams(self):
        '''
        
        '''
        delays = {}
        #mean delays
        loc = np.zeros((len(self.y), len(self.X)))
        loc[:, 0] = self.delays[0]
        loc[:, 1::2] = self.delays[0]
        loc[:, 2::2] = self.delays[1]
        #standard deviations
        scale = loc * self.delay_rel_sd
        
        #prepare output
        delay_loc = {}
        for i, y in enumerate(self.y):
            delay_loc.update({y : loc[i]})
        
        delay_scale = {}
        for i, y in enumerate(self.y):
            delay_scale.update({y : scale[i]})
                
        return delay_loc, delay_scale
735952
def draw(self):
        '''
        
        '''
        observed_arr = self.extract_conditions()
        conv_arr = self.inference(observed_arr)

        if self.__conditon_noise_sampler is not None:
            self.__conditon_noise_sampler.output_shape = conv_arr.shape
            noise_arr = self.__conditon_noise_sampler.generate()
            conv_arr += noise_arr

        deconv_arr = self.__deconvolution_model.inference(conv_arr)
        return np.concatenate((deconv_arr, observed_arr), axis=1)
567051
def _get_warped_mean(self, mean, std, pred_init=None, deg_gauss_hermite=20):
        
        gh_samples, gh_weights = np.polynomial.hermite.hermgauss(deg_gauss_hermite)
        gh_samples = gh_samples[:, None]
        gh_weights = gh_weights[None, :]
        return gh_weights.dot(self._get_warped_term(mean, std, gh_samples)) / np.sqrt(np.pi)
731013
def getCovariance(self,normalize=True,i0=None,i1=None,pos0=None,pos1=None,chrom=None,center=True,unit=True,pos_cum0=None,pos_cum1=None,blocksize=None,X=None,**kw_args):
        """"""
        if X is not None:
            K=X.dot(X.T)
            Nsnp=X.shape[1]
        else:
            if (i0 is None) and (i1 is None) and ((pos0 is not None) & (pos1 is not None) & (chrom is not None)) or ((pos_cum0 is not None) & (pos_cum1 is not None)):
                i0,i1=self.getGenoIndex(pos0=pos0,pos1=pos1,chrom=chrom,pos_cum0=pos_cum0,pos_cum1=pose_cum1)

            [N,M]=self.genoM.shape
            if blocksize is None:
                blocksize=M
            if i0 is None:
                i0=0
            if i1 is None:
                i1=M
            nread = i0
            K=None
            Nsnp=i1-i0
            while nread<i1:
                thisblock=min(blocksize,i1-nread)
                X=self.getGenotypes(i0=nread,i1=(nread+thisblock),center=center,unit=unit,**kw_args)    
                if K is None:
                    K=X.dot(X.T)
                else:
                    K+=X.dot(X.T)
                nread+=thisblock
        if normalize:
            K/=(K.diagonal().mean())
        else:#divide by number of SNPs in K
            K/=Nsnp
        return K
459067
def _Aterm(N, h, m, k, dW):
    """"""
    sqrt2h = np.sqrt(2.0/h)
    Xk = np.random.normal(0.0, 1.0, (N, m, 1))
    Yk = np.random.normal(0.0, 1.0, (N, m, 1))
    term1 = _dot(Xk, _t(Yk + sqrt2h*dW))
    term2 = _dot(Yk + sqrt2h*dW, _t(Xk))
    return (term1 - term2)/k
832942
def gen_mu(K, delta, c):
    

    S = c * log(K/delta) * sqrt(K) 
    tau = gen_tau(S, K, delta)
    rho = gen_rho(K)
    normalizer = sum(rho) + sum(tau)
    return [(rho[d] + tau[d])/normalizer for d in range(K)]
297607
def lognormcdf(x, mu, tau):
    """"""
    x = np.atleast_1d(x)
    return np.array(
        [0.5 * (1 - flib.derf(-(np.sqrt(tau / 2)) * (np.log(y) - mu))) for y in x])
998991
def mean_fill(adf):
        
        ordpt = adf.values[0]
        if not pd.isnull(ordpt):
            return ordpt

        fdmn = adf.iloc[1:-1].mean()
        if not pd.isnull(fdmn):
            return fdmn

        flspt = adf.values[-1]
        if not pd.isnull(flspt):
            return flspt

        return nan
256818
def _get_mean(self, sites, C, ln_y_ref, exp1, exp2):
        
        # we do not support estimating of basin depth and instead
        # rely on it being available (since we require it).
        # centered_z1pt0
        centered_z1pt0 = self._get_centered_z1pt0(sites)
        # we consider random variables being zero since we want
        # to find the exact mean value.
        eta = epsilon = 0.

        ln_y = (
            # first line of eq. 12
            ln_y_ref + eta
            # second line
            + C['phi1'] * np.log(sites.vs30 / 1130).clip(-np.inf, 0)
            # third line
            + C['phi2'] * (exp1 - exp2)
            * np.log((np.exp(ln_y_ref) * np.exp(eta) + C['phi4']) / C['phi4'])
            # fourth line
            + C['phi5']
            * (1.0 - np.exp(-1. * centered_z1pt0 / C['phi6']))
            # fifth line
            + epsilon
        )

        return ln_y
272257
def update_H(self, mean_field, l):
        """"""

        self.H_s = self.spin_hamiltonian(mean_field, l)
        try:
            self.eig_energies, self.eig_states = diagonalize(self.H_s)
        except np.linalg.linalg.LinAlgError:
            np.savez('errorhamil', H=self.H_s, fiel=mean_field, lamb=l)
            raise
        except ValueError:
            np.savez('errorhamil', H=self.H_s, fiel=mean_field, lamb=l)
            print(mean_field, l)
            raise
980034
@Reference(authors = "D. Pelleg, A. Moore", //
      title = "X-means: Extending K-means with Efficient Estimation on the Number of Clusters", //
      booktitle = "Proc. 17th Int. Conf. on Machine Learning (ICML 2000)", //
      url = "http://www.pelleg.org/shared/hp/download/xmeans.ps", //
      bibkey = "DBLP:conf/icml/PellegM00")
  public static <V extends NumberVector> double logLikelihood(Relation<V> relation, Clustering<? extends MeanModel> clustering, NumberVectorDistanceFunction<? super V> distanceFunction) {
    List<? extends Cluster<? extends MeanModel>> clusters = clustering.getAllClusters();
    // number of clusters
    final int m = clusters.size();

    // number of objects in the clustering
    int n = 0;
    // cluster sizes
    int[] n_i = new int[m];
    // total variance
    double d = 0.;
    // variances
    double[] d_i = new double[m];

    // Iterate over clusters:
    Iterator<? extends Cluster<? extends MeanModel>> it = clusters.iterator();
    for(int i = 0; it.hasNext(); ++i) {
      Cluster<? extends MeanModel> cluster = it.next();
      n += n_i[i] = cluster.size();
      d += d_i[i] = varianceOfCluster(cluster, distanceFunction, relation);
    }

    // No remaining variance, if every point is on its own:
    if(n <= m) {
      return Double.NEGATIVE_INFINITY;
    }

    // Total variance (corrected for bias)
    final double logv = FastMath.log(d / (n - m));

    final int dim = RelationUtil.dimensionality(relation);
    // log likelihood of this clustering
    double logLikelihood = 0.;

    // Aggregate
    for(int i = 0; i < m; i++) {
      logLikelihood += n_i[i] * FastMath.log(n_i[i]) // Post. entropy Rn log Rn
          - n_i[i] * .5 * MathUtil.LOGTWOPI // Rn/2 log2pi
          - n_i[i] * dim * .5 * logv // Rn M/2 log sigma^2
          - (d_i[i] - m) * .5; // (Rn-K)/2
    }
    logLikelihood -= n * FastMath.log(n); // Prior entropy, sum_i Rn log R
    return logLikelihood;
  }
316518
def mse(X):
    '''  '''
    return np.mean(np.square(np.abs(np.fft.fft(X, axis=1))), axis=1)