Non-central Student's t random variates.
546561
def rnoncentral_t(mu, lam, nu, size=None):
    
    tau = rgamma(nu / 2., nu / (2. * lam), size)
    return rnormal(mu, tau)
275293
def inference(self, kern, X, Z, likelihood, Y, indexD, output_dim, Y_metadata=None, Lm=None, dL_dKmm=None, Kuu_sigma=None):
        

        input_dim = Z.shape[0]

        uncertain_inputs = isinstance(X, VariationalPosterior)

        beta = 1./likelihood.variance
        if len(beta)==1:
            beta = np.zeros(output_dim)+beta

        beta_exp = np.zeros(indexD.shape[0])
        for d in range(output_dim):
            beta_exp[indexD==d] = beta[d]

        psi0, psi1, psi2 = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

        psi2_sum = (beta_exp[:,None,None]*psi2).sum(0)/output_dim

        #======================================================================
        # Compute Common Components
        #======================================================================

        Kmm = kern.K(Z).copy()
        if Kuu_sigma is not None:
            diag.add(Kmm, Kuu_sigma)
        else:
            diag.add(Kmm, self.const_jitter)
        Lm = jitchol(Kmm)

        logL = 0.
        dL_dthetaL = np.zeros(output_dim)
        dL_dKmm = np.zeros_like(Kmm)
        dL_dpsi0 = np.zeros_like(psi0)
        dL_dpsi1 = np.zeros_like(psi1)
        dL_dpsi2 = np.zeros_like(psi2)
        wv = np.empty((Kmm.shape[0],output_dim))

        for d in range(output_dim):
            idx_d = indexD==d
            Y_d = Y[idx_d]
            N_d = Y_d.shape[0]
            beta_d = beta[d]

            psi2_d = psi2[idx_d].sum(0)*beta_d
            psi1Y = Y_d.T.dot(psi1[idx_d])*beta_d
            psi0_d = psi0[idx_d].sum()*beta_d
            YRY_d = np.square(Y_d).sum()*beta_d

            LmInvPsi2LmInvT = backsub_both_sides(Lm, psi2_d, 'right')

            Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
            LL = jitchol(Lambda)
            LmLL = Lm.dot(LL)

            b  = dtrtrs(LmLL, psi1Y.T)[0].T
            bbt = np.square(b).sum()
            v = dtrtrs(LmLL, b.T, trans=1)[0].T
            LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)

            tmp = -backsub_both_sides(LL, LLinvPsi1TYYTPsi1LLinvT)
            dL_dpsi2R = backsub_both_sides(Lm, tmp+np.eye(input_dim))/2

            logL_R = -N_d*np.log(beta_d)
            logL += -((N_d*log_2_pi+logL_R+psi0_d-np.trace(LmInvPsi2LmInvT))+YRY_d- bbt)/2.

            dL_dKmm +=  dL_dpsi2R - backsub_both_sides(Lm, LmInvPsi2LmInvT)/2

            dL_dthetaL[d:d+1] = (YRY_d*beta_d + beta_d*psi0_d - N_d*beta_d)/2. - beta_d*(dL_dpsi2R*psi2_d).sum() - beta_d*np.trace(LLinvPsi1TYYTPsi1LLinvT)

            dL_dpsi0[idx_d] = -beta_d/2.
            dL_dpsi1[idx_d] = beta_d*np.dot(Y_d,v)
            dL_dpsi2[idx_d] = beta_d*dL_dpsi2R
            wv[:,d] = v

        LmInvPsi2LmInvT = backsub_both_sides(Lm, psi2_sum, 'right')

        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LmLL = Lm.dot(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        dL_dpsi2R_common = dpotri(LmLL)[0]/-2.
        dL_dpsi2 += dL_dpsi2R_common[None,:,:]*beta_exp[:,None,None]

        for d in range(output_dim):
            dL_dthetaL[d] += (dL_dpsi2R_common*psi2[indexD==d].sum(0)).sum()*-beta[d]*beta[d]

        dL_dKmm += dL_dpsi2R_common*output_dim

        logL += -output_dim*logdet_L/2.

        #======================================================================
        # Compute dL_dKmm
        #======================================================================

        # dL_dKmm =  dL_dpsi2R - output_dim* backsub_both_sides(Lm, LmInvPsi2LmInvT)/2 #LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

        #======================================================================
        # Compute the Posterior distribution of inducing points p(u|Y)
        #======================================================================

        LLInvLmT = dtrtrs(LL, Lm.T)[0]
        cov = tdot(LLInvLmT.T)

        wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
        post = Posterior(woodbury_inv=wd_inv, woodbury_vector=wv, K=Kmm, mean=None, cov=cov, K_chol=Lm)

        #======================================================================
        # Compute dL_dthetaL for uncertian input and non-heter noise
        #======================================================================

        # for d in range(output_dim):
        #     dL_dthetaL[d:d+1] += - beta[d]*beta[d]*(dL_dpsi2R[None,:,:] * psi2[indexD==d]/output_dim).sum()
        # dL_dthetaL += - (dL_dpsi2R[None,:,:] * psi2_sum*D beta*(dL_dpsi2R*psi2).sum()

        #======================================================================
        # Compute dL_dpsi
        #======================================================================

        if not uncertain_inputs:
            dL_dpsi1 += (psi1[:,None,:]*dL_dpsi2).sum(2)*2.

        if uncertain_inputs:
            grad_dict = {'dL_dKmm': dL_dKmm,
                         'dL_dpsi0':dL_dpsi0,
                         'dL_dpsi1':dL_dpsi1,
                         'dL_dpsi2':dL_dpsi2,
                         'dL_dthetaL':dL_dthetaL}
        else:
            grad_dict = {'dL_dKmm': dL_dKmm,
                         'dL_dKdiag':dL_dpsi0,
                         'dL_dKnm':dL_dpsi1,
                         'dL_dthetaL':dL_dthetaL}

        return post, logL, grad_dict
624290
def sample_posterior(x,post,nsamples=1):
    
    cdf = post.cumsum()
    cdf /= cdf.max()
    u = rand.random(size=nsamples)
    inds = np.digitize(u,cdf)
    return x[inds]
942048
def _rand_init(x_bounds, x_types, selection_num_starting_points):
    '''
    
    '''
    return [lib_data.rand(x_bounds, x_types) for i \
                    in range(0, selection_num_starting_points)]
1222157
def RandomStandardNormal(shape, dtype, seed):
    
    if seed:
        np.random.seed(seed)
    return np.random.normal(size=reduce(mul, shape)).reshape(shape).astype(dtype_map[dtype]),
419856
def _distance_sqr_stats_naive_generic(x, y, matrix_centered, product,
                                      exponent=1):
    """"""
    a = matrix_centered(x, exponent=exponent)
    b = matrix_centered(y, exponent=exponent)

    covariance_xy_sqr = product(a, b)
    variance_x_sqr = product(a, a)
    variance_y_sqr = product(b, b)

    denominator_sqr = np.absolute(variance_x_sqr * variance_y_sqr)
    denominator = _sqrt(denominator_sqr)

    # Comparisons using a tolerance can change results if the
    # covariance has a similar order of magnitude
    if denominator == 0.0:
        correlation_xy_sqr = 0.0
    else:
        correlation_xy_sqr = covariance_xy_sqr / denominator

    return Stats(covariance_xy=covariance_xy_sqr,
                 correlation_xy=correlation_xy_sqr,
                 variance_x=variance_x_sqr,
                 variance_y=variance_y_sqr)
1194758
def build_sample_ace_problem_wang04(N=100):
    """"""
    x = [numpy.random.uniform(-1, 1, size=N)
         for _i in range(0, 5)]
    noise = numpy.random.standard_normal(N)
    y = numpy.log(4.0 + numpy.sin(4 * x[0]) + numpy.abs(x[1]) + x[2] ** 2 +
                  x[3] ** 3 + x[4] + 0.1 * noise)
    return x, y
1044394
def build_sample_ace_problem_breiman85(N=200):
    """"""
    x_cubed = numpy.random.standard_normal(N)
    x = scipy.special.cbrt(x_cubed)
    noise = numpy.random.standard_normal(N)
    y = numpy.exp((x ** 3.0) + noise)
    return [x], y
946009
def _random_x(self):
        """"""
        return (tuple(random.random() for _ in range(self.fmodel.dim_x)),)
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
1164754
def K(self, X, X2, target):
        """"""
        if X2 is None: X2 = X
       # i1 = X[:,1]
       # i2 = X2[:,1]
       # X = X[:,0].reshape(-1,1)
       # X2 = X2[:,0].reshape(-1,1)
        dist = np.abs(X - X2.T)
        
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU

        k1 = np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
        k2 = (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
        k3 = np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
943878
function fuzzy (n, pct) {
  var variance = n * pct * 2
  return Math.floor(n + (Math.random() * variance - variance / 2))
}
317853
def point_predict(f, x, size=1, nugget=None):
    
    orig_shape = x.shape
    x = regularize_array(x)

    mu = f.M_internal(x, regularize=False)
    V = f.C_internal(x, regularize=False)
    if nugget is not None:
        V += nugget
    out = random.normal(size=(size, x.shape[0])) * sqrt(V) + mu
    return out.reshape((size,) + orig_shape[:-1]).squeeze()
753850
def _random_stochastic_matrix(m, n, k=None, sparse=False, format='csr',
                              random_state=None):
    
    if k is None:
        k = n
    # m prob vectors of dimension k, shape (m, k)
    probvecs = probvec(m, k, random_state=random_state)

    if k == n:
        P = probvecs
        if sparse:
            return scipy.sparse.coo_matrix(P).asformat(format)
        else:
            return P

    # if k < n:
    rows = np.repeat(np.arange(m), k)
    cols = \
        sample_without_replacement(
            n, k, num_trials=m, random_state=random_state
        ).ravel()
    data = probvecs.ravel()

    if sparse:
        P = scipy.sparse.coo_matrix((data, (rows, cols)), shape=(m, n))
        return P.asformat(format)
    else:
        P = np.zeros((m, n))
        P[rows, cols] = data
        return P
164269
public void computeNccDescriptor( NccFeature f , float x0 , float y0 , float x1 , float y1 ) {
		double mean = 0;
		float widthStep = (x1-x0)/15.0f;
		float heightStep = (y1-y0)/15.0f;

		// compute the mean value
		int index = 0;
		for( int y = 0; y < 15; y++ ) {
			float sampleY = y0 + y*heightStep;
			for( int x = 0; x < 15; x++ ) {
				mean += f.value[index++] = interpolate.get_fast(x0 + x * widthStep, sampleY);
			}
		}
		mean /= 15*15;

		// compute the variance and save the difference from the mean
		double variance = 0;
		index = 0;
		for( int y = 0; y < 15; y++ ) {
			for( int x = 0; x < 15; x++ ) {
				double v = f.value[index++] -= mean;
				variance += v*v;
			}
		}
		variance /= 15*15;
		f.mean = mean;
		f.sigma = Math.sqrt(variance);
	}
37372
def rand(x_bounds, x_types):
    '''
    
    '''
    outputs = []

    for i, _ in enumerate(x_bounds):
        if x_types[i] == "discrete_int":
            temp = x_bounds[i][random.randint(0, len(x_bounds[i]) - 1)]
            outputs.append(temp)
        elif x_types[i] == "range_int":
            temp = random.randint(x_bounds[i][0], x_bounds[i][1])
            outputs.append(temp)
        elif x_types[i] == "range_continuous":
            temp = random.uniform(x_bounds[i][0], x_bounds[i][1])
            outputs.append(temp)
        else:
            return None

    return outputs
111797
def _init_random_gaussians(self, X):
        """"""

        n_samples = np.shape(X)[0]
        self.priors = (1 / self.k) * np.ones(self.k)
        for _ in range(self.k):
            params = {}
            params["mean"] = X[np.random.choice(range(n_samples))]
            params["cov"] = calculate_covariance_matrix(X)
            self.parameters.append(params)