Skip to content

Hessian PDF Covariance Matrix for theory predictions #1743

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
comane merged 6 commits into from
Jun 9, 2023
Merged

Hessian PDF Covariance Matrix for theory predictions #1743

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
e31db1f
added hessian covariance matrix for pdf errors
comane May 25, 2023
ebc95f8
run black on covmats.py+ add test
comane May 30, 2023
7751162
added new check for either replicas or symmhessian type
comane May 30, 2023
13896bf
added check to covmats
comane May 30, 2023
7efa7fd
added regression table for hessian test
comane May 31, 2023
88df97f
use for checks
comane Jun 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions validphys2/src/validphys/checks.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,16 @@ def check_pdf_is_montecarlo(ns, **kwargs):
raise CheckError(f"Error type of PDF {pdf} must be 'replicas' and not {etype}")


@make_argcheck
def check_pdf_is_montecarlo_or_symmhessian(ns, **kwargs):
pdf = ns['pdf']
etype = pdf.error_type
check(
etype in {'replicas', 'symhessian'},
f"Error type of PDF {pdf} must be either 'replicas' or 'symmhessian' and not {etype}",
)


@make_check
def check_know_errors(ns, **kwargs):
pdf = ns['pdf']
Expand Down
29 changes: 24 additions & 5 deletions validphys2/src/validphys/covmats.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
check_data_cuts_match_theorycovmat,
check_dataset_cuts_match_theorycovmat,
check_norm_threshold,
check_pdf_is_montecarlo,
check_pdf_is_montecarlo_or_symmhessian,
check_speclabels_different,
)
from validphys.commondata import loaded_commondata_with_cuts
Expand Down Expand Up @@ -677,11 +677,15 @@ def groups_corrmat(groups_covmat):
return mat


@check_pdf_is_montecarlo
Copy link
Copy Markdown
Contributor

@Zaharid Zaharid May 26, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should still be some check, for the error member types we don't support.

@check_pdf_is_montecarlo_or_symmhessian
def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
"""For a given `dataset`, returns the sum of the covariance matrix given by
`covmat_t0_considered` and the PDF error: a covariance matrix estimated from the
replica theory predictions from a given monte carlo `pdf`
`covmat_t0_considered` and the PDF error:
- If the PDF error_type is 'replicas', a covariance matrix is estimated from
the replica theory predictions
- If the PDF error_type is 'symmhessian', a covariance matrix is estimated using
formulas from (mc2hessian) https://arxiv.org/pdf/1505.06736.pdf


Parameters
----------
Expand Down Expand Up @@ -716,7 +720,22 @@ def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
True
"""
th = ThPredictionsResult.from_convolution(pdf, dataset)
pdf_cov = np.cov(th.error_members, rowvar=True)

if pdf.error_type == 'replicas':
pdf_cov = np.cov(th.error_members, rowvar=True)

elif pdf.error_type == 'symmhessian':
rescale_fac = pdf._rescale_factor()
hessian_eigenvectors = th.error_members
central_predictions = th.central_value

# need to subtract the central set which is not the same as the average of the
# Hessian eigenvectors.
X = hessian_eigenvectors - central_predictions.reshape((central_predictions.shape[0], 1))
# need to rescale the Hessian eigenvectors in case the eigenvector confidence interval is not 68%
X = X / rescale_fac
pdf_cov = X @ X.T

return pdf_cov + covmat_t0_considered


Expand Down
Loading