Python numpy.random.poisson() Examples

def generate(shp=(30, 30, 20, 10), K=5, alpha=0.1, beta=0.1):
    """Generate a count tensor from the BPTF model.

    PARAMS:
    shp -- (tuple) shape of the generated count tensor
    K -- (int) number of latent components
    alpha -- (float) shape parameter of gamma prior over factors
    beta -- (float) rate parameter of gamma prior over factors

    RETURNS:
    Mu -- (np.ndarray) true Poisson rates
    Y -- (np.ndarray) generated count tensor
    """
    Theta_DK_M = [rn.gamma(alpha, 1./beta, size=(D, K)) for D in shp]
    Mu = parafac(Theta_DK_M)
    assert Mu.shape == shp
    Y = rn.poisson(Mu)
    return Mu, Y 

def corrupt(Y, p=0.05):
    """Corrupt a count tensor with anomalies.

    The corruption noise model is:

        corrupt(y) = y * g, where g ~ Gamma(10, 2)

    PARAMS:
    p -- (float) proportion of tensor entries to corrupt

    RETURNS:
    out -- (np.ndarray) corrupted count tensor
    mask -- (np.ndarray) boolean array, same shape as count tensor
                         True means that entry was corrupted.
    """
    out = Y.copy()
    mask = (rn.random(size=out.shape) < p).astype(bool)
    out[mask] = rn.poisson(out[mask] * rn.gamma(10., 2., size=out[mask].shape))
    return out, mask 

def rvs(self, size=None):
        return random.poisson(self.rate, size=size) 

def logpdf(self, x):
        return stats.poisson.logpmf(x, self.rate) 

def ppf(self, u):
        return stats.poisson.ppf(u, self.rate) 

def generate_spot(p):
    position = rand(p['dimension'])
    gene = random.choice(range(len(codebook)))
    barcode = array(codebook[gene])
    photons = [poisson(p['N_photons_per_flour'])*poisson(p['N_flour'])*b for b in barcode]
    return DataFrame({'position': [position], 'barcode': [barcode], 'photons': [photons], 'gene':gene})

# right now there is no jitter on positions of the spots, we might want to make it a vector
# EPY: END code

# EPY: START code 

def generate_events(self,start_date,end_date):
        '''
        Generate a sequence of events at the customers daily rates.  Each count for an event on a day is droing from
        a poisson distribution with the customers average rate.  If the number is greater than zero, that number of events
        are created as tuples of time stamps and the event index (which is the database type id).  The time of the
        event is randomly set to anything on the 24 hour range.
        :param start_date: datetime.date for start of simulation
        :param end_date: datetime.date for end of simulation
        :return: The total count of each event, the list of all of the event tuples
        '''

        delta = end_date - start_date

        events=[]
        counts=[0]*len(self.behave_per_day)
        for i in range(delta.days):
            the_date = start_date + timedelta(days=i)
            if the_date in Customer.date_multipliers:
                multiplier = Customer.date_multipliers[the_date]
            else:
                if the_date.weekday() >= 4:
                    multiplier = random.uniform(1.00,1.2)
                else:
                    multiplier = random.uniform(0.825,1.025)
                Customer.date_multipliers[the_date]=multiplier
            for event_idx,rate in  enumerate(self.behave_per_day):
                new_count= int(round(multiplier*random.poisson(rate)))
                counts[event_idx] += new_count
                for n in range(0,new_count):
                    event_time=datetime.combine(the_date,time(randrange(24),randrange(60),randrange(60)))
                    new_event=(event_time,event_idx)
                    events.append(new_event )

        self.events.extend(events)

        return counts 

def seed_links(self):
        # bias with non-random seed links
        self.bias_links()
        if self.link_count() < self.max_links_cached:
            num_words = max(1,npr.poisson(1.33)+1)  # mean of 1.33 words per search
            if num_words == 1:
                word = ' '.join(random.sample(self.words,num_words))
            else:
                if npr.uniform() < 0.5:
                    word = ' '.join(random.sample(self.words,num_words))
                else:      # quote the first two words together
                    word = ' '.join(['"{}"'.format(' '.join(random.sample(self.words, 2))),
                                     ' '.join(random.sample(self.words, num_words-2))])
            if self.debug: print(f'Seeding with search for \'{word}\'…')
            self.get_websearch(word) 

def __call__(self):
    return nr.poisson(lam=self.s, size=np.size(self.s)) - self.s 

def detect(Y, K=5, alpha=0.1, thresh=1e-5):
    """Detect anomalies using BPTF.

    This method fits BPTF to Y and obtains Mu, which is the model's
    reconstruction of Y (computed from the inferred latent factors).
    Anomalies are then all entries of Y whose probability given Mu
    is less than a given threshold.

        If P(y | mu) < thresh ==> y is  anomaly!

        Here P(y | mu) = Pois(y; mu), the PMF of the Poisson distribution.

    PARAMS:
    Y -- (np.ndarray) data count tensor
    K -- (int) number of latent components
    alpha -- (float) shape parameter of gamma prior over factors
    thresh -- (float) anomaly threshold (between 0 and 1).
    """
    bptf = BPTF(n_modes=Y.ndim,
                n_components=K,
                max_iter=100,
                tol=1e-4,
                smoothness=100,
                verbose=False,
                alpha=alpha,
                debug=False)
    bptf.fit(Y)
    Mu = bptf.reconstruct()

    return st.poisson.pmf(Y, Mu) < thresh 

def setup_class(self):
        from scitbx.array_family import flex

        spot = flex.double(flex.grid(11, 11))
        for j in range(11):
            for i in range(11):
                spot[j, i] = exp(-((j - 5) ** 2 + (i - 5) ** 2) / 2 ** 2)

        self.image = flex.double(flex.grid(2000, 2000))
        for n in range(200):
            x = randint(0, 2000)
            y = randint(0, 2000)
            for j in range(0, 11):
                for i in range(0, 11):
                    xx = x + i - 5
                    yy = y + j - 5
                    if xx >= 0 and yy >= 0 and xx < 2000 and yy < 2000:
                        self.image[yy, xx] = poisson(100 * spot[j, i])

        background = flex.double(list(poisson(5, 2000 * 2000)))
        background.reshape(flex.grid(2000, 2000))

        self.image += background

        # Create an image
        self.mask = flex.random_bool(2000 * 2000, 0.99)
        self.mask.reshape(flex.grid(2000, 2000))
        self.gain = flex.random_double(2000 * 2000) + 1.0
        self.gain.reshape(flex.grid(2000, 2000))
        self.size = (3, 3)
        self.min_count = 2 

def RobotStep(args):
  global t

  # Input from Sensors
  SL, SR = Env.GetSensors(x[t], y[t], w[t])

  RL_spikes = 0.0
  RR_spikes = 0.0
  for t2 in xrange(dt):

    # Deliver stimulus as a Poisson spike stream to the sensor neurons and
    # noisy base current to the motor neurons
    I = np.hstack([rn.poisson(SL*15, N1), rn.poisson(SR*15, N2),
                   5*rn.randn(N3), 5*rn.randn(N4)])

    # Update network
    net.setCurrent(I)
    fired = net.update()

    RL_spikes += np.sum(np.logical_and(fired > (N1+N2), fired < N1+N2+N3))
    RR_spikes += np.sum(fired > (N1+N2+N3))

    # Maintain record of membrane potential
    v[t2,:],_ = net.getState()

  # Output to motors
  # Calculate motor firing rates in Hz
  RL = 1.0*RL_spikes/(dt*N3)*1000.0
  RR = 1.0*RR_spikes/(dt*N4)*1000.0

  # Set wheel velocities (as fractions of Umax)
  UL = (Umin/Umax + RL/Rmax*(1 - Umin/Umax))
  UR = (Umin/Umax + RR/Rmax*(1 - Umin/Umax))

  # Update Environment
  x[t+1], y[t+1], w[t+1] = RobotUpdate(x[t], y[t], w[t], UL, UR,
                                       Umax, dt, xmax, ymax)

  ## PLOTTING
  for i in range(Ns):
    pl11[i].set_data(range(dt), v[:,i])
    pl12[i].set_data(range(dt), v[:,i+Ns])

  for i in range(Nm):
    pl21[i].set_data(range(dt), v[:,2*Ns+i])
    pl22[i].set_data(range(dt), v[:,2*Ns+Nm+i])

  ax2.scatter(x, y)
  manager1.canvas.draw()
  manager2.canvas.draw()

  t += 1

  if t == len(x)-1:
    print 'Terminating simulation'
    StopSimulation()

# Get the thing going