Random multivariate hypergeometric variates.

    Parameters:
      - `n` : Number of draws.
      - `m` : Number of items in each categoy.
547818
def rmultivariate_hypergeometric(n, m, size=None):
    

    N = len(m)
    urn = np.repeat(np.arange(N), m)

    if size:
        draw = np.array([[urn[i] for i in np.random.permutation(len(urn))[:n]]
                         for j in range(size)])

        r = [[np.sum(draw[j] == i) for i in range(len(m))]
             for j in range(size)]
    else:
        draw = np.array([urn[i] for i in np.random.permutation(len(urn))[:n]])

        r = [np.sum(draw == i) for i in range(len(m))]
    return np.asarray(r)
1228456
def fill_noncontinous_variables(self, samples):
        
        init_points_count = samples.shape[0]

        for (idx, var) in enumerate(self.space.space_expanded):
            if isinstance(var, DiscreteVariable) or isinstance(var, CategoricalVariable) :
                sample_var = np.atleast_2d(np.random.choice(var.domain, init_points_count))
                samples[:,idx] = sample_var.flatten()

            # sample in the case of bandit variables
            elif isinstance(var, BanditVariable):
                # Bandit variable is represented by a several adjacent columns in the samples array
                idx_samples = np.random.randint(var.domain.shape[0], size=init_points_count)
                bandit_idx = np.arange(idx, idx + var.domain.shape[1])
                samples[:, bandit_idx] = var.domain[idx_samples,:]
863231
public static Matrix3D random(int M, int N, int K) {
        Matrix3D A = new Matrix3D(M, N, K);
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
            	for (int l = 0; l < K; l++)
            		A.data[i][j][l] = Math.random();
        return A;
    }
669111
def MeanVarEstEQ(y, x, covariates, tol=1e-8):
    

    pcount = covariates.shape[0] + 2
    N = y.shape[0]
    beta_count = pcount * 2
    X = numpy.ones((N, pcount))
    X[:, 1] = x

    aprime = [1, 0]
    idx = 2
    for cov in covariates:
        aprime.append(-numpy.mean(cov)/numpy.std(cov))
        X[:, idx] = cov
        idx += 1
    aprime = numpy.matrix(aprime)

    # http://stackoverflow.com/questions/7267226/range-for-floats
    def frange(x, y, jump):
        while x < y:
            yield x
            x += jump
        yield 0

    class PhiReturn(object):
        def __init__(self, phi, dtheta):
            self.phi = phi
            self.dtheta = dtheta

    def dot_diag(x, y):
        """This is a conveniance function to perform the following \
        in a more efficient manner:

        :param x: arr1
        :param y: arr2

        x.dot(numpy.diag(y))

        y must be a single dimensioned array

        """
        if len(y.shape) != 1:
            print >> sys.stderr, "You can't pass an array of shape %s to dot_diag" % (y.shape)
            sys.exit(1)

        result = numpy.empty(x.shape)
        for i in range(0, x.shape[0]):
            result[i] = x[i] * y

        return result

    def diag_dot(x, y):
        """diagonal multiply
        """
        result = numpy.empty(y.shape)
        for i in range(0, y.shape[0]):
            result[i] = y[i] * x[i]
        return result

    def Phi(theta):
        MM = y - numpy.dot(X, theta[0:pcount])
        SS = numpy.exp(numpy.dot(-X, theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        Phi = numpy.hstack((DD1.dot(X), (DD2-0.5).dot(X)))

        tX = X.transpose()

        t1 = numpy.empty(tX.shape)
        t2 = numpy.empty(tX.shape)
        t3 = numpy.empty(tX.shape)
        for i in range(0,tX.shape[0]):
            t1[i] = tX[i]*-SS
            t2[i] = tX[i]*-DD1
            t3[i] = tX[i]*-DD2
        t1 = t1.dot(X)
        t2 = t2.dot(X)
        t3 = t3.dot(X)

        dtheta = numpy.hstack((numpy.vstack((t1, t2)), numpy.vstack((t2, t3))))
        return PhiReturn(Phi, dtheta)



    class MvSolveReturn(object):
        def __init__(self, theta, dtheta):
            self.theta = theta
            self.dtheta = dtheta

    def MVsolve(theta_new):
        solution_found = False

        mvsolve_iterations = 0

        while not solution_found:
            theta_old = theta_new.copy()
            tmp = Phi(theta_old)
            #print "ITR", mvsolve_iterations, numpy.mean(theta_new),numpy.sum(numpy.absolute(theta_new - theta_old))
            theta_new = theta_old - scipy.linalg.solve(tmp.dtheta, tmp.phi)
            mvsolve_iterations += 1

            if (numpy.sum(numpy.absolute(theta_new - theta_old)) < tol):
                tmp = Phi(theta_new)
                solution_found = True
            else:
                if mvsolve_iterations > 25000:
                    #print >> sys.stderr, mvsolve_iterations, "failures"
                    raise UnsolvedLocus("")

        return MvSolveReturn(theta_new, tmp.dtheta), mvsolve_iterations

    def MVcalcB(theta):
        MM = y - X.dot(theta[0:pcount])
        SS = numpy.exp(-X.dot(theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        AA = numpy.hstack((diag_dot(DD1, X), diag_dot(DD2-0.5,X)))

        return numpy.transpose(AA).dot(AA)/ N


    mod = None
    itr = 0
    total_iterations = 0
    for i in frange(0.00, 1.0, 0.05):
        theta=numpy.empty((beta_count))
        theta[:] = i

        try:
            mod, iterations = MVsolve(theta)
            total_iterations += iterations
            if i > 0.05:
                print >> sys.stderr, "Completed: ", total_iterations, itr
            break
        except exceptions.ValueError as e:
            pass
        except numpy.linalg.linalg.LinAlgError as e:
            pass
        except Exception as inst:
            #print type(inst)
            pass
        itr += 1

    if not mod:
        raise UnsolvedLocus("")

    try:
        ainv = scipy.linalg.inv(mod.dtheta) * N
    except:
        raise ValueError("Singular Matrix Encountered")
    B = MVcalcB(mod.theta)
    V = ainv.dot(B).dot(ainv.transpose())



    # Focus on the two parameters of interest
    theta2 = numpy.array([mod.theta[1], mod.theta[pcount+1]])
    V2 = V[1:beta_count:pcount,1:beta_count:pcount]
    pvalt = 1 - scipy.stats.chi2.cdf(theta2.dot(scipy.linalg.inv(V2)).dot(theta2) * N, 2)


    ## From Chun's updated code:
    theta= mod.theta
    se = numpy.sqrt(numpy.diag(V)/N)
    pval = 2*scipy.stats.norm.cdf(-numpy.absolute(mod.theta/se))
    return pvalt, theta, pval, se, V/N
1257331
def _rand_sparse(m, n, density, format='csr'):
    """"""
    nnz = max(min(int(m*n*density), m*n), 0)

    row = np.random.randint(low=0, high=m-1, size=nnz)
    col = np.random.randint(low=0, high=n-1, size=nnz)
    data = np.ones(nnz, dtype=float)

    # duplicate (i,j) entries will be summed together
    return sp.sparse.csr_matrix((data, (row, col)), shape=(m, n))
473511
def multivariate_hypergeometric_expval(n, m):
    
    m = np.asarray(m, float)
    return n * (m / m.sum())
220928
def gaussian_polygons(points, n=10):
    
    gdf = gpd.GeoDataFrame(data={'cluster_number': classify_clusters(points, n=n)}, geometry=points)
    polygons = []
    for i in range(n):
        sel_points = gdf[gdf['cluster_number'] == i].geometry
        polygons.append(shapely.geometry.MultiPoint([(p.x, p.y) for p in sel_points]).convex_hull)
    polygons = [p for p in polygons if
                (not isinstance(p, shapely.geometry.Point)) and (not isinstance(p, shapely.geometry.LineString))]
    return gpd.GeoSeries(polygons)
1141526
def get_2D_samples_gauss(n, m, sigma, random_state=None):
    
    return make_2D_samples_gauss(n, m, sigma, random_state=None)
469101
def get_samples(self, n):
        

        normalized_w = self.weights / np.sum(self.weights)
        get_rand_index = st.rv_discrete(values=(range(self.N),
                                        normalized_w)).rvs(size=n)
        samples = np.zeros(n)
        k = 0
        j = 0
        while (k < n):
            i = get_rand_index[j]
            j = j + 1
            if (j == n):
                get_rand_index = st.rv_discrete(values=(range(self.N),
                                                normalized_w)).rvs(size=n)
                j = 0
            v = np.random.normal(loc=self.points[i], scale=self.sigma[i])
            if (v > self.max_limit or v < self.min_limit):
                continue
            else:
                samples[k] = v
                k = k + 1
                if (k == n):
                    break
        return samples
957636
def generate_covariates(n, d, n_binary=0, p=0.5):
    
    # pylint: disable=chained-comparison
    assert n_binary >= 0 and n_binary <= d, "binary must be between 0 and d"
    covariates = np.zeros((n, d + 1))
    covariates[:, : d - n_binary] = random.exponential(1, size=(n, d - n_binary))
    covariates[:, d - n_binary : -1] = random.binomial(1, p, size=(n, n_binary))
    covariates[:, -1] = np.ones(n)
    return covariates
52361
def random_samples(self, n):
        
        samples = self.samples
        inds = rand.randint(len(samples),size=int(n))

        newsamples = samples.iloc[inds]
        newsamples.reset_index(inplace=True)
        return newsamples
1182854
def _kmeans_init(X, n_clusters, method='balanced', rng=None):
    
    n_samples = X.shape[0]
    if rng is None:
        cent_idx = np.random.choice(n_samples, replace=False, size=n_clusters)
    else:
        #print('Generate random centers using RNG')
        cent_idx = rng.choice(n_samples, replace=False, size=n_clusters)
    
    centers = X[cent_idx,:]
    mean_X = np.mean(X, axis=0)
    
    if method == 'balanced':
        centers[n_clusters-1] = n_clusters*mean_X - np.sum(centers[:(n_clusters-1)], axis=0)
    
    return (centers)
541750
def random_draw(self, size=None):
        
        return scipy.asarray([p.rvs(size=size) for p in self.univariate_priors])
355928
def _simulate_coef_from_bootstraps(
            self, n_draws, coef_bootstraps, cov_bootstraps):
        """"""
        # Sample indices uniformly from {0, ..., n_bootstraps - 1}
        # (Wood pg. 199 step 6)
        random_bootstrap_indices = np.random.choice(
            np.arange(len(coef_bootstraps)), size=n_draws, replace=True)

        # Simulate `n_draws` many random coefficient vectors from a
        # multivariate normal distribution with mean and covariance given by
        # the bootstrap samples (indexed by `random_bootstrap_indices`) of
        # `coef_bootstraps` and `cov_bootstraps`. Because it's faster to draw
        # many samples from a certain distribution all at once, we make a dict
        # mapping bootstrap indices to draw indices and use the `size`
        # parameter of `np.random.multivariate_normal` to sample the draws
        # needed from that bootstrap sample all at once.
        bootstrap_index_to_draw_indices = defaultdict(list)
        for draw_index, bootstrap_index in enumerate(random_bootstrap_indices):
            bootstrap_index_to_draw_indices[bootstrap_index].append(draw_index)

        coef_draws = np.empty((n_draws, len(self.coef_)))

        for bootstrap, draw_indices in bootstrap_index_to_draw_indices.items():
            coef_draws[draw_indices] = np.random.multivariate_normal(
                coef_bootstraps[bootstrap], cov_bootstraps[bootstrap],
                size=len(draw_indices))

        return coef_draws
1226916
def Sample(self, n):
        
        size = n,
        return numpy.random.beta(self.alpha, self.beta, size)
936302
def make_data(n):
    
    p,r,d,w = {},{},{},{}

    J = range(1,n+1)

    for j in J:
        p[j] = random.randint(1,4)
        w[j] = random.randint(1,3)

    T = sum(p)
    for j in J:
        r[j] = random.randint(0,5)
        d[j] = r[j] + random.randint(0,5)

    return J,p,r,d,w
594047
def random_pairs_with_replacement(n, shape, random_state=None):
    """"""

    if not isinstance(random_state, np.random.RandomState):
        random_state = np.random.RandomState(random_state)

    n_max = max_pairs(shape)

    if n_max <= 0:
        raise ValueError('n_max must be larger than 0')

    # make random pairs
    indices = random_state.randint(0, n_max, n)

    if len(shape) == 1:
        return _map_tril_1d_on_2d(indices, shape[0])
    else:
        return np.unravel_index(indices, shape)