[docs]defset_initial(self,distribution=None):ifdistributionisNone:# assume standard normal by defaultifself.domainisnotNone:# assume uniform if domain specifiedmn=np.min(self.domain,axis=1)mx=np.max(self.domain,axis=1)distribution=dist.uniform(loc=mn,scale=mx-mn)else:distribution=dist.norm()initial_dist=distributionself._in=initial_dist.pdf(self.X).prod(axis=1)self._up=Noneself._pr=None
[docs]defset_predicted(self,distribution=None,**kwargs):ifdistributionisNone:# Reweight kde of predicted by weights from previous iteration if presentdistribution=gkde(self.y.T,**kwargs,weights=self._weights)pred_pdf=distribution.pdf(self.y.T).Telse:pred_pdf=distribution.pdf(self.y,**kwargs)self._pr=pred_pdfself._up=None
[docs]deffit(self,**kwargs):ifself._inisNone:self.set_initial()self._pr=Noneifself._prisNone:self.set_predicted(**kwargs)ifself._obisNone:self.set_observed()# Store ratio of observed/predicted# e.g. to comptue E(r) and to pass on to future iterationsself._r=np.divide(self._ob,self._pr)# If present, multiply by weights (e.g. from a previous iteration)ifself._weightsisnotNone:self._r*=self._weights# Multiply by initial to get updated pdfup_pdf=np.multiply(self._in,self._r)self._up=up_pdf
[docs]classBayesProblem(object):""" Sets up Bayesian Inverse Problem for parameter identification Example Usage ------------- >>> from mud.base import BayesProblem >>> import numpy as np >>> from scipy.stats import distributions as ds >>> X = np.random.rand(100,1) >>> num_obs = 50 >>> Y = np.repeat(X, num_obs, 1) >>> y = np.ones(num_obs)*0.5 + np.random.randn(num_obs)*0.05 >>> B = BayesProblem(X, Y, np.array([[0,1], [0,1]])) >>> B.set_likelihood(ds.norm(loc=y, scale=0.05)) >>> np.round(B.map_point()[0],1) 0.5 """def__init__(self,X,y,domain=None):self.X=Xify.ndim==1:y=y.reshape(-1,1)self.y=yself.domain=domainself._ps=Noneself._pr=Noneself._ll=None
[docs]defset_likelihood(self,distribution,log=False):iflog:self._log=Trueself._ll=distribution.logpdf(self.y).sum(axis=1)# below is an equivalent evaluation (demonstrating the expected symmetry)# std, mean = distribution.std(), distribution.mean()# self._ll = dist.norm(self.y, std).logpdf(mean).sum(axis=1)else:self._log=Falseself._ll=distribution.pdf(self.y).prod(axis=1)# equivalent# self._ll = dist.norm(self.y).pdf(distribution.mean())/distribution.std()# self._ll = self._ll.prod(axis=1)self._ps=None
[docs]defset_prior(self,distribution=None):ifdistributionisNone:# assume standard normal by defaultifself.domainisnotNone:# assume uniform if domain specifiedmn=np.min(self.domain,axis=1)mx=np.max(self.domain,axis=1)distribution=dist.uniform(loc=mn,scale=mx-mn)else:distribution=dist.norm()prior_dist=distributionself._pr=prior_dist.pdf(self.X).prod(axis=1)self._ps=None