1- """An example implementation of a Gaussian noise model ."""
1+ """Example implementations of Gaussian noise models ."""
22
33from functools import partial
44
55import numpy as np
66import scipy .stats as ss
77
88import elfi
9+ from elfi .examples .gnk import euclidean_multidim
910
1011
11- def Gauss (mu , sigma , n_obs = 50 , batch_size = 1 , random_state = None ):
12- """Sample the Gaussian distribution.
12+ def gauss (mu , sigma , n_obs = 50 , batch_size = 1 , random_state = None ):
13+ """Sample the 1-D Gaussian distribution.
1314
1415 Parameters
1516 ----------
1617 mu : float, array_like
1718 sigma : float, array_like
1819 n_obs : int, optional
1920 batch_size : int, optional
20- random_state : RandomState, optional
21+ random_state : np.random.RandomState, optional
22+
23+ Returns
24+ -------
25+ y_obs : array_like
26+ 1-D observation.
27+
28+ """
29+ # Handling batching.
30+ batches_mu = np .asanyarray (mu ).reshape ((- 1 , 1 ))
31+ batches_sigma = np .asanyarray (sigma ).reshape ((- 1 , 1 ))
32+
33+ # Sampling observations.
34+ y_obs = ss .norm .rvs (loc = batches_mu , scale = batches_sigma ,
35+ size = (batch_size , n_obs ), random_state = random_state )
36+ return y_obs
37+
38+
39+ def gauss_nd_mean (* mu , cov_matrix , n_obs = 15 , batch_size = 1 , random_state = None ):
40+ """Sample an n-D Gaussian distribution.
41+
42+ Parameters
43+ ----------
44+ *mu : array_like
45+ Mean parameters.
46+ cov_matrix : array_like
47+ Covariance matrix.
48+ n_obs : int, optional
49+ batch_size : int, optional
50+ random_state : np.random.RandomState, optional
51+
52+ Returns
53+ -------
54+ y_obs : array_like
55+ n-D observation.
2156
2257 """
23- # Standardising the parameter's format.
24- mu = np .asanyarray (mu ).reshape ((- 1 , 1 ))
25- sigma = np .asanyarray (sigma ).reshape ((- 1 , 1 ))
26- y = ss .norm .rvs (loc = mu , scale = sigma , size = (batch_size , n_obs ), random_state = random_state )
27- return y
58+ n_dim = len (mu )
59+
60+ # Handling batching.
61+ batches_mu = np .zeros (shape = (batch_size , n_dim ))
62+ for idx_dim , param_mu in enumerate (mu ):
63+ batches_mu [:, idx_dim ] = param_mu
64+
65+ # Sampling the observations.
66+ y_obs = np .zeros (shape = (batch_size , n_obs , n_dim ))
67+ for idx_batch in range (batch_size ):
68+ y_batch = ss .multivariate_normal .rvs (mean = batches_mu [idx_batch ],
69+ cov = cov_matrix ,
70+ size = n_obs ,
71+ random_state = random_state )
72+ if n_dim == 1 :
73+ y_batch = y_batch [:, np .newaxis ]
74+ y_obs [idx_batch , :, :] = y_batch
75+ return y_obs
2876
2977
3078def ss_mean (x ):
@@ -39,36 +87,71 @@ def ss_var(x):
3987 return ss
4088
4189
42- def get_model (n_obs = 50 , true_params = None , seed_obs = None ):
43- """Return a complete Gaussian noise model.
90+ def get_model (n_obs = 50 , true_params = None , seed_obs = None , nd_mean = False , cov_matrix = None ):
91+ """Return a Gaussian noise model.
4492
4593 Parameters
4694 ----------
4795 n_obs : int, optional
48- the number of observations
4996 true_params : list, optional
50- true_params[0] corresponds to the mean,
51- true_params[1] corresponds to the standard deviation
97+ Default parameter settings.
5298 seed_obs : int, optional
53- seed for the observed data generation
99+ Seed for the observed data generation.
100+ nd_mean : bool, optional
101+ Option to use an n-D mean Gaussian noise model.
102+ cov_matrix : None, optional
103+ Covariance matrix, a requirement for the nd_mean model.
54104
55105 Returns
56106 -------
57107 m : elfi.ElfiModel
58108
59109 """
110+ # Defining the default settings.
60111 if true_params is None :
61- true_params = [10 , 2 ]
112+ if nd_mean :
113+ true_params = [4 , 4 ] # 2-D mean.
114+ else :
115+ true_params = [4 , .4 ] # mean and standard deviation.
62116
63- y_obs = Gauss (* true_params , n_obs = n_obs , random_state = np .random .RandomState (seed_obs ))
64- sim_fn = partial (Gauss , n_obs = n_obs )
117+ # Choosing the simulator for both observations and simulations.
118+ if nd_mean :
119+ sim_fn = partial (gauss_nd_mean , cov_matrix = cov_matrix , n_obs = n_obs )
120+ else :
121+ sim_fn = partial (gauss , n_obs = n_obs )
122+
123+ # Obtaining the observations.
124+ y_obs = sim_fn (* true_params , n_obs = n_obs , random_state = np .random .RandomState (seed_obs ))
65125
66126 m = elfi .ElfiModel ()
67- elfi .Prior ('uniform' , - 10 , 50 , model = m , name = 'mu' )
68- elfi .Prior ('truncnorm' , 0.01 , 5 , model = m , name = 'sigma' )
69- elfi .Simulator (sim_fn , m ['mu' ], m ['sigma' ], observed = y_obs , name = 'Gauss' )
70- elfi .Summary (ss_mean , m ['Gauss' ], name = 'S1' )
71- elfi .Summary (ss_var , m ['Gauss' ], name = 'S2' )
72- elfi .Distance ('euclidean' , m ['S1' ], m ['S2' ], name = 'd' )
127+
128+ # Initialising the priors.
129+ eps_prior = 5 # The longest distance from the median of an initialised prior's distribution.
130+ priors = []
131+ if nd_mean :
132+ n_dim = len (true_params )
133+ for i in range (n_dim ):
134+ name_prior = 'mu_{}' .format (i )
135+ prior_mu = elfi .Prior ('uniform' , true_params [i ] - eps_prior ,
136+ 2 * eps_prior , model = m , name = name_prior )
137+ priors .append (prior_mu )
138+ else :
139+ priors .append (elfi .Prior ('uniform' , true_params [0 ] - eps_prior ,
140+ 2 * eps_prior , model = m , name = 'mu' ))
141+ priors .append (elfi .Prior ('truncnorm' ,
142+ np .amax ([.01 , true_params [1 ] - eps_prior ]),
143+ 2 * eps_prior , model = m , name = 'sigma' ))
144+ elfi .Simulator (sim_fn , * priors , observed = y_obs , name = 'gauss' )
145+
146+ # Initialising the summary statistics.
147+ sumstats = []
148+ sumstats .append (elfi .Summary (ss_mean , m ['gauss' ], name = 'ss_mean' ))
149+ sumstats .append (elfi .Summary (ss_var , m ['gauss' ], name = 'ss_var' ))
150+
151+ # Choosing the discrepancy metric.
152+ if nd_mean :
153+ elfi .Discrepancy (euclidean_multidim , * sumstats , name = 'd' )
154+ else :
155+ elfi .Distance ('euclidean' , * sumstats , name = 'd' )
73156
74157 return m
0 commit comments