我们从Python开源项目中，提取了以下13个代码示例，用于说明如何使用numpy.random.gamma()。
def sample_dirichlet(alpha, n_samples=1): """ Sample points from a dirichlet distribution with parameter alpha. @param alpha: alpha parameter of a dirichlet distribution @type alpha: array """ from numpy import array, sum, transpose, ones from numpy.random import gamma alpha = array(alpha, ndmin=1) X = gamma(alpha, ones(len(alpha)), [n_samples, len(alpha)]) return transpose(transpose(X) / sum(X, -1))
def gen_inv_gaussian(a, b, p, burnin=10): """ Sampler based on Gibbs sampling. Assumes scalar p. """ from numpy.random import gamma from numpy import sqrt s = a * 0. + 1. if p < 0: a, b = b, a for i in range(burnin): l = b + 2 * s m = sqrt(l / a) x = inv_gaussian(m, l, shape=m.shape) s = gamma(abs(p) + 0.5, x) if p >= 0: return x else: return 1 / x
def optimize_hyper_hdp(self): # Optimize \alpha_0 m_dot = self.msampler.m_dotk.sum() alpha_0 = self.zsampler.alpha_0 n_jdot = np.array(self.zsampler.data_dims) #norm = np.linalg.norm(n_jdot/alpha_0) #u_j = binomial(1, n_jdot/(norm* alpha_0)) u_j = binomial(1, n_jdot/(n_jdot + alpha_0)) v_j = beta(alpha_0 + 1, n_jdot) new_alpha0 = gamma(self.a_alpha + m_dot - u_j.sum(), 1/( self.b_alpha - np.log(v_j).sum()), size=5).mean() self.zsampler.alpha_0 = new_alpha0 # Optimize \gamma K = self.zsampler._K gmma = self.betasampler.gmma #norm = np.linalg.norm(m_dot/gmma) #u = binomial(1, m_dot / (norm*gmma)) u = binomial(1, m_dot / (m_dot + gmma)) v = beta(gmma + 1, m_dot) new_gmma = gamma(self.a_gmma + K -1 + u, 1/(self.b_gmma - np.log(v)), size=5).mean() self.betasampler.gmma = new_gmma print('alpha a, b: %s, %s ' % (self.a_alpha + m_dot - u_j.sum(), 1/( self.b_alpha - np.log(v_j).sum()))) print( 'hyper sample: alpha_0: %s gamma: %s' % (new_alpha0, new_gmma)) return
def truncated_gamma(shape=None, alpha=1., beta=1., x_min=None, x_max=None): """ Generate random variates from a lower-and upper-bounded gamma distribution. @param shape: shape of the random sample @param alpha: shape parameter (alpha > 0.) @param beta: scale parameter (beta >= 0.) @param x_min: lower bound of variate @param x_max: upper bound of variate @return: random variates of lower-bounded gamma distribution """ from scipy.special import gammainc, gammaincinv from numpy.random import gamma from numpy import inf if x_min is None and x_max is None: return gamma(alpha, 1 / beta, shape) elif x_min is None: x_min = 0. elif x_max is None: x_max = inf x_min = max(0., x_min) x_max = min(1e300, x_max) a = gammainc(alpha, beta * x_min) b = gammainc(alpha, beta * x_max) return probability_transform(shape, lambda x, alpha=alpha: gammaincinv(alpha, x), a, b) / beta
def _create_prices(t): last_average = 100 if t==0 else source.data['average'][-1] returns = asarray(lognormal(mean.value, stddev.value, 1)) average = last_average * cumprod(returns) high = average * exp(abs(gamma(1, 0.03, size=1))) low = average / exp(abs(gamma(1, 0.03, size=1))) delta = high - low open = low + delta * uniform(0.05, 0.95, size=1) close = low + delta * uniform(0.05, 0.95, size=1) return open[0], high[0], low[0], close[0], average[0]
def _create_prices(t): global last_average returns = asarray(lognormal(mean, stddev, 1)) average = last_average * cumprod(returns) last_average = average high = average * exp(abs(gamma(1, 0.03, size=1))) low = average / exp(abs(gamma(1, 0.03, size=1))) delta = high - low open = low + delta * uniform(0.05, 0.95, size=1) close = low + delta * uniform(0.05, 0.95, size=1) return open[0], high[0], low[0], close[0], average[0]
def fit(self): #self.fdebug() print( '__init__ Entropy: %f' % self.entropy()) for _it in range(self.iterations): self._iteration = _it ### Sampling begin_it = datetime.now() self.sample() self.time_it = (datetime.now() - begin_it).total_seconds() / 60 if _it >= self.iterations - self.burnin: if _it % self.thinning == 0: self.samples.append([self._theta, self._phi]) self.compute_measures() print('.', end='') lgg.info('Iteration %d, Entropy: %f \t\t K=%d alpha: %f gamma: %f' % (_it, self._entropy, self._K,self._alpha, self._gmma)) if self.write: self.write_it_step(self) if _it > 0 and _it % self.snapshot_freq == 0: self.save(silent=True) sys.stdout.flush() ### Clean Things if not self.samples: self.samples.append([self._theta, self._phi]) self._reduce_latent() self.samples = None # free space return
def __init__(self, sampler, data_t=None, **kwargs): self.burnin = kwargs.get('burnin', 0.05) # Ratio of iteration self.thinning = kwargs.get('thinning', 1) self.comm = dict() # Empty dict to store communities and blockmodel structure self.data_t = data_t self._csv_typo = '# it it_time entropy_train entropy_test K alpha gamma alpha_mean delta_mean alpha_var delta_var' self.fmt = '%d %.4f %.8f %.8f %d %.8f %.8f %.4f %.4f %.4f %.4f' #self.fmt = '%s %s %s %s %s %s %s %s %s %s %s' GibbsSampler.__init__(self, sampler, **kwargs)
def optimize_hyper_hdp(self): # Optimize \alpha_0 m_dot = self.msampler.m_dotk.sum() alpha_0 = self.zsampler.alpha_0 n_jdot = np.array(self.zsampler.data_dims, dtype=float) # @debug add row count + line count for masked ! #p = np.power(n_jdot / alpha_0, np.arange(n_jdot.shape[0])) #norm = np.linalg.norm(p) #u_j = binomial(1, p/norm) u_j = binomial(1, alpha_0/(n_jdot + alpha_0)) #u_j = binomial(1, n_jdot/(n_jdot + alpha_0)) try: v_j = beta(alpha_0 + 1, n_jdot) except: #n_jdot[n_jdot == 0] = np.finfo(float).eps lgg.warning('Unable to optimize MMSB parameters, possible empty sequence...') return shape_a = self.a_alpha + m_dot - u_j.sum() if shape_a <= 0: lgg.warning('Unable to optimize MMSB parameters, possible empty sequence...') return new_alpha0 = gamma(shape_a, 1/( self.b_alpha - np.log(v_j).sum()), size=3).mean() self.zsampler.alpha_0 = new_alpha0 # Optimize \gamma K = self.zsampler._K #m_dot = self.msampler.m_dotk gmma = self.betasampler.gmma #p = np.power(m_dot / gmma, np.arange(m_dot.shape[0])) #norm = np.linalg.norm(p) #u = binomial(1, p/norm) u = binomial(1, gmma / (m_dot + gmma)) #u = binomial(1, m_dot / (m_dot + gmma)) v = beta(gmma + 1, m_dot) new_gmma = gamma(self.a_gmma + K -1 + u, 1/(self.b_gmma - np.log(v)), size=3).mean() self.betasampler.gmma = new_gmma #print 'm_dot %d, alpha a, b: %s, %s ' % (m_dot, self.a_alpha + m_dot - u_j.sum(), 1/( self.b_alpha - np.log(v_j).sum())) #print 'gamma a, b: %s, %s ' % (self.a_gmma + K -1 + u, 1/(self.b_gmma - np.log(v))) lgg.debug('hyper sample: alpha_0: %s gamma: %s' % (new_alpha0, new_gmma)) return
def __init__(self, expe, frontend): self.comm = dict() # Empty dict to store communities and blockmodel structure #self._csv_typo = '# it it_time likelihood likelihood_t K alpha gamma alpha_mean delta_mean alpha_var delta_var' self.fmt = '%d %.4f %.8f %.8f %d %.8f %.8f %.4f %.4f %.4f %.4f' #self.fmt = '%s %s %s %s %s %s %s %s %s %s %s' GibbsSampler.__init__(self, expe, frontend)
def __init__(self, n=None, method='stub', **distribution): """ Possible arguments for the distribution are: - network_type: specify the type of network that should be constructed (THIS IS MANDATORY). It can either be the name of a distribution or of a certain network type. ['l_partition', 'poisson', 'normal', 'binomial', 'exponential', 'geometric', 'gamma', 'power', 'weibull'] For specific parameters of the distributions, see: http://docs.scipy.org/doc/numpy/reference/routines.random.html - method: The probabilistic framework after which the network will be constructed. - distribution specific arguments. Check out the description of the specific numpy function. Or just give the argument network_type and look at what the error tells you. see self._create_graph for more information """ _Graph.__init__(self) self.is_static = True self._rewiring_attempts = 100000 self._stub_attempts = 100000 self.permitted_types = allowed_dists + ["l_partition", 'full'] self.is_directed = False # to do: pass usefull info in here. self._info = {} #for now only undirected networks if n is not None: self.n = n if method in ['proba', 'stub']: self.method = method else: raise ValueError(method + ' is not a permitted method! Chose either "proba" or "stub"') try: self.nw_name = distribution.pop('network_type') empty_graph = False except KeyError: self.nn = [] self._convert_to_array() empty_graph = True #create an empty graph if network_type is not given if not empty_graph: if self.nw_name not in self.permitted_types: raise ValueError( "The specified network type \"%s\" is not permitted. \ Please chose from " % self.nw_name + '[' + ', '.join(self.permitted_types) + ']') self.distribution = distribution self._create_graph(**self.distribution)
