Python numpy.random.multinomial() Examples

def draw_links(self,n=1,log_sampling=False):
        """ Draw multiple random links. """
        urls = []
        domain_array = np.array([dmn for dmn in self.domain_links])
        domain_count = np.array([len(self.domain_links[domain_array[k]]) for k in range(domain_array.shape[0])])
        p = np.array([np.float(c) for c in domain_count])
        count_total = p.sum()
        if log_sampling:  # log-sampling [log(x+1)] to bias lower count domains
            p = np.fromiter((np.log1p(x) for x in p), dtype=p.dtype)
        if count_total > 0:
            p = p/p.sum()
            cnts = npr.multinomial(n, pvals=p)
            if n > 1:
                for k in range(cnts.shape[0]):
                    domain = domain_array[k]
                    cnt = min(cnts[k],domain_count[k])
                    for url in random.sample(self.domain_links[domain],cnt):
                        urls.append(url)
            else:
                k = int(np.nonzero(cnts)[0])
                domain = domain_array[k]
                url = random.sample(self.domain_links[domain],1)[0]
                urls.append(url)
        return urls 

def test_basic(self):
        random.multinomial(100, [0.2, 0.8]) 

def test_multinomial(self):
        def gen_random(state, out):
            out[...] = state.multinomial(10, [1/6.]*6, size=10000)
        self.check_function(gen_random, sz=(10000, 6))

# See Issue #4263 

def test_multinomial(self):
        np.random.seed(self.seed)
        actual = np.random.multinomial(20, [1/6.]*6, size=(3, 2))
        desired = np.array([[[4, 3, 5, 4, 2, 2],
                             [5, 2, 8, 2, 2, 1]],
                            [[3, 4, 3, 6, 0, 4],
                             [2, 1, 4, 3, 6, 4]],
                            [[4, 4, 2, 5, 2, 3],
                             [4, 3, 4, 2, 3, 4]]])
        assert_array_equal(actual, desired) 

def test_zero_probability(self):
        random.multinomial(100, [0.2, 0.8, 0.0, 0.0, 0.0]) 

def test_size(self):
        # gh-3173
        p = [0.5, 0.5]
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, [2, 2]).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, (2, 2)).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, np.array((2, 2))).shape,
                     (2, 2, 2))

        assert_raises(TypeError, np.random.multinomial, 1, p,
                      np.float(1)) 

def test_basic(self):
        random.multinomial(100, [0.2, 0.8]) 

def test_size(self):
        # gh-3173
        p = [0.5, 0.5]
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, [2, 2]).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, (2, 2)).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, np.array((2, 2))).shape,
                     (2, 2, 2))

        assert_raises(TypeError, np.random.multinomial, 1, p,
                      float(1)) 

def test_zero_probability(self):
        random.multinomial(100, [0.2, 0.8, 0.0, 0.0, 0.0]) 

def test_basic(self):
        random.multinomial(100, [0.2, 0.8]) 

def test_multinomial(self):
        def gen_random(state, out):
            out[...] = state.multinomial(10, [1/6.]*6, size=10000)
        self.check_function(gen_random, sz=(10000, 6))

# See Issue #4263 

def test_multinomial(self):
        np.random.seed(self.seed)
        actual = np.random.multinomial(20, [1/6.]*6, size=(3, 2))
        desired = np.array([[[4, 3, 5, 4, 2, 2],
                             [5, 2, 8, 2, 2, 1]],
                            [[3, 4, 3, 6, 0, 4],
                             [2, 1, 4, 3, 6, 4]],
                            [[4, 4, 2, 5, 2, 3],
                             [4, 3, 4, 2, 3, 4]]])
        assert_array_equal(actual, desired) 

def test_multinomial(self):
        def gen_random(state, out):
            out[...] = state.multinomial(10, [1/6.]*6, size=10000)
        self.check_function(gen_random, sz=(10000, 6))

# See Issue #4263 

def test_multinomial(self):
        np.random.seed(self.seed)
        actual = np.random.multinomial(20, [1/6.]*6, size=(3, 2))
        desired = np.array([[[4, 3, 5, 4, 2, 2],
                             [5, 2, 8, 2, 2, 1]],
                            [[3, 4, 3, 6, 0, 4],
                             [2, 1, 4, 3, 6, 4]],
                            [[4, 4, 2, 5, 2, 3],
                             [4, 3, 4, 2, 3, 4]]])
        assert_array_equal(actual, desired) 

def test_size(self):
        # gh-3173
        p = [0.5, 0.5]
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, [2, 2]).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, (2, 2)).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, np.array((2, 2))).shape,
                     (2, 2, 2))

        assert_raises(TypeError, np.random.multinomial, 1, p,
                      np.float(1)) 

def test_zero_probability(self):
        random.multinomial(100, [0.2, 0.8, 0.0, 0.0, 0.0]) 

def sample_random_read(gene, true_psi, read_len, overhang_len):
    """
    Sample a random read (not taking into account overhang) from the
    given set of exons and the true Psi value.

    Note that if we're given a gene that has only two isoforms, the 'align' function of
    gene will return a read summary in the form of (NI, NE, NB) rather than an alignment to
    the two isoforms (which is a pair (0/1, 0/1)).
    """
    iso_lens = [iso.len for iso in gene.isoforms]
    num_positions = array([(l - read_len + 1) for l in iso_lens])
    # probability of sampling a particular position from an isoform -- assume uniform for now
    iso_probs = [1/float(n) for n in num_positions]
    psi_frag_denom = sum(num_positions * array(true_psi))
    psi_frags = [(num_pos * curr_psi)/psi_frag_denom for num_pos, curr_psi \
                 in zip(num_positions, true_psi)]
    # Choose isoform to sample read from
    chosen_iso = list(multinomial(1, psi_frags)).index(1)
    isoform_position_prob = ones(num_positions[chosen_iso]) * iso_probs[chosen_iso]
    sampled_read_start = list(multinomial(1, isoform_position_prob)).index(1)
    sampled_read_end = sampled_read_start + read_len - 1
#    seq = gene.isoforms[chosen_iso].seq[sampled_read_start:sampled_read_end]
#    alignment, category = gene.align(seq, overhang=overhang_len)
    ##
    ## Trying out new alignment method
    ##
    # convert coordinates to genomic
    genomic_read_start, genomic_read_end = \
			gene.isoforms[chosen_iso].isoform_coords_to_genomic(sampled_read_start,
									    sampled_read_end)
    alignment, category = gene.align_read(genomic_read_start, genomic_read_end, overhang=overhang_len)
    return (tuple(alignment), [sampled_read_start, sampled_read_end], category, chosen_iso) 

def draw_domain(self,log_sampling=False):
        """ Draw a single, random domain. """
        domain = None
        domain_array = np.array([dmn for dmn in self.domain_links])
        domain_count = np.array([len(self.domain_links[domain_array[k]]) for k in range(domain_array.shape[0])])
        p = np.array([np.float(c) for c in domain_count])
        count_total = p.sum()
        if log_sampling:  # log-sampling [log(x+1)] to bias lower count domains
            p = np.fromiter((np.log1p(x) for x in p), dtype=p.dtype)
        if count_total > 0:
            p = p/p.sum()
            cnts = npr.multinomial(1, pvals=p)
            k = int(np.nonzero(cnts)[0])
            domain = domain_array[k]
        return domain 

def test_multinomial(self):
        def gen_random(state, out):
            out[...] = state.multinomial(10, [1/6.]*6, size=10000)
        self.check_function(gen_random, sz=(10000,6)) 

def test_multinomial(self):
        np.random.seed(self.seed)
        actual = np.random.multinomial(20, [1/6.]*6, size=(3, 2))
        desired = np.array([[[4, 3, 5, 4, 2, 2],
                             [5, 2, 8, 2, 2, 1]],
                            [[3, 4, 3, 6, 0, 4],
                             [2, 1, 4, 3, 6, 4]],
                            [[4, 4, 2, 5, 2, 3],
                             [4, 3, 4, 2, 3, 4]]])
        np.testing.assert_array_equal(actual, desired) 

def test_size(self):
        # gh-3173
        p = [0.5, 0.5]
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, [2, 2]).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, (2, 2)).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, np.array((2, 2))).shape,
                     (2, 2, 2))

        assert_raises(TypeError, np.random.multinomial, 1, p,
                      np.float(1)) 

def test_zero_probability(self):
        random.multinomial(100, [0.2, 0.8, 0.0, 0.0, 0.0]) 

def test_basic(self):
        random.multinomial(100, [0.2, 0.8]) 

def test_multinomial(self):
        def gen_random(state, out):
            out[...] = state.multinomial(10, [1/6.]*6, size=10000)
        self.check_function(gen_random, sz=(10000, 6))

# See Issue #4263 

def test_multinomial(self):
        np.random.seed(self.seed)
        actual = np.random.multinomial(20, [1/6.]*6, size=(3, 2))
        desired = np.array([[[4, 3, 5, 4, 2, 2],
                             [5, 2, 8, 2, 2, 1]],
                            [[3, 4, 3, 6, 0, 4],
                             [2, 1, 4, 3, 6, 4]],
                            [[4, 4, 2, 5, 2, 3],
                             [4, 3, 4, 2, 3, 4]]])
        assert_array_equal(actual, desired) 

def test_size(self):
        # gh-3173
        p = [0.5, 0.5]
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, np.uint32(1)).shape, (1, 2))
        assert_equal(np.random.multinomial(1, p, [2, 2]).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, (2, 2)).shape, (2, 2, 2))
        assert_equal(np.random.multinomial(1, p, np.array((2, 2))).shape,
                     (2, 2, 2))

        assert_raises(TypeError, np.random.multinomial, 1, p,
                      np.float(1)) 

def test_zero_probability(self):
        random.multinomial(100, [0.2, 0.8, 0.0, 0.0, 0.0]) 

def test_basic(self):
        random.multinomial(100, [0.2, 0.8]) 

def read_counts_to_read_list(ni, ne, nb):
    """
    Convert a set of read counts for a two-isoform gene (NI, NE, NB) to a list of reads.
    """
    reads = []
    reads.extend(ni * [[1, 0]])
    reads.extend(ne * [[0, 1]])
    reads.extend(nb * [[1, 1]])
    return array(reads)

# def sample_random_read(gene, true_psi, read_len, overhang_len):
#     """
#     Sample a random read (not taking into account overhang) from the
#     given set of exons and the given true Psi value.
#     """
#     iso1_len = gene.isoforms[0]['len']
#     iso2_len = gene.isoforms[1]['len']
#     num_inc = 0
#     num_exc = 0
#     num_both = 0
#     num_positions_iso1 = iso1_len - read_len + 1
#     num_positions_iso2 = iso2_len - read_len + 1
#     p1 = 1/float(num_positions_iso1)
#     p2 = 1/float(num_positions_iso2)
#     psi_frag = (num_positions_iso1*true_psi)/((num_positions_iso1*true_psi + num_positions_iso2*(1-true_psi)))
#     # Choose isoform to sample read from
#     if rand() < psi_frag:
#         isoform_position_prob = ones(num_positions_iso1)*p1
#         sampled_read_start = list(multinomial(1, isoform_position_prob)).index(1)
#         sampled_read_end = sampled_read_start + read_len
#         seq = gene.isoforms[0]['seq'][sampled_read_start:sampled_read_end]
#         [n1, n2, nb], category = gene.align_two_isoforms(seq, overhang=overhang_len)
#         return [[n1, n2, nb], [sampled_read_start, sampled_read_end], category]
#     else:
#         isoform_position_prob = ones(num_positions_iso2)*p2
#         sampled_read_start = list(multinomial(1, isoform_position_prob)).index(1)
#         sampled_read_end = sampled_read_start + read_len
#         seq = gene.isoforms[1]['seq'][sampled_read_start:sampled_read_end]
#         [n1, n2, nb], category = gene.align_two_isoforms(seq, overhang=overhang_len)
#         return [[n1, n2, nb], [sampled_read_start, sampled_read_end], category] 

def test_basic(self):
        random.multinomial(100, [0.2, 0.8])