<figure><img src="image_1.png" /></figure>

HAL
open science

# Bootstrap Method with Fisher Formalism in Cosmology

Fabien Dournac

► To cite this version:

Fabien Dournac. Bootstrap Method with Fisher Formalism in Cosmology. 2024. hal-04723110

HAL Id: hal-04723110

https://hal.science/hal-04723110v1

Preprint submitted on 27 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Bootstrap Method with Fisher Formalism in Cosmology

Fabien Dournac

October 2024

## 1 Introduction

This report details the use of the Bootstrap method in a cosmology context ($w_0 w_a$ CDM model investigated for example with EUCLID survey where we have a spectroscopic probe and a photometric one) to combine two Fisher matrices (F1.txt and F2.txt) in order to analyze the associated uncertainties. Each Fisher matrix contains information of 7 cosmological parameters ($\Omega_m, \Omega_b, w_0, w_a, h, n_s, \sigma_8$). Bootstrap samples are generated using two different methods: the first is based on the combinations of the Fisher matrices' eigenvectors, and the second uses linear combinations of eigenvectors with Dirichlet weights.

## 2 Definition of Figure-of-Merit (FoM)

In cosmology, the Figure of Merit (FoM) is a metric used to quantify the constraining power of a survey or experiment, particularly when analyzing dark energy parameters such as $w_0$ and $w_a$ in the $w_0 w_a$ CDM model. Here, $w_0$ represents the present-day value of the dark energy equation of state, and $w_a$ describes its time evolution.

The FoM is defined as the inverse of the square root of the determinant of the covariance matrix of $w_0$ and $w_a$, and is given by the following equation:

$$ \text{FoM} = \frac{1}{\sqrt{\det(\text{Cov}(w_0, w_a))}} $$

Where: - Cov($w_0, w_a$) is the covariance matrix of the dark energy parameters $w_0$ and $w_a$. - The determinant of the covariance matrix, det(Cov($w_0, w_a$)), corresponds to the area of the error ellipse in the $w_0-w_a$ plane.

A smaller determinant corresponds to a smaller error ellipse, indicating tighter constraints on the parameters. Hence, a higher FoM value signifies stronger constraints on the dark energy parameters, offering a more precise understanding of the role of dark energy in cosmic expansion.

The FoM is widely used to compare the effectiveness of different cosmological surveys, helping to assess and optimize the precision of parameter estimation.

## 3 Description of the Bootstrap Algorithm for Fisher Matrices

### 3.1 Introduction to Fisher's Formalism

The Fisher Information Matrix (FIM) is a positive-definite symmetric matrix that quantifies the amount of information a set of observations contains about the parameters of a statistical model. It is denoted $I(\theta)$, where $\theta$ represents the parameter vector.

Mathematically, $I_{ij} = -E\left[\frac{\partial^2 \ln f_{\theta_0}(X)}{\partial \theta_i \partial \theta_j}\right]$, where $f$ is the likelihood function. The $ij^{th}$ entry of the FIM is:

$$ I_{ij} = \text{Cov} \left[ \frac{\partial \ln f_{\theta_0}(X)}{\partial \theta_i}, \frac{\partial \ln f_{\theta_0}(X)}{\partial \theta_j} \right] $$

Consequently, the FIM is the covariance matrix between score functions (derivatives of the log-likelihood function with respect to parameters $\theta_i$ and $\theta_j$). As a covariance matrix, the FIM's eigenvalues and eigenvectors indicate the principal components and main directions of variance. The primary directions provide the $\theta^*$ values, representing high-value data points where the model should focus. Conversely, the inverse of the Fisher matrix, denoted $I^{-1}(\theta) = \text{Cov}[\hat{\theta}_i, \hat{\theta}_j]$, provides estimates of parameter variances and covariances.

In applied statistics, the Bootstrap is a tool that enables inference by evaluating the variability of the statistics of interest directly from the data without relying on asymptotic theory. In scenarios where asymptotic theoretical derivations are complex or only feasible under restrictive assumptions, the Bootstrap is particularly attractive. In the context of the Bootstrap for Fisher matrices, it is crucial to use techniques that preserve their inherent characteristics,
such as positive definiteness and symmetry, during resampling. For the specific scenario where we have two Fisher matrices, A and B, corresponding to identical parameters, the Bootstrap can be achieved by manipulating the eigenvalues and eigenvectors. These manipulations involve using the original matrices, A and B, as sources for modifications.

## 3.2 Algorithm Description

### 1. Eigenvalue Decomposition:

Perform eigenvalue decomposition on both matrices A and B. This yields: For A: Eigenvalues $\Lambda_A = (\lambda_A^1, \lambda_A^2, ..., \lambda_A^n)$ and corresponding eigenvectors $V_A = (v_A^1, v_A^2, ..., v_A^n)$. For B: Eigenvalues $\Lambda_B = (\lambda_B^1, \lambda_B^2, ..., \lambda_B^n)$ and corresponding eigenvectors $V_B = (v_B^1, v_B^2, ..., v_B^n)$.

### 2. Generation of Permutation Masks:

To define the order in which the eigenvalues will be swapped, we create a vector of masks m of size n (where n is the number of eigenvalues). When n is small, these mask vectors can be efficiently generated by converting a sequence of integers from 0 to $2^n - 1$ into binary masks.

### 3. Swapping Eigenvalues and Eigenvectors:

Loop through the mask vector m:

For each index i in m: Select the eigenvalue: $\lambda_C^i = \begin{cases} \lambda_A^i & m(i) = 0 \\ \lambda_B^i & m(i) = 1 \end{cases}$ (where m(i) is the value at index i in the permutation vector). Select the corresponding eigenvectors: $v_C^i = \begin{cases} v_A^i & m(i) = 0 \\ v_B^i & m(i) = 1 \end{cases}$.

### 4. Reconstruction:

After manipulating the eigenvalues and eigenvectors, the modified matrix (denoted as C) can be reconstructed as:

$$C = V_C \Lambda_C V_C^T$$

where $V_C$ represents the combined eigenvectors and $\Lambda_C$ is the diagonal matrix containing the altered eigenvalues.

### 5. Aggregation:

The Bootstrap replicates generated, denoted as $C_0, C_1, ..., C_{2^{n-1}}$, serve as the basis for aggregation. A simple approach is to compute the average of these replicates:

$$\hat{C} = \frac{1}{2^n} \sum_{i=0}^{2^n-1} C_i.$$

# 4 Methodology

## 4.1 Loading the Matrices

The Fisher matrices are loaded from text files:

```python
def load_matrix(path):
    return np.loadtxt(path, delimiter="-")
```

## 4.2 QR Decomposition for Orthonormalization

A QR decomposition is performed to orthonormalize the matrices:

```python
def gram_schmidt(A):
    Q, _ = np.linalg.qr(A)
    return Q
```
## 4.3 Checking for Positive Definiteness

Cholesky decomposition is used to verify if a matrix is positive definite:

def is_posixitive_definite (matrix):

try:

- = np. linalg. cholesky (matrix)

return True

except np. linalg. LinAlgError:

return False

## 4.4 Generating Bootstrap Samples Using Eigenvector Combinations

Bootstrap samples are generated by combining the eigenvectors of the Fisher matrices. This method uses all possible combinations of the eigenvectors of the two matrices:

def build_Bootstrap_sample (qd_decompositions, orthonormal_basis = True ):

n = qd_decompositions [0][0]. size

n Combinations = 2n

samples = []

for k in range (n Combinations ):

binary_str = bin (k)[2:] .zfill(n)

mask = np. array([int(bit) for bit in binary_str])

eigenvecs = np. column_stack(

[qd_decompositions [mask[i] ][1][ : , i]

for i in range (len (mask)) ])

)

eigenvals = [qd_decompositions [mask[i] ][0][i]

for i in range (len (mask)) ])

if orthonormal_basis :

eigenvecs = gram_schmidt(eigenvecs)

sample = np. dot(np. dot(eigenvecs, np. diag(eigenvals)),

np. transpose(eigenvecs))

if is_posixitive_definite (sample):

samples. append(sample)

return samples

There are 2⁷ = 128 samples generated using this initial method, which is based on all possible combinations of the eigenvectors of the two Fisher matrices.

## 4.5 Generating Additional Bootstrap Samples with Dirichlet Weights

Additional samples are generated by combining the initial samples with Dirichlet weights. The Dirichlet weights are
used because they preserve the norm at 1, allowing for convex combinations of the initial samples:

def generate_additionalSamples(initialSamples, n_additionalSamples):

additionalSamples = []

n_initial_sample = len (initial_sample)

for _ in range (n_additional_sample):

weights = np. random. dirichlet(

np. ones(n_initial_sample), size =1)[0]

additional_sample = np. average(

initial_sample, axis =0, weights =weights)

if is_posixitive_definite (additional_sample):

additional_sample. append (additional_sample)

return additional_sample

In this analysis, 10,000 additional samples are generated, bringing the total number of samples to 128 + 10,000 =
10,128.

5 Calculating the Mean of Covariance Matrices in the Bootstrap

In the context of estimating constraints on the parameters, we used the Bootstrap method to generate multiple
resampled samples. For each sample, we calculated the associated covariance matrix.
The mean of the covariance matrices is calculated as follows:

$$ \overline{\text{Cov}} = \frac{1}{n} \sum_{i=1}^{n} \text{Cov}_i $$

where n represents the number of Bootstrap resamples.

As we seek to obtain an arithmetic mean of the covariances, we divide by n. This approach is consistent with standard statistical methods for the mean of estimates.

# 6 Results

## 6.1 Histogram of the Figure of Merit (FoM)

The histogram below shows the distribution of the FoM values after removing outliers.

<figure><img src="image_2.png" /><figcaption>Figure 1: Histogram of the FoM after removing outliers - FoM = 3463 ± 388</figcaption></figure>

## 6.2 Sample Statistics

The Bootstrap sample statistics are calculated, including covariance matrices and the FoM. Outliers are removed to obtain a cleaner distribution.

average_matrix, covariance, FoM_stat = collect_statistics(allSamples)
FoM_stat_cleaned = remove_outliers(np.array(FoM_stat))

## 6.3 Variance of Variances

The variance of variances for each parameter is calculated from a large number of Bootstrap samples.

mean_variances, std_variances, error_on_error = calculate_variances_of_variances(covariance, n_Bootstrap=1000)
<table><caption>Table 1: Estimated parameters and errors (1 sigma) obtained through Bootstrap - FoM = 3463 ± 388</caption><thead><tr><th>Parameter</th><th>Estimated Variance ± Error (1 sigma)</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.0017607253886199 ± 5.55984110e-08</td></tr><tr><td>Ω<sub>b</sub></td><td>0.0008017141842887 ± 6.49595749e-08</td></tr><tr><td>w<sub>0</sub></td><td>0.0124673873467692 ± 1.97120496e-07</td></tr><tr><td>w<sub>a</sub></td><td>0.0265941394218139 ± 2.66356392e-07</td></tr><tr><td>h</td><td>0.0011305980414546 ± 1.89673603e-07</td></tr><tr><td>n<sub>s</sub></td><td>0.0023873741304972 ± 1.20929837e-07</td></tr><tr><td>σ<sub>8</sub></td><td>0.0015490070479626 ± 1.03860327e-07</td></tr></tbody></table>

# 7 Table of Constraints and Variances on Variances

# 8 Conclusion

The use of Dirichlet weights to generate Bootstrap samples helps preserve the norm at 1, ensuring an appropriate linear combination of the eigenvectors. The results show that removing outliers leads to a more centered and more Gaussian-like distribution of the FoM. The initial analysis based on the $2^7 = 128$ samples provides a solid basis, enriched by the generation of 10,000 additional samples.

# 9 Appendices

## 9.1 Complete Code

import os, sys

import numpy as np

import matplotlib.pylab as plt

import subprocess

def load_matrix(path):

"""Load a matrix from a text file."""

return np.loadtxt(path, delimiter="-")

def gram_schmidt(A):

"""Perform QR decomposition for orthonormalization of matrix A."""

Q, _ = np Gram .qr(A)

return Q

def is_posi tive_definite (matrix):

"""Check if a matrix is positive definite using Cholesky decomposition."""

try:

_ = np linalg cholesky (matrix)

return True

except np linalg LinAlgError:

return False

def build_Bootstrap_sample(qd_decompositions, orthonormal_basis = True):

"""Generate Bootstrap samples from given eigen decompositions."""

n = qd_decompositions[0][0].size

ncombinations = 2n

samples = []

for k in range(ncombinations):

binary_str = bin(k)[2:].zfill(n)

mask = np .array([int(bit) for bit in binary_str])

eigenvecs = np .column_stack(

[ qd_decompositions [mask [i] ] [1] [: , i]

for i in range(len (mask)) ])

)

eigenvals = [ qd_decompositions [mask [i] ] [0] [i]
for i in range(len mask)]

if orthonormal:

eigenvecs = gram_schmidt(eigenvecs)

sample = np.dot(np.dot(eigenvecs, np Eag( eigenvals )),
np transpose eigenvects)

if is-positive-definite (sample):

samples.append(sample)

return samples

def generate_additionalSamples(initialSamples, nAdditionalSamples):

"""Generate additional Bootstrap samples by combining initial samples."""

additionalSamples = []

nInitialSamples = len(initialSamples)

for in range(nAdditionalSamples):

weights = np.random. dirichlet(
np. ones(nInitialSamples), size =1)[0]

additional_sample = np. average(
initial Samples, axis =0, weights =weights)

if is-positive-definite (additional_sample):

additionalSamples.append(additional_sample)

return additionalSamples

def collect_statistics (samples):

"""Calculate statistics on samples, including covariance matrices."""

covariance matrices = [np. lINALg. inv(sample) for sample in samples]

FoM_stat = [
np. sqrt(1.0 / (np. lINALg. det(cov [2:4, 2:4]))))
for cov in covariances

}

average_matrix = np. average samples, axis =0)

np. savetxt("Fisher_average_FoM.txt", average_matrix)

return average_matrix, covariances, FoM_stat

def calculate_variances_of_variances (covariance佔, n_Bootstrap =1000):

"""Calculate the variance of variances for each parameter

with a large number of Bootstrap samples."""

variances = [np. diag(cov) for cov in covariance佔s]

variances_array = np. array(variances)

Bootstrap Means = []

Bootstrap_vars = []

for in range(n_Bootstrap):

sampled indices = np.random.choice(
len (variances_array), len (variances_array), replace =True)

sampled_variances = variances_array[sampled_indices]

Bootstrap Means. append(np. mean(sampled_variances, axis =0))

Bootstrap_vars. append(np. var(sampled_variances, axis =0))

Bootstrap Means = np. array(Bootstrap Means)

Bootstrap_vars = np. array(Bootstrap_vars)

mean_variances = np. mean(Bootstrap Means, axis =0)

std_variances = np. std(Bootstrap Means, axis =0)

error_on_error = std_variances / np. sqrt(n_Bootstrap)

return mean_variances, std_variances, error_on_error

def plot_FoM_histogram (FoM_stat):

"""Plot a histogram of the Figure of Merit."""

plt. hist(FoM_stat, bins =100)

plt.xlabel("Value")
```python
plt.ylabel("Frequency")
plt.title("FoM-Histogram")
plt.show()
```

**def save_results(F_final):**
"""Save the final Fisher matrix and its covariance."""
F_finalInverse = np凝.lndalg.凝inv(F_final)
np savetxt("F_final.txt", F_final)
np savetxt("F_finalInverse.txt", F_finalInverse)

**def remove_outliers(data, m=2):**
"""Remove outliers from a dataset using the Z-score method."""
return data [np凝.abs(data - np凝.mean(data)) < m * np凝.std(data)]

# Main processing sequence

matrix_paths = ["./F1.txt", "/F2.txt"]
matrices = [load_matrix(path) for path in matrix_paths]
qd_decompositions = [np凝.lndalg.凝eig(matrix) for matrix in matrices]
initialSamples = build_Bootstrap_sample(qd_decompositions)
additionalSamples = generate_additionalSamples(
    initialSamples, n_additionalSamples=10000)
allSamples = initialSamples + additionalSamples
average_matrix, covariance佔mrices, FoM_stat = collect_statistics(allSamples)

# Remove outliers from FoM_stat

FoM_stat_cleaned = remove_outliers(np凝.array(FoM_stat))

covariance_average = np凝.lndalg.凝inv(average佔matrix)
final_variances = np凝 sqrt(np凝.diag(covariance_average))
mean_variances, std_variances, error_on_error = calculate_variances_of_variances(
    covariance佔mrices, n_Bootstrap=1000)

# Display the results

parameter_names = ["wm", "wb", "w0", "wa", "h", "ns", "s8"]
Bootstrap_constraints = np凝 sqrt(np凝.diag(np凝.lndalg.凝inv(average佔matrix)))

print("Estimated variance +/-- Error (1sigma) for each parameter (Bootstrap):")
for i, name in enumerate(parameter_names):
    print(f"{name} +/-- {Bootstrap_constraints[i] : .16 f} +/-- {error_on_error[i] : .8 e}")

plot_FoM_histogram(FoM_stat_cleaned)
save_results(average佔matrix)

# Checking in case - constraints + FoM

# os.system("printFoM F_final.txt OPT GCph 0 0 F N")
FoM_final = float(subprocess.凝check_output(
    'printFoM F_final.txt OPT GCph 0 0 F N | grep FoM',
    |
    | awk '{F=' + '}' | cut -f2 | 'python -m foM'
    '))
print('\n-FoM=' , FoM_final, '+/-' , np.凝std(FoM_stat), '\n')

# 10 Why Use the Bootstrap on Fisher Rather Than Covariance Matrices?

The Bootstrap is a resampling method that allows estimating the distribution of an estimator or a statistic using samples generated from the observed data. In the context of Fisher tests, Bootstrap is often preferred over methods based on covariance matrices for several reasons:

1. Robustness to Assumptions: The Fisher test relies on strict assumptions, such as homogeneity of variance-covariances and multivariate normality of data. If these assumptions are not met, the test results may be biased. Bootstrap, however, does not require such assumptions, making it more robust.
2. Adaptation to Small Sample Sizes: When the sample size is small, estimates based on covariance matrices can be unreliable. Bootstrap, by directly resampling observations, provides more stable and reliable estimates even for small samples.

3. Better Estimation of the Statistic's Distribution: Bootstrap allows for directly estimating the empirical distribution of the test statistic (in this case, the Fisher criterion), whereas the classical approach relies on asymptotic approximations that may not be accurate under certain conditions.

These advantages make Bootstrap a particularly suitable method for evaluating Fisher's statistic, especially in contexts where classical assumptions are difficult to verify.

# 11 Dirichlet Weights: Analytical Formula and Usage

In the context of Bootstrap, Dirichlet weights play a crucial role in generating weighted samples. The Dirichlet distribution is a probability distribution over vectors of positive numbers that sum to one. It is used to model the weights assigned to each observation during resampling.

## 11.1 Analytical Formula of the Dirichlet Distribution

The Dirichlet distribution with parameter $\alpha = (\alpha_1, \alpha_2, \dots, \alpha_k)$, where $\alpha_i > 0$, is defined over a vector of probabilities $(p_1, p_2, \dots, p_k)$ such that $p_i > 0$ and $\sum_{i=1}^k p_i = 1$. Its probability density function is given by:

$$f(p_1, p_2, \dots, p_k; \alpha_1, \alpha_2, \dots, \alpha_k) = \frac{\Gamma\left(\sum_{i=1}^k \alpha_i\right)}{\prod_{i=1}^k \Gamma(\alpha_i)} \prod_{i=1}^k p_i^{\alpha_i - 1}$$

where $\Gamma(\cdot)$ is the gamma function.

## 11.2 Using Dirichlet Weights in Bootstrap

In the context of Bootstrap, Dirichlet weights are used to weight observations during resampling. Suppose we have a sample of size $n$. To generate a weighted Bootstrap sample, we draw a probability vector $(p_1, p_2, \dots, p_n)$ according to a Dirichlet distribution with parameter $\alpha = (1, 1, \dots, 1)$, i.e., with $\alpha_i = 1$ for all $i$. This choice of parameter generates weights that are uniform on average, but with variability due to the resampling process.

The new weighted observations $X_i^*$ are then obtained by assigning each initial observation $X_i$ a weight $p_i$. This process can be formally written as:

$$X_i^* = p_i \cdot X_i \quad \text{for } i = 1, 2, \dots, n$$

where the $p_i$ follow the Dirichlet distribution mentioned above. This mechanism captures the uncertainty of resampling while accounting for potential variations in the contributions of different observations.