Add to MCMC

NSSEA/__constraints.py

@@ -94,6 +94,8 @@
 import scipy.stats as sc
 import pandas as pd
 import xarray as xr
+from cmdstanpy import CmdStanModel
+
 
 from.__multi_model import MultiModel
 from .__tools import matrix_squareroot
@@ -328,6 +330,46 @@ def _constrain_law_keep( climIn , Yo , keep , n_mcmc_drawn_min , n_mcmc_drawn_ma
 	pass
 ##}}}
 
+def _constrain_law_ess( climIn , Yo , n_mcmc_drawn_min , n_mcmc_drawn_max , verbose , **kwargs ):##{{{
+	#Warning : Sample are saved given another coordinate : "sample_MCMC" with X sample attached
+	clim = climIn.copy()
+
+	pb = ProgressBar( clim.n_sample + 1 , "constrain_law" , verbose )
+
+	n_ess = kwargs.get("n_ess")
+	if n_ess is None:
+		n_ess = 10	
+	min_rate_accept = 0
+	ns_params_names = clim.ns_law.get_params_names()
+	n_features = clim.ns_law.n_ns_params
+	## Define prior
+	prior_mean   = clim.data["mm_mean"][-clim.n_coef:].values
+	prior_cov    = clim.data["mm_cov"][-clim.n_coef:,-clim.n_coef:].values
+	prior_law    = sc.multivariate_normal( mean = prior_mean , cov = prior_cov , allow_singular = True )
+
+	#OUtput
+	results=np.array([])
+	sample_names =[s+"_"+str(i) for i in range(n_ess) for s in clim.sample[1:]]+["BE"]
+
+	law_coef_bay   = xr.DataArray( np.zeros( (n_features,(clim.n_sample)*n_ess + 1,1) ) , coords = [ ns_params_names , sample_names , ["Multi_Synthesis"] ] , dims = ["coef","sample_MCMC","model"] )
+	## And now MCMC loop
+	for s in clim.sample[1:]:
+		pb.print()
+		X   = clim.X.loc[Yo.index,s,"F","Multi_Synthesis"].values.squeeze()
+		n_mcmc_drawn = np.random.randint( n_mcmc_drawn_min , n_mcmc_drawn_max )
+		draw = clim.ns_law.drawn_bayesian( Yo.values.squeeze() , X , n_mcmc_drawn , prior_law , min_rate_accept , **kwargs )
+		n_tirage=(len(draw)//n_ess)
+		select=draw[0::n_tirage][:n_ess]
+		law_coef_bay.loc[:,[s+"_"+str(i) for i in range(n_ess)],"Multi_Synthesis"]=select.T
+
+	clim.law_coef=law_coef_bay
+	clim.law_coef.loc[:,"BE",:] = clim.law_coef[:,1:,:].median( dim = "sample_MCMC" )
+	clim.BE_is_median = True
+
+	pb.end()
+
+	return clim
+
 def constrain_law( climIn , Yo , keep = "all" , n_mcmc_drawn_min = 5000 , n_mcmc_drawn_max = 10000 , verbose = False , **kwargs ):##{{{
 	"""
 	NSSEA.constrain_law
@@ -338,7 +380,7 @@ def constrain_law( climIn , Yo , keep = "all" , n_mcmc_drawn_min = 5000 , n_mcmc
 	---------
 	climIn : [NSSEA.Climatology] clim variable
 	Yo       : [pandas.DataFrame] Observations of ns_law
-	keep     : [ "all" or a float between 0 and 1] If keep < 1, only a ratio of 
+	keep     : [ "all", "ess", or a float between 0 and 1] If keep < 1, only a ratio of 
 	          keep covariates is used, and many coefficients are drawn for the
 	          same covariate. Faster, but can reduce confidence interval
 	          uncertainty.
@@ -350,8 +392,10 @@ def constrain_law( climIn , Yo , keep = "all" , n_mcmc_drawn_min = 5000 , n_mcmc
 	------
 	clim : [NSSEA.Climatology] A copy is returned
 	"""
-
-	if keep == "all" or not keep < 1:
+	if keep == "ess":
+		return _constrain_law_ess( climIn , Yo , n_mcmc_drawn_min , n_mcmc_drawn_max , verbose , **kwargs )	
+	elif keep == "all" or not keep < 1:
+		#Maybe avoid a possible comparison between str and int ? put percentage as a kwarg ?
 		return _constrain_law_all( climIn , Yo , n_mcmc_drawn_min , n_mcmc_drawn_max , verbose , **kwargs )
 	else:
 		return _constrain_law_keep( climIn , Yo , keep , n_mcmc_drawn_min , n_mcmc_drawn_max , verbose , **kwargs )
@@ -625,5 +669,79 @@ def constraint_C0( climIn , Yo , verbose = False ): ##{{{
 	return climIn.copy()
 ##}}}
 
+# Version originale Contrainte STAN
+#Only for GEV original (covariable for mu and sigma, exponential sigma)
+def stan_constrain(climIn,Yo,stan_file, **kwargs):
+	"""
+	NSSEA.constrain_law
+	===================
+	Constrain the law_coef of the clim with a MCMC approach.
+
+	Arguments
+	---------
+	climIn : [NSSEA.Climatology] clim variable
+	Yo       : [pandas.DataFrame] Observations of ns_law
+	keep     : [ "all", "ess", or a float between 0 and 1] If keep < 1, only a ratio of 
+	          keep covariates is used, and many coefficients are drawn for the
+	          same covariate. Faster, but can reduce confidence interval
+	          uncertainty.
+	n_mcmc_draw_min: [integer] Minimum number of coef to draw for each covariate
+	n_mcmc_draw_max: [integer] Maximum number of coef to draw for each covariate
+	verbose  : [bool] Print (or not) state of execution
+
+	Return
+	------
+	clim : [NSSEA.Climatology] A copy is returned
+	"""
+	clim      = climIn.copy()
+
+	#data
+	X=clim.X.loc[Yo.index,'BE',"F","Multi_Synthesis"].values.squeeze()
+	N_X=len(X)
+	Y=Yo.values.squeeze()
+	N=len(Y)
+	p_m   = clim.data["mm_mean"][-clim.n_coef:].values
+	p_cov    = clim.data["mm_cov"][-clim.n_coef:,-clim.n_coef:].values
+
+	u,s,v = np.linalg.svd(p_cov)
+	p_std = u @ np.sqrt(np.diag(s)) @ v.T
+
+	#clim para
+	n_features = clim.ns_law.n_ns_params
+	ns_params_names = clim.ns_law.get_params_names()
+	n_ess = kwargs.get("n_ess")
+	if n_ess is None:
+        	n_ess = 10
+
+	sample_names =[s+"_"+str(i) for i in range(n_ess) for s in clim.sample[1:]]+["BE"]
+	law_coef_bay   = xr.DataArray( np.zeros( (n_features,(clim.n_sample)*n_ess + 1,1) ) ,
+		coords = [ ns_params_names , sample_names , ["Multi_Synthesis"] ] ,
+		dims = ["coef","sample_MCMC","model"] )
+
+
+
+	newDF = pd.DataFrame() #creates a new dataframe that's empty
+	samples = clim.X.loc[Yo.index,:,"F","Multi_Synthesis"].sample.values.squeeze()[1::]
 
+	#compile stan file
+	model = CmdStanModel(stan_file=stan_file)
+	#Run constraint with varying covariable
+	for s in samples:
+		print(s)
+		X=clim.X.loc[Yo.index,s,"F","Multi_Synthesis"].values.squeeze()
+		data = {"N":N,
+		"Y": Y,
+		"N_X": N_X,
+		"X": X,
+		"p_m":p_m,
+		"p_std":p_std
+		}
+		fit = model.sample(data=data,iter_sampling=int(n_ess/4),show_progress=False)
+		df = fit.draws_pd()
+		sub_def=df[['para[1]','para[2]','para[3]','para[4]','para[5]']]
+		law_coef_bay.loc[:,[s+"_"+str(i) for i in range(n_ess)], "Multi_Synthesis"] = sub_def.T
+	clim.law_coef=law_coef_bay
+	clim.law_coef.loc[:,"BE",:] = clim.law_coef[:,1:,:].median( dim = "sample_MCMC" )
+	clim.BE_is_median = True
+	return clim
 

NSSEA/__covariates.py

@@ -261,7 +261,7 @@ def error_distribution(self):
 
 ##}}}
 
-def covariates_FC_GAM( clim , lX , XN , dof = 7 , method = "pygam" , verbose = False ):##{{{
+def covariates_FC_GAM( clim , lX , XN , dof = 7 , method = "pygam" , verbose = False, light=False ):##{{{
 	"""
 	NSSEA.covariates_FC_GAM
 	=======================
@@ -288,15 +288,16 @@ def covariates_FC_GAM( clim , lX , XN , dof = 7 , method = "pygam" , verbose = F
 	method : [str] "statsmodels" or "pygam", select the package used to solve
 	         GAM model
 	verbose : [bool] If we print the progress of the fit or not.
+	light: [bool] Only Best estimate or several samples. Only Best Estimate only possible for pygam
 
 	"""
 	if method == "pygam":
-		return _covariates_FC_GAM_pygam( clim , lX , XN , dof , verbose )
+		return _covariates_FC_GAM_pygam( clim , lX , XN , dof , verbose, light )
 	else:
 		return _covariates_FC_GAM_statsmodels( clim , lX , XN , dof , verbose )
 ##}}}
 
-def _covariates_FC_GAM_pygam( clim , lX , XN , dof = 7 , verbose = False ):##{{{
+def _covariates_FC_GAM_pygam( clim , lX , XN , dof = 7 , verbose = False, light=False ):##{{{
 	"""
 	NSSEA.covariates_FC_GAM_pygam
 	=============================
@@ -340,20 +341,20 @@ def _covariates_FC_GAM_pygam( clim , lX , XN , dof = 7 , verbose = False ):##{{{
 		## Distribution of GAM coefficients
 		gam_law = gam_model.error_distribution()
 		coefs_  = gam_law.rvs(n_sample)
-
-		for i,s in enumerate(samples[1:]):
-			pb.print()
+		if not light:
+			for i,s in enumerate(samples[1:]):
+				pb.print()
 
-			xn = XN.values[:,i+1]
-			XF = np.stack( (time  ,xn) , -1 )
-			XC = np.stack( (time_C,xn) , -1 )
+				xn = XN.values[:,i+1]
+				XF = np.stack( (time  ,xn) , -1 )
+				XC = np.stack( (time_C,xn) , -1 )
 
-			## Perturbation
-			gam_model.coef_ = coefs_[i,:]
+				## Perturbation
+				gam_model.coef_ = coefs_[i,:]
 
-			## Final decomposition
-			dX.loc[:,s,"F",model] = gam_model.predict( XF )
-			dX.loc[:,s,"C",model] = gam_model.predict( XC )
+				## Final decomposition
+				dX.loc[:,s,"F",model] = gam_model.predict( XF )
+				dX.loc[:,s,"C",model] = gam_model.predict( XC )
 
 	clim.X = dX
 	pb.end()

NSSEA/__init__.py

@@ -116,6 +116,7 @@
 from .__nsstats     import statistics_fixed_IF
 
 from .__nsstats     import build_params_along_time
+from .__nsstats     import build_params_along_time_ess
 from .__nsstats     import add_return_time
 from .__nsstats     import add_FAR
 from .__nsstats     import add_bias
@@ -124,7 +125,7 @@
 from .__constraints import constrain_covariate
 from .__constraints import constrain_law
 from .__constraints import constraint_C0
-
+from .__constraints import stan_constrain
 
 
 

NSSEA/__multi_model.py

@@ -129,7 +129,6 @@ def _fit( self , mm_matrix ):##{{{
 		cov_S = np.zeros( (n_params,n_params) )
 		for i in range(n_models):
 			cov_S += np.cov( mm_matrix[:,1:,i] )
-
 		SSM     = np.cov( mm_matrix[:,0,:] ) * ( n_models - 1 )
 		cov_CMU = matrix_positive_part( SSM / ( n_models - 1 ) - cov_S / n_models )
 		self.cov  = ( n_models + 1 ) / n_models * cov_CMU + cov_S / n_models**2
@@ -152,13 +151,20 @@ def fit( self , mm_matrix ):##{{{
 		self._fit(mm_matrix)
 	##}}}
 
+	def rvs_old(self):##{{{
+		"""
+		Return a random sample from multi model
+		"""
+		### Depreciated for scipy 1.11.1
+		return self.mean + self.std @ np.random.normal(size = self.mean.size) 
+	##}}}
 	def rvs(self):##{{{
 		"""
 		Return a random sample from multi model
 		"""
-		return self.mean + self.std @ np.random.normal(size = self.mean.size)
+		### Added for scipy 1.11.1
+		return np.random.default_rng().multivariate_normal( mean=self.mean, cov = self.cov)
 	##}}}
-
 	## Properties {{{
 
 	@property

NSSEA/__nsfit.py

@@ -94,11 +94,12 @@
 from .__tools import ProgressBar
 
 
+
 ###############
 ## Functions ##
 ###############
 
-def nslaw_fit( lY , clim , verbose = False ):
+def nslaw_fit( lY , clim , verbose = False, light=False ):
 	"""
 	NSSEA.nslaw_fit
 	===============
@@ -139,21 +140,22 @@ def nslaw_fit( lY , clim , verbose = False ):
 		law = clim.ns_law
 		law.fit(Y.values,X.values)
 		law_coef.loc[:,"BE",model] = law.get_params()
-
-		for s in sample:
-			pb.print()
+		coef_be = law_coef.loc[:,"BE",model].values.ravel()
+		if not light:
+			for s in sample:
+				pb.print()
 
-			fit_is_valid = False
-			while not fit_is_valid:
+				fit_is_valid = False
+				while not fit_is_valid:
 
-				idx = np.random.choice( tY.size , tY.size , replace = True )
+					idx = np.random.choice( tY.size , tY.size , replace = True )
 
-				tYs = tY.values[idx]
-				Ys = Y.iloc[idx].values
-				Xs = clim.X.loc[tYs,s,"F",model].values
-				law.fit(Ys,Xs)
-				fit_is_valid = law.check( Y.values.squeeze() , X.values.squeeze() , np.arange( 0 , tY.size , 1 ) )
-			law_coef.loc[:,s,model] = law.get_params()
+					tYs = tY.values[idx]
+					Ys = Y.iloc[idx].values
+					Xs = clim.X.loc[tYs,s,"F",model].values
+					law.fit(Ys,Xs, init = coef_be)
+					fit_is_valid = law.check( Y.values.squeeze() , X.values.squeeze() , np.arange( 0 , tY.size , 1 ) )
+				law_coef.loc[:,s,model] = law.get_params()
 
 	clim.law_coef = law_coef
 	pb.end()

NSSEA/__nsstats.py

@@ -555,7 +555,7 @@ def build_params_along_time( clim , verbose = False ):##{{{
 	"""
 	NSSEA.extremes_stats
 	====================
-	Build trajectories of params alon time
+	Build trajectories of params along time
 
 	Arguments
 	---------
@@ -596,4 +596,49 @@ def build_params_along_time( clim , verbose = False ):##{{{
 ##}}}
 
 
+def build_params_along_time_ess( clim , verbose = False ):##{{{
+	"""
+	NSSEA.extremes_stats
+	====================
+	Build trajectories of params along time
+	Adaptated for MCMC sample changes with ess
+	Beware: Only function I adapted
+
+	Arguments
+	---------
+	clim : NSSEA.Climatology
+		A clim variable
+	verbose: bool
+		Print state of execution or not
+
+	Return
+	------
+	params : xr.DataArray
+		An array containing params along time
+
+	"""
+	ns_law = clim.ns_law
+
+	l_params = [k for k in clim.ns_law.lparams]
+	xrdims   = ["time","sample","forcing","param","model"]
+	xrcoords = [clim.time,clim.law_coef.sample_MCMC,["F","C"],l_params,clim.model]
+	s_params = xr.DataArray( np.zeros( (clim.n_time,len(clim.law_coef.sample_MCMC),2,len(l_params),clim.n_model) ) , dims = xrdims , coords = xrcoords )
+
+
+	pb = ProgressBar(  clim.n_model * (len(clim.law_coef.sample_MCMC)) , "build_params_along_time" , verbose = verbose )
+	for m in clim.model:
+		for s in s_params.sample:
+			pb.print()
+			s_X=str(s.values).split("_")[0]
+			clim.ns_law.set_params( clim.law_coef.loc[:,s,m].values )
+			for f in s_params.forcing:
+				clim.ns_law.set_covariable( clim.X.loc[clim.time,s_X,f,m].values , clim.time )
+				for p in l_params:
+					s_params.loc[:,s,f,p,m] = clim.ns_law.lparams[p](clim.time)
+
+
+	if verbose: pb.end()
+
+	return s_params
+