Random normal variates.
329934
def rnormal(mu, tau, size=None):
    
    return np.random.normal(mu, 1. / np.sqrt(tau), size)
1029800
def normal_h2(size: int = 10000) -> Histogram2D:
    
    data1 = np.random.normal(0, 1, (size,))
    data2 = np.random.normal(0, 1, (size,))
    return h2(data1, data2, name="normal", axis_names=tuple("xy"), title="2D normal distribution")
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]),
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;
}
799499
def _log_multivariate_normal_density_spherical(X, means, covars):
    """"""
    cv = covars.copy()
    if covars.ndim == 1:
        cv = cv[:, np.newaxis]
    if cv.shape[1] == 1:
        cv = np.tile(cv, (1, X.shape[-1]))
    return _log_multivariate_normal_density_diag(X, means, cv)
357596
def resample(self,N,**kwargs):
        
        lovals = self.mu - np.absolute(rand.normal(size=N)*self.siglo)
        hivals = self.mu + np.absolute(rand.normal(size=N)*self.sighi)

        u = rand.random(size=N)
        hi = (u < float(self.sighi)/(self.sighi + self.siglo))
        lo = (u >= float(self.sighi)/(self.sighi + self.siglo))

        vals = np.zeros(N)
        vals[hi] = hivals[hi]
        vals[lo] = lovals[lo]
        return vals
709417
def _optm(x, alpha, mu, sigma):
    '''
    '''
    return ss.skewnorm.pdf(x, alpha, mu, sigma)
283293
public void setRandomStandardNormal() {
        FnIntDoubleToDouble lambda = new FnIntDoubleToDouble() {
            @Override
            public double call(int idx, double val) {
                return Gaussian.nextDouble(0.0, 1.0);
            }
        };
        apply(lambda);
    }
379393
def _gauss(mean: int, sigma: int) -> int:
        
        return int(random.gauss(mean, sigma))
407197
def random_normal(attrs, inputs, proto_obj):
    """"""
    try:
        from onnx.mapping import TENSOR_TYPE_TO_NP_TYPE
    except ImportError:
        raise ImportError("Onnx and protobuf need to be installed. "
                          "Instructions to install - https://github.com/onnx/onnx")
    new_attr = translation_utils._remove_attributes(attrs, ['seed'])
    new_attr = translation_utils._fix_attribute_names(new_attr, {'mean': 'loc'})
    new_attr['dtype'] = TENSOR_TYPE_TO_NP_TYPE[int(new_attr.get('dtype', 1))]
    return 'random_normal', new_attr, inputs
712330
def build_sample_ace_problem_breiman2(N=500):
    """"""
    x = numpy.linspace(0, 1, N)
    # x = numpy.random.uniform(0, 1, size=N)
    noise = numpy.random.standard_normal(N)
    y = numpy.exp(numpy.sin(2 * numpy.pi * x)) + 0.0 * noise
    return [x], y
786159
def _sample_weight(kappa, dim):
    
    dim = dim - 1  # since S^{n-1}
    b = dim / (np.sqrt(4. * kappa ** 2 + dim ** 2) + 2 * kappa)
    x = (1. - b) / (1. + b)
    c = kappa * x + dim * np.log(1 - x ** 2)

    while True:
        z = np.random.beta(dim / 2., dim / 2.)
        w = (1. - (1. + b) * z) / (1. - (1. - b) * z)
        u = np.random.uniform(low=0, high=1)
        if kappa * w + dim * np.log(1. - x * w) - c >= np.log(u):
            return w
1209201
def _rndPointDisposition(dx, dy):
        """"""
        x = int(random.uniform(-dx, dx))
        y = int(random.uniform(-dy, dy))
        return (x, y)
170225
def normal_equation(x, y):
    
    num_samples = y.size
    x = np.insert(x, 0, np.ones(num_samples), axis=1)

    return pinv(x)*y
478140
def simulation(n_obs=300, n_vars=None, alpha=None, beta=None, gamma=None, alpha_=None, t_max=None,
               noise_model='normal', noise_level=1, switches=None):
    
    from .tools.dynamical_model_utils import unspliced, spliced, vectorize

    def draw_poisson(n):
        from random import uniform  # draw from poisson
        t = np.cumsum([- .1 * np.log(uniform(0, 1)) for _ in range(n - 1)])
        return np.insert(t, 0, 0)  # prepend t0=0

    def simulate_dynamics(tau, alpha, beta, gamma, u0, s0, noise_model, noise_level):
        ut = unspliced(tau, u0, alpha, beta)
        st = spliced(tau, s0, u0, alpha, beta, gamma)

        if noise_model is 'normal':  # add noise
            ut += np.random.normal(scale=noise_level * np.percentile(ut, 99) / 10, size=len(ut))
            st += np.random.normal(scale=noise_level * np.percentile(st, 99) / 10, size=len(st))
        ut, st = np.clip(ut, 0, None), np.clip(st, 0, None)
        return ut, st

    def simulate_gillespie(alpha, beta, gamma):
        # update rules: transcription (u+1,s), splicing (u-1,s+1), degradation (u,s-1), nothing (u,s)
        update_rule = np.array([[1, 0], [-1, 1], [0, -1], [0, 0]])

        def update(props):
            if np.sum(props) > 0:
                props /= np.sum(props)
            p_cumsum = props.cumsum()
            p = np.random.rand()
            i = 0
            while p > p_cumsum[i]:
                i += 1
            return update_rule[i]

        u, s = np.zeros(len(alpha)), np.zeros(len(alpha))
        for i, alpha_i in enumerate(alpha):
            u_, s_ = (u[i - 1], s[i - 1]) if i > 0 else (0, 0)
            du, ds = update(props=np.array([alpha_i, beta * u_, gamma * s_]))
            u[i], s[i] = (u_ + du, s_ + ds)
        return u, s

    alpha = 5 if alpha is None else alpha
    beta = .5 if beta is None else beta
    gamma = .3 if gamma is None else gamma
    alpha_ = 0 if alpha_ is None else alpha_

    t = draw_poisson(n_obs)
    if t_max is not None:
        t *= t_max / np.max(t)
    t_max = np.max(t)

    def cycle(array, n_vars=None):
        if isinstance(array, (list, tuple)):
            return array if n_vars is None else array * int(np.ceil(n_vars / len(array)))
        else:
            return [array] if n_vars is None else [array] * n_vars

    # switching time point obtained as fraction of t_max rounded down
    switches = cycle([.4, .7, 1, .1], n_vars) if switches is None else cycle(switches, n_vars)
    t_ = np.array([np.max(t[t < t_i * t_max]) for t_i in switches])

    noise_level = cycle(noise_level, len(switches) if n_vars is None else n_vars)

    n_vars = max(len(switches), len(noise_level))
    U = np.zeros(shape=(len(t), n_vars))
    S = np.zeros(shape=(len(t), n_vars))
    for i in range(n_vars):
        alpha_i = alpha[i] if isinstance(alpha, (tuple, list, np.ndarray)) else alpha
        beta_i = beta[i] if isinstance(beta, (tuple, list, np.ndarray)) else beta
        gamma_i = gamma[i] if isinstance(gamma, (tuple, list, np.ndarray)) else gamma
        tau, alpha_vec, u0_vec, s0_vec = vectorize(t, t_[i], alpha_i, beta_i, gamma_i, alpha_=alpha_, u0=0, s0=0)
        beta_i
        if noise_model is 'gillespie':
            U[:, i], S[:, i] = simulate_gillespie(alpha_vec, beta, gamma)
        else:
            U[:, i], S[:, i] = simulate_dynamics(tau, alpha_vec, beta_i, gamma_i,
                                                 u0_vec, s0_vec, noise_model, noise_level[i])

    obs = {'true_t': t.round(2)}
    var = {'true_t_': t_[:n_vars],
           'true_alpha': np.ones(n_vars) * alpha,
           'true_beta': np.ones(n_vars) * beta,
           'true_gamma': np.ones(n_vars) * gamma,
           'true_scaling': np.ones(n_vars)}
    layers = {'unspliced': U, 'spliced': S}

    return AnnData(S, obs, var, layers=layers)
576499
def rmv_normal(mu, tau, size=1):
    

    sig = np.linalg.cholesky(tau)
    mu_size = np.shape(mu)

    if size == 1:
        out = np.random.normal(size=mu_size)
        try:
            flib.dtrsm_wrap(sig, out, 'L', 'T', 'L', 1.)
        except:
            out = np.linalg.solve(sig, out)
        out += mu
        return out
    else:
        if not hasattr(size, '__iter__'):
            size = (size,)
        tot_size = np.prod(size)
        out = np.random.normal(size=(tot_size,) + mu_size)
        for i in xrange(tot_size):
            try:
                flib.dtrsm_wrap(sig, out[i, :], 'L', 'T', 'L', 1.)
            except:
                out[i, :] = np.linalg.solve(sig, out[i,:])
            out[i, :] += mu
        return out.reshape(size + mu_size)
53586
def _mean_square_error(y, y_pred, w):
    """"""
    return np.average(((y_pred - y) ** 2), weights=w)