Python scipy.stats 模块，rv_discrete() 实例源码

def sampledEF(dv,sv,N):
    #rv_discrete only handles integers... making a dict to translate the range into the central bin value
    vals = np.arange(0.5,599.5,1.0)
    valDict = {i:vals[i] for i in range(len(vals))}
    sv_rand = stats.rv_discrete(name='sv_rand',values=(range(len(vals)),sv))  #creating a RNG using the PDF as a weighting function
    dv_rand = stats.rv_discrete(name='dv_rand',values=(range(len(vals)),dv))
    randSVs = sv_rand.rvs(size=N) #sampling 10k systolic volumes 
    randDVs = dv_rand.rvs(size=N)
    randSVs = [valDict[x] for x in randSVs] # converting integer to central bin values
    randDVs = [valDict[x] for x in randDVs]
    randEFs= [(float(x)-float(y))/float(x) for x,y in zip(randDVs,randSVs)]  # array of EFs generated from the SV and DV samples
    (binCont,bins)=np.histogram(randEFs,bins=101,range=(0,1.0))   # histogramming the EFs
    EF_pdf = np.array(binCont,dtype=np.float)  
    EF_bin_cent = [(bins[i]+bins[i+1])/2.0 for i in range(len(bins)-1)]  # central values of EF bins
    EF_pdf = EF_pdf/np.sum(EF_pdf)    # normalizing the EF PDF
    # Calculating Expectation Value of EF  EV_EF
    mean_EF = 0.0
    for i in range(len(EF_bin_cent)):
        mean_EF += EF_bin_cent[i]*EF_pdf[i]
    return mean_EF

def act(self, action):
        try:
            assert not self.terminal()
        except AssertionError as e:
            e.args += ('Further action not permitted: terminal state ' +
                       ' reached. Episode is over.',)
            raise

        probs = self.T[action][self.state, :]
        pmf = stats.rv_discrete(name='pmf',
                                values=(self.states, probs))
        successor_state = pmf.rvs()
        self.state = successor_state

        r = self.rewards[successor_state]
        return r

def act(self, action):
        try:
            assert not self.terminal()
        except AssertionError as e:
            e.args += ('Further action not permitted: terminal state ' +
                       ' reached. Episode is over.',)
            raise

        probs = self.T[action][self.state, :]
        pmf = stats.rv_discrete(name='pmf',
                                values=(self.states, probs))
        successor_state = pmf.rvs()
        self.state = successor_state

        r = self.rewards[successor_state]
        return r

def plot_categorical(X, output_dir, parameter_name):
    try:
        os.makedirs(output_dir)
    except FileExistsError:
        pass

    X_prime = collections.OrderedDict()
    for value in X:
        if value not in X_prime:
            X_prime[value] = 0
        X_prime[value] += (1.0 / len(X))
    distrib = rv_discrete(values=(list(range(len(X_prime))), list(X_prime.values())))

    fig, ax = plt.subplots()
    # TODO: resampled from dist, but will do.
    ax.hist(distrib.rvs(size=100), range=(0, len(X_prime)))
    ax.legend(loc='upper left')

    plt.savefig(output_dir + parameter_name + '.png', bbox_inches='tight')
    plt.close()

def get_samples(self, n):
        """Sample the GMM distribution.

        Arguments
        ---------
        n : int
            Number of samples needed

        Returns
        -------
        1D array
            Samples from the distribution
        """

        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

def from_discrete_dist(dist, sample_num=1):
    return rv_discrete(values=(np.arange(len(dist)), dist)).rvs(size=sample_num)

def from_discrete_dist(dist, sample_num=1):
    return rv_discrete(values=(np.arange(len(dist)), dist)).rvs(size=sample_num)

def inverseTransformSample(distribution, uniformSamples):
    ''' This function lets you convert from a standard uniform sample [0,1] to
    a sample from an arbitrary distribution. This is done by taking the cdf [0,1] of 
    the arbitrary distribution, and calculating its inverse to picking the sample."
    '''
    assert (isinstance(distribution, rv_continuous) or isinstance(distribution, rv_discrete) or isinstance(distribution,rv_frozen))
    newSamples = distribution.ppf(uniformSamples)
    return newSamples

def binary():
    return rv_discrete(name="binary", values=[(0, 1), (0.5, 0.5)])

def integers(min, max):
    min = int(min)
    max = int(max)
    count = max - min + 1
    nums  = range(min, max + 1)
    probs = [1.0/count] * count
    return rv_discrete(name='integers', values=[nums, probs])

def __init__(self, param_name, X):
        self.param_name = param_name
        self.X_prime = collections.OrderedDict()
        for value in X:
            if value not in self.X_prime:
                self.X_prime[value] = 0
            self.X_prime[value] += (1.0 / len(X))
        self.distrib = rv_discrete(values=(list(range(len(self.X_prime))), list(self.X_prime.values())))

def inverseTransformSample(distribution, uniformSamples):
    ''' This function lets you convert from a standard uniform sample [0,1] to
    a sample from an arbitrary distribution. This is done by taking the cdf [0,1] of 
    the arbitrary distribution, and calculating its inverse to picking the sample."
    '''
    assert (isinstance(distribution, rv_continuous) or isinstance(distribution, rv_discrete) or isinstance(distribution,rv_frozen))
    newSamples = distribution.ppf(uniformSamples)
    return newSamples

def __call__(self, *args, **kwargs):
    frozen = self.underlying(*args, **kwargs)

    return stats.rv_discrete(
      a = 0, b = self.max_value,
      values = (np.arange(self.max_value + 1), truncated_pmf(frozen, self.max_value + 1))
    )

def __call__(self, *args, **kwargs):
    underlying_frozen = self.underlying(*args, **kwargs)
    pmf = binarize_pdf(underlying_frozen, limit=self.max_value + 1)

    return stats.rv_discrete(
      a=0,
      b=self.max_value,
      values=(np.arange(self.max_value + 1), pmf)
    )

def __init__(self, category_priors, compounds, image_shape, n_frames, max_value=1024, dtype='uint16'):
    self.category_priors = np.array(category_priors) / float(np.sum(category_priors))
    self.compounds = compounds
    self.image_shape = image_shape
    self.n_pixels = np.prod(image_shape, dtype='uint64')
    self.n_frames = n_frames
    self.max_value = max_value
    self.dtype = dtype

    self.category_distribution = stats.rv_discrete(
        name="pixel category distribution",
        values=(np.arange(self.category_priors.shape[0]), self.category_priors)
      )

def get_area_distribution(tracks, fit=False):
  area = np.sum(tracks > 0, axis=(1, 2))

  if not fit:
    count = np.bincount(area)
    probability = count / float(np.sum(count))
    return stats.rv_discrete(
      a=0,
      b=np.max(probability.shape[0]),
      name='signal distribution',
      values=(np.arange(count.shape[0]), probability)
    )
  else:
    exp_params = stats.expon.fit(area)
    return stats.expon(*exp_params)

def get_signal_distribution(tracks):
  singal = tracks[tracks > 0].ravel()
  count = np.bincount(singal)
  probability = count / float(np.sum(count))

  return stats.rv_discrete(
    a = 0,
    b = np.max(probability.shape[0]),
    name='signal distribution',
    values = (np.arange(count.shape[0]), probability)
  )

def iteration(self, node_status=True):
        """
        Execute a single model iteration

        :return: Iteration_id, Incremental node status (dictionary node->status)
        """
        self.clean_initial_status(self.available_statuses.values())
        actual_status = {node: nstatus for node, nstatus in future.utils.iteritems(self.status)}

        if self.actual_iteration == 0:
            self.actual_iteration += 1
            delta, node_count, status_delta = self.status_delta(actual_status)
            if node_status:
                return {"iteration": 0, "status": actual_status.copy(),
                        "node_count": node_count.copy(), "status_delta": status_delta.copy()}
            else:
                return {"iteration": 0, "status": {},
                        "node_count": node_count.copy(), "status_delta": status_delta.copy()}

        for u in self.graph.nodes():
            if actual_status[u] != 0:
                continue

            if self.params['model']['adopter_rate'] > 0:
                xk = (0, 1)
                pk = (1 - self.params['model']['adopter_rate'], self.params['model']['adopter_rate'])
                probability = stats.rv_discrete(name='probability', values=(xk, pk))
                number_probability = probability.rvs()

                if number_probability == 1:
                    actual_status[u] = 1
                    continue

            neighbors = self.graph.neighbors(u)
            if isinstance(self.graph, nx.DiGraph):
                neighbors = self.graph.predecessors(u)

            infected = 0
            for v in neighbors:
                infected += self.status[v]

            if infected > 0 and actual_status[u] == 0:
                eventp = np.random.random_sample()
                if eventp >= self.params['nodes']['profile'][u]:
                    actual_status[u] = 1
                else:
                    if self.params['model']['blocked'] != 0:
                        blip = np.random.random_sample()
                        if blip > self.params['model']['blocked']:
                            actual_status[u] = -1

        delta, node_count, status_delta = self.status_delta(actual_status)
        self.status = actual_status
        self.actual_iteration += 1

        if node_status:
            return {"iteration": self.actual_iteration - 1, "status": delta.copy(),
                    "node_count": node_count.copy(), "status_delta": status_delta.copy()}
        else:
            return {"iteration": self.actual_iteration - 1, "status": {},
                    "node_count": node_count.copy(), "status_delta": status_delta.copy()}