Fixed bug in T^2 likelihood; updated figures

CorrMatrix_ABC/covest.py

@@ -188,7 +188,9 @@ def par_fish_SH(n, n_D, par):
     """Parameter RMS from Fisher matrix esimation of SH likelihood.
     """
 
-    return [np.sqrt(alpha_new(n_S, n_D) * 2.0 * n_S / (n_S - 1.0)) * par for n_S in n]
+    #return [np.sqrt(alpha_new(n_S, n_D) * 2.0 * n_S / (n_S - 1.0)) * par for n_S in n]
+    return [np.sqrt(alpha_new(n_S, n_D) * (n_S - 1.0) / n_S) * par for n_S in n]
+    #return [np.sqrt((n_S - 1.0) / n_S) * par for n_S in n]
 
 
 def coeff_TJ14(n_S, n_D, n_P):
@@ -697,7 +699,7 @@ def __init__(self, par_name, n_n_S, n_R, file_base='mean_std', yscale='linear',
             self.std[p] = np.zeros(shape = (n_n_S, n_R))
 
         self.F  = np.zeros(shape = (n_n_S, n_R, 2, 2))
-        self.fs = 12
+        self.fs = 14
 
         self.n_D = n_D
         self.n_S_arr = n_S_arr
@@ -1025,7 +1027,7 @@ def plot_mean_std(self, par=None, boxwidth=None, xlog=False, ylim=None, model='a
         linestyle = ['-', '--']
 
         plot_sth = False
-        fs = 13
+        fs = 14
         plot_init(self.n_D, n_R, fs=fs)
 
         box_width = set_box_width(boxwidth, xlog, n)
@@ -1140,10 +1142,10 @@ def plot_mean_std(self, par=None, boxwidth=None, xlog=False, ylim=None, model='a
                 ylabel = stat_notation(which)
                 plt.ylabel(ylabel)
                 ax.set_yscale(self.yscale[j])
-                leg = ax.legend(loc=f'{loc[which]} right', frameon=False, handlelength=1.3)
+                leg = ax.legend(loc=f'{loc[which]} right', frameon=False, handlelength=1.3, fontsize=12)
                 plt.gca().add_artist(leg)
                 if which in leg2:
-                    ax.legend(leg2[which], lab2[which], loc=f'{loc[which]} left', frameon=False, handlelength=0.9)
+                    ax.legend(leg2[which], lab2[which], loc=f'{loc[which]} left', frameon=False, handlelength=0.9, fontsize=12)
 
                 plt.xlim(xmin, xmax)
 
@@ -1276,7 +1278,7 @@ def plot_std_var(self, par=None, sig_var_noise=None, xlog=False, title=False, mo
                         yy = self.fct['std_var_TJ14'](n_fine, self.n_D, par[i])
                         p_sym = par_symbol(p, eq=False)
                         plot_add_legend(True, n_fine, yy, linestyle=linestyle[i],
-                                        color=color[i], label=fr'$\mathbf{{\hat F}}({p_sym})$', linewidth=2)
+                                        color=color[i], label=fr'$\mathbf{{\hat F}}({p_sym})$', linewidth=12)
                         cols.append(self.fct['std_var_TJ14'](n, self.n_D, par[i]))
                         names.append('TJ14({})'.format(p))
 
@@ -1306,7 +1308,7 @@ def plot_std_var(self, par=None, sig_var_noise=None, xlog=False, title=False, mo
         plt.xlabel('$n_{\\rm s}$')
         plt.ylabel(stat_notation('std_var'))
         ax.set_yscale('log')
-        ax.legend(loc='best', numpoints=1, frameon=False, handlelength=1.3)
+        ax.legend(loc='best', numpoints=1, frameon=False, handlelength=1.3, fontsize=12)
         #ax.set_aspect(aspect=1)
 
 	    # x-ticks

results/Example_2/Euclid_dist_acf2lin/std_2mean_std_ABC.txt

@@ -1,6 +1,6 @@
-# n_S sigma(sigma^2_tilt) sigma(sigma^2_ampl) TJ14(tilt) TJ14(ampl)
-2 2.859957682e-06 1.847028056e-05 nan nan
-5 2.467277298e-06 1.186425654e-05 nan nan
-10 1.876307812e-06 8.776322349e-06 3.937050307e-07 1.100917147e-05
-20 1.509240552e-06 1.437147627e-05 3.264429663e-07 9.128322758e-06
-40 1.762607642e-06 8.17724968e-06 1.565754885e-07 4.378319471e-06
+# n_S sigma(sigma^2_tilt) sigma(sigma^2_ampl)
+2 2.859957682e-06 1.847028056e-05
+5 2.467277298e-06 1.186425654e-05
+10 1.876307812e-06 8.776322349e-06
+20 1.509240552e-06 1.437147627e-05
+40 1.762607642e-06 8.17724968e-06

results/Example_3/Euclid_dist_acf2lin/std_2mean_std_ABC.txt

@@ -1,6 +1,6 @@
-# n_S sigma(sigma^2_Omegam) sigma(sigma^2_sigma8) TJ14(Omegam) TJ14(sigma8)
-2 1.171630964e-05 2.55991623e-05 nan nan
-5 7.495316861e-06 1.374675137e-05 nan nan
-10 1.442707666e-05 3.860273721e-05 7.495544903e-06 2.503703555e-05
-20 1.173516019e-05 1.668958924e-05 6.214977511e-06 2.07596132e-05
-40 1.665798398e-05 3.413116483e-05 2.980959128e-06 9.957165308e-06
+# n_S sigma(sigma^2_Omegam) sigma(sigma^2_sigma8)
+2 1.171630964e-05 2.55991623e-05
+5 7.495316861e-06 1.374675137e-05
+10 1.442707666e-05 3.860273721e-05
+20 1.173516019e-05 1.668958924e-05
+40 1.665798398e-05 3.413116483e-05

scripts/ABC/job_ABC.py

@@ -649,8 +649,9 @@ def main(argv=None):
         yscale=['linear', 'log'],
         n_D=param.n_D,
         n_S_arr=n_S_arr,
-        fct= {'std': par_fish_SH, 'std_var_TJ14' : std_fish_deb_TJ14}
+        fct={},
     )
+       # fct={'std': par_fish_SH, 'std_var_TJ14' : std_fish_deb_TJ14}
 
 
     # Create simulations/write observation

scripts/cov_matrix_definition/cm_likelihood.py

@@ -482,7 +482,7 @@ def fit_corr_inv_true(x1, cov_true, sig2, n_jobs=3):
     toy_data['nobs'] = len(x1)     # sample size = n_D
     toy_data['x'] = x1             # explanatory variable
 
-    # cov = covariance of the data!
+    # cov = covariance of the data
     y = multivariate_normal.rvs(mean=x1, cov=cov_true, size=1)
     toy_data['y'] = y              # response variable, here one realisation
 
@@ -491,7 +491,6 @@ def fit_corr_inv_true(x1, cov_true, sig2, n_jobs=3):
     toy_data['cov_inv'] = cov_inv
 
     # STAN code
-    # the fitting code does not believe that observations are correlated!
     stan_code = """
     data {
         int<lower=0> nobs;
@@ -516,18 +515,11 @@ def fit_corr_inv_true(x1, cov_true, sig2, n_jobs=3):
     """
 
     import pystan
-    start = time.time()
     fit = pystan.stan(model_code=stan_code, data=toy_data, iter=2000, chains=n_jobs, verbose=False, n_jobs=n_jobs)
 
-    end = time.time()
-
-    #elapsed = end - start
-    #print 'elapsed time = ' + str(elapsed)
-
     return fit
 
 
-
 def fit_corr(x1, cov_true, cov_est, n_jobs=3):
     """
     Generates one draw from a multivariate normal distribution
@@ -577,21 +569,12 @@ def fit_corr(x1, cov_true, cov_est, n_jobs=3):
     """
 
     import pystan
-    start = time.time()
     fit = pystan.stan(model_code=stan_code, data=toy_data, iter=2000, chains=n_jobs, verbose=False, n_jobs=n_jobs)
 
-    # Testing: fast call to pystan
-    #fit = pystan.stan(model_code=stan_code, data=toy_data, iter=1, chains=1, verbose=False, n_jobs=n_jobs)
-
-    end = time.time()
-
-    #elapsed = end - start
-    #print 'elapsed time = ' + str(elapsed)
-
     return fit
 
 
-def fit_corr_SH(x1, cov_true, cov_est_inv, n_jobs=3):
+def fit_corr_SH(x1, cov_true, cov_est_inv, n_S, n_jobs=3):
     """
     Generates one draw from a multivariate student-t distribution
     (see Sellentin & Heavens 2015)
@@ -601,6 +584,7 @@ def fit_corr_SH(x1, cov_true, cov_est_inv, n_jobs=3):
     input:  x1, mean of multivariate normal distribution - vector of floats
             cov_true, square covariance matrix for the simulation
             cov_est_inv, inverse square covariance matrix for the fitting
+            n_S, number of simulations used for covariance
             n_jobs, number of parallel jobs
 
     output: fit, Stan fitting object
@@ -610,19 +594,20 @@ def fit_corr_SH(x1, cov_true, cov_est_inv, n_jobs=3):
     toy_data = {}                  # build data dictionary
     toy_data['nobs'] = len(x1)     # sample size = n_D
     toy_data['x'] = x1             # explanatory variable
+    toy_data['ns'] = n_S           # number of simulations used for covariance
 
-    # cov = covariance of the data!
+    # cov = covariance of the data
     y = multivariate_normal.rvs(mean=x1, cov=cov_true, size=1)
     toy_data['y'] = y              # response variable, here one realisation
 
     # set estimated inverse covariance matrix for fitting
     toy_data['cov_est_inv'] = cov_est_inv
 
     # STAN code
-    # the fitting code does not believe that observations are correlated!
     stan_code = """
     data {
-        int<lower=0> nobs;                                 
+        int<lower=0> nobs;
+        int<lower=0> ns;
         vector[nobs] x;                       
         vector[nobs] y;   
         matrix[nobs, nobs] cov_est_inv;                    
@@ -637,30 +622,24 @@ def fit_corr_SH(x1, cov_true, cov_est_inv, n_jobs=3):
 
         mu = b + a * x;
 
-        # Sellentin & Heavens debiasing scheme:
-        # Replace normal likelihood with t-distribution
-        # y ~ (1 + (log_normal(mu, cov_est)) / (1 + n_S))^(-n_S/2)  
+        // Sellentin & Heavens debiasing scheme:
+        // Replace normal likelihood with t-distribution
+        // y ~ (1 + (log_normal(mu, cov_est)) / (n_S - 1))^(-n_S/2)  
 
         chi2 = (y - mu)' * cov_est_inv * (y - mu);
 
-        # target += log-likelihood
-        target += log(1.0 + chi2/(1.0 + nobs)) * -nobs/2.0;
+        // target += log-likelihood
+        target += log(1.0 + chi2/(ns - 1.0)) * -ns/2.0;
     }
     """
 
 
     import pystan
-    start = time.time()
     fit = pystan.stan(model_code=stan_code, data=toy_data, iter=2000, chains=n_jobs, verbose=False, n_jobs=n_jobs)
-    end = time.time()
-
-    #elapsed = end - start
-    #print 'elapsed time = ' + str(elapsed)
 
     return fit
 
 
-
 def set_fit_MCMC(fit, res, i, run):
     """Set mean and std of fit according to stan output res.
 
@@ -774,7 +753,7 @@ def simulate(x1, yreal, n_S_arr, sigma_ML, sigma_m1_ML, sigma_m1_ML_deb, fish_an
                 if re.search('SH', options.likelihood) is not None:
                     if options.verbose == True:
                         print('Running MCMC with Sellentin&Heavens (SH) likelihood')
-                    res = fit_corr_SH(x1, cov, cov_est_inv, n_jobs=options.n_jobs)
+                    res = fit_corr_SH(x1, cov, cov_est_inv, n_S, n_jobs=options.n_jobs)
                     set_fit_MCMC(fit_SH, res, i, run)
 
 
@@ -940,22 +919,22 @@ def main(argv=None):
 
 
     # Initialisation of results
-    sigma_ML = Results(tr_name, n_n_S, options.n_R, file_base='sigma_ML')
-    sigma_m1_ML = Results(tr_name, n_n_S, options.n_R, file_base='sigma_m1_ML', yscale='log', fct={'mean': tr_N_m1_ML})
-    sigma_m1_ML_deb = Results(tr_name, n_n_S, options.n_R, file_base='sigma_m1_ML_deb', fct={'mean': no_bias})
+    sigma_ML = Results(tr_name, n_n_S, options.n_R, file_base='sigma_ML', n_S_arr=n_S_arr)
+    sigma_m1_ML = Results(tr_name, n_n_S, options.n_R, file_base='sigma_m1_ML', yscale='log', fct={'mean': tr_N_m1_ML}, n_S_arr=n_S_arr)
+    sigma_m1_ML_deb = Results(tr_name, n_n_S, options.n_R, file_base='sigma_m1_ML_deb', fct={'mean': no_bias}, n_S_arr=n_S_arr)
 
-    fish_ana = Results(par_name, n_n_S, options.n_R, file_base='std_Fisher_ana', yscale='log', fct={'std': par_fish})
+    fish_ana = Results(par_name, n_n_S, options.n_R, file_base='std_Fisher_ana', yscale='log', fct={'std': par_fish}, n_S_arr=n_S_arr)
 
     fish_num = Results(par_name, n_n_S, options.n_R, file_base='std_Fisher_num', yscale='log', \
-                       fct={'std': par_fish, 'std_var_TJK13': std_fish_biased_TJK13, 'std_var_TJ14': std_fish_biased_TJ14})
+                       fct={'std': par_fish, 'std_var_TJK13': std_fish_biased_TJK13, 'std_var_TJ14': std_fish_biased_TJ14}, n_S_arr=n_S_arr)
     fish_deb = Results(par_name, n_n_S, options.n_R, file_base='std_Fisher_deb', yscale='log', \
-                       fct={'std': no_bias, 'std_var_TJK13': std_fish_deb, 'std_var_TJ14': std_fish_deb_TJ14})
+                       fct={'std': no_bias, 'std_var_TJK13': std_fish_deb, 'std_var_TJ14': std_fish_deb_TJ14}, n_S_arr=n_S_arr)
     fit_norm_num = Results(par_name, n_n_S, options.n_R, file_base='mean_std_fit_norm', yscale=['linear', 'log'],
-                       fct={'std': par_fish, 'std_var_TJK13': std_fish_biased_TJK13, 'std_var_TJ14': std_fish_biased_TJ14})
+                       fct={'std': par_fish, 'std_var_TJK13': std_fish_biased_TJK13, 'std_var_TJ14': std_fish_biased_TJ14}, n_S_arr=n_S_arr)
     fit_norm_deb = Results(par_name, n_n_S, options.n_R, file_base='mean_std_fit_norm_deb', yscale=['linear', 'log'],
-                       fct={'std': no_bias, 'std_var_TJ14': std_fish_deb_TJ14})
+                       fct={'std': no_bias, 'std_var_TJ14': std_fish_deb_TJ14}, n_S_arr=n_S_arr)
     fit_SH = Results(par_name, n_n_S, options.n_R, file_base='mean_std_fit_SH', yscale=['linear', 'log'],
-                       fct={'std': par_fish_SH})
+                       fct={'std': par_fish_SH}, n_S_arr=n_S_arr)
 
     # Data
     x1 = uniform.rvs(loc=-delta/2, scale=delta, size=options.n_D)        # exploratory variable
@@ -994,17 +973,17 @@ def main(argv=None):
         print('Creating plots')
 
     if options.do_fish_ana == True:
-        fish_ana.plot_mean_std(n_S_arr, options.n_D, par={'std': dpar_exact})
-    fish_num.plot_mean_std(n_S_arr, options.n_D, par={'std': dpar_exact})
-    fish_deb.plot_mean_std(n_S_arr, options.n_D, par={'std': dpar_exact})
+        fish_ana.plot_mean_std(par={'std': dpar_exact})
+    fish_num.plot_mean_std(par={'std': dpar_exact})
+    fish_deb.plot_mean_std(par={'std': dpar_exact})
 
     if options.do_fit_stan == True:
         if re.search('norm_biased', options.likelihood) is not None:
-            fit_norm_num.plot_mean_std(n_S_arr, options.n_D, par={'mean': options.par, 'std': dpar_exact}, boxwidth=param.boxwidth)
+            fit_norm_num.plot_mean_std(par={'mean': options.par, 'std': dpar_exact}, boxwidth=param.boxwidth)
         if re.search('norm_deb', options.likelihood) is not None:
-            fit_norm_deb.plot_mean_std(n_S_arr, options.n_D, par={'mean': options.par, 'std': dpar_exact}, boxwidth=param.boxwidth)
+            fit_norm_deb.plot_mean_std(par={'mean': options.par, 'std': dpar_exact}, boxwidth=param.boxwidth)
         if re.search('SH', options.likelihood) is not None:
-            fit_SH.plot_mean_std(n_S_arr, options.n_D, par={'mean': options.par, 'std': dpar_exact}, boxwidth=param.boxwidth)
+            fit_SH.plot_mean_std(par={'mean': options.par, 'std': dpar_exact}, boxwidth=param.boxwidth)
 
     dpar2 = dpar_exact**2
 
@@ -1013,24 +992,24 @@ def main(argv=None):
     else:
         title = False
 
-    fish_num.plot_std_var(n_S_arr, options.n_D, par=dpar2, title=title)
+    fish_num.plot_std_var(par=dpar2, title=title)
 
-    fish_deb.plot_std_var(n_S_arr, options.n_D, par=dpar2, title=title)
+    fish_deb.plot_std_var(par=dpar2, title=title)
     if options.do_fit_stan == True:
         if re.search('norm_biased', options.likelihood) is not None:
-            fit_norm_num.plot_std_var(n_S_arr, options.n_D, par=dpar2, sig_var_noise=param.sig_var_noise, title=title)
+            fit_norm_num.plot_std_var(par=dpar2, sig_var_noise=param.sig_var_noise, title=title)
         if re.search('norm_deb', options.likelihood) is not None:
-            fit_norm_deb.plot_std_var(n_S_arr, options.n_D, par=dpar2, sig_var_noise=param.sig_var_noise, title=title)
+            fit_norm_deb.plot_std_var(par=dpar2, sig_var_noise=param.sig_var_noise, title=title)
         if re.search('SH', options.likelihood) is not None:
-            fit_SH.plot_std_var(n_S_arr, options.n_D, par=dpar2, sig_var_noise=param.sig_var_noise, title=title)
+            fit_SH.plot_std_var(par=dpar2, sig_var_noise=param.sig_var_noise, title=title)
 
-    sigma_ML.plot_mean_std(n_S_arr, options.n_D, par={'mean': [options.sig2]}, boxwidth=param.boxwidth, ylim=[4.9, 5.1])
+    sigma_ML.plot_mean_std(par={'mean': [options.sig2]}, boxwidth=param.boxwidth, ylim=[4.9, 5.1])
 
     # Precision matrix trace
     f_mean = 1/options.sig2
 
-    sigma_m1_ML.plot_mean_std(n_S_arr, options.n_D, par={'mean': [f_mean]}, boxwidth=param.boxwidth)
-    sigma_m1_ML_deb.plot_mean_std(n_S_arr, options.n_D, par={'mean': [f_mean]}, boxwidth=param.boxwidth)
+    sigma_m1_ML.plot_mean_std(par={'mean': [f_mean]}, boxwidth=param.boxwidth)
+    sigma_m1_ML_deb.plot_mean_std(par={'mean': [f_mean]}, boxwidth=param.boxwidth)
 
     ### End main program
 

scripts/plot_many.py

@@ -317,9 +317,10 @@ def plot_box(fits, n_D, par, dy, which, boxwidth=None, xlog=False, ylog=False,
 
             # xticks
             for n_S in n:
-                x_loc.append(n_S)
-                lab = '{}'.format(n_S)
-                x_lab.append(lab)
+                if n_S not in (974, 2094, 4502):
+                    x_loc.append(n_S)
+                    lab = '{}'.format(n_S)
+                    x_lab.append(lab)
 
             # Curves for predicted values
             if fit.fct is not None and which in fit.fct:
@@ -398,6 +399,8 @@ def plot_box(fits, n_D, par, dy, which, boxwidth=None, xlog=False, ylog=False,
             frac = float(n_D) / float(n_S)
             if frac > 100:
                 lab = '{:.0f}'.format(frac)
+            elif frac < 1:
+                lab = '{:.1g}'.format(frac)
             else:
                 lab = '{:.2g}'.format(frac)
             x2_lab.append(lab)
@@ -429,15 +432,13 @@ def main(argv=None):
     # Read sampling results
     fit_ABC = read_fit(param.par_name, param.ABC, 'ABC',
                        fct={'mean': no_bias}, flab={'mean': 'true value'}, verbose=param.verbose)
-                             # fct={'std': no_bias, 'std_var': std_fish_deb_TJ14},
     fit_MCMC_norm = read_fit(param.par_name, param.MCMC_norm, 'MCMC $N$',
                              fct={'std': no_bias, 'std_var': std_fish_Gupta}, #std_fish_deb},
                              flab={'std': r'$\mathbf{\hat F}$', 'std_var': r'$\mathbf{\hat F}$'},
                              verbose=param.verbose)
-                             #flab={'std': r'$\hat{\bm{\mathrm{F}}}$', 'std_var': r'$\hat{\bm{\mathrm{F}}}$'},
     fit_MCMC_T2 = read_fit(param.par_name, param.MCMC_T2, 'MCMC $T^2$',
-                           fct={'std': par_fish_SH}, flab={'std': r'$\mathbf{\hat F}_{T^2}$'}, verbose=param.verbose)
-                           #fct={'std': par_fish_SH}, flab={'std': r'$\hat{\bm{\mathrm{F}}}_{T^2}$'}, verbose=param.verbose)
+                           fct={}, flab={}, verbose=param.verbose)
+                           #fct={'std': par_fish_SH}, flab={'std': r'$\mathbf{\hat F}_{T^2}$'}, verbose=param.verbose)
 
     mode = -1
     delta = 200