# Licensed under a 3-clause BSD style license - see LICENSE.rstimporthtmlimportnumpyasnpfrom.utilsimportget_random_state__all__=["InverseCDFSampler"]
[docs]classInverseCDFSampler:"""Inverse CDF sampler. It determines a set of random numbers and calculate the cumulative distribution function. Parameters ---------- pdf : `~gammapy.maps.Map` Map of the predicted source counts. axis : int Axis along which sampling the indexes. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. """def__init__(self,pdf,axis=None,random_state=0):self.random_state=get_random_state(random_state)self.axis=axisifaxisisnotNone:self.cdf=np.cumsum(pdf,axis=self.axis)self.cdf/=self.cdf[:,[-1]]else:self.pdf_shape=pdf.shapepdf=pdf.ravel()/pdf.sum()self.sortindex=np.argsort(pdf,axis=None)self.pdf=pdf[self.sortindex]self.cdf=np.cumsum(self.pdf)def_repr_html_(self):try:returnself.to_html()exceptAttributeError:returnf"<pre>{html.escape(str(self))}</pre>"
[docs]defsample_axis(self):"""Sample along a given axis. Returns ------- index : tuple of `~numpy.ndarray` Coordinates of the drawn sample. """choices=self.random_state.uniform(high=1,size=len(self.cdf))shape_cdf=self.cdf.shapecdf_all=np.insert(self.cdf,0,0,axis=1)edges=np.arange(shape_cdf[1]+1)-0.5pix_coords=[]forcdf,choiceinzip(cdf_all,choices):pix=np.interp(choice,cdf,edges)pix_coords.append(pix)returnnp.array(pix_coords)
[docs]defsample(self,size):"""Draw sample from the given PDF. Parameters ---------- size : int Number of samples to draw. Returns ------- index : tuple of `~numpy.ndarray` Coordinates of the drawn sample. """# pick numbers which are uniformly random over the cumulative distribution functionchoice=self.random_state.uniform(high=1,size=size)# find the indices corresponding to this point on the CDFindex=np.searchsorted(self.cdf,choice)index=self.sortindex[index]# map back to multi-dimensional indexingindex=np.unravel_index(index,self.pdf_shape)index=np.vstack(index)index=index+self.random_state.uniform(low=-0.5,high=0.5,size=index.shape)returnindex