Fisher information Matrix and Principal Component Analysis

Caroline B. Owen

This notebook gives and overview of Fisher information matrices and principal component analysis in the context of gravitational wave astronomy. 

A Fisher information matrix (FIM) is a way to approximate the shape of a probability distribution as a Gaussian distribution. We use it to make forecasts about  about what kinds of measurements we will be able to make with different models and detectors. 

See arXiv:gr-qc/0703086  or arXiv:1502.05735 for an overviews. 
See arXiv:2510.11861 for an example of using a Fisher analysis to make predictions about what we can learn from future detectors.  
(This paper makes use of the python package GWBench for computing  FIMs. See arXiv:2010.15202)
See arXiv:1808.05286 for an example of using a Fisher analysis to make predictions about the constraints we will be able to place on dark matter.

Set the working directory

Set the working directory to the git repo. This notebook will access files in the L05_bayes and L07_FIM_PCA directories

SetDirectory@"/Users/cbo/Documents/Work/Classes/CFCGWA_UNIMIB";

Multivariate Gaussian distributions

The probability distribution function (PDF) on a  k-dimentional random vector θ with  mean vector μ and covarriance matrix  Σ  is

𝒩[θ|μ,Σ] =((Exp[ -(1/2)[(((θ-μ))^T)(Σ^(-1))(θ-μ)] ])/( (Sqrt[(((2π))^k)Det[Σ]])))

the diagonal elements of Σ contain the variances ((Subscript[σ, i])^2) of the marginal distributions

(Subscript[Σ, ii]) = ((Subscript[σ, i])^2)

and the off-diagonal elements are proportional to the correlation coefficients

(Subscript[Σ, ij]) = (Subscript[ρ, ij ])(Subscript[σ, i])(Subscript[σ, j])

The second derivatives of the natural log of the PDF with respect to the components of θ give us the components of (Σ^(-1))

(∂/(∂(θ^i)))(∂/(∂(θ^j)))log 𝒩[θ|μ,Σ] = -(Subscript[((Σ^(-1))), ij])

Below, we confirm this property for a 2 D distribution

Σ = {{((σ1)^2),ρ σ1 σ2},{ρ σ1 σ2,((σ2)^2)}};
μ = {μ1, μ2};
θ = {θ1,θ2};

pdf = PDF[MultinormalDistribution[μ,Σ],θ]

(*Take the log and second derivatives*)
Table[D[Log[pdf],θ[[aa]],θ[[bb]]],{aa,1,2},{bb,1,2}]//Simplify
(* Compare to (Σ^(-1))*)
%/Inverse[Σ]//Simplify

We can use this property to estimate the variances and pairwise correlations of distributions that are approximately Gaussian.

L05 toy problem

In the toy problem in L05, our log likelihood was

log ℒ[d|θ]=(∑_(i=1)^N)log 𝒩[(Subscript[x, i])|0,σ]

with θ = {σ}

Assuming the likelihood is approximately Gaussian, we can  estimate the standard deviation on the model parameters θ

-(∂/(∂(θ^i)))(∂/(∂(θ^j)))log ℒ[d|θ]  = -(∑_(i=1)^N)(∂/(∂(θ^i)))(∂/(∂(θ^j)))log 𝒩[(Subscript[x, i])|0,σ]

We only have a single parameter, so

-((∂^2)/(∂(σ^2)))log ℒ[d|θ]  = -(∑_(i=1)^N)((∂^2)/(∂(σ^2)))log 𝒩[(Subscript[x, i])|0,σ]

The estimated standard deviation on our measurement of σ is

(Subscript[σ, σ]) = (([-(∑_(i=1)^N)((∂^2)/(∂(σ^2)))log 𝒩[(Subscript[x, i])|0,σ]])^(-1/2))

We can compute this estimate for the specific data realization we considered

data = Flatten[Import["L05_bayes/data-1.csv"]];

(* define the PDF*)
 pdf = PDF[NormalDistribution[0,σ],x]
(* take the log and simplify *)
 FullSimplify[Log[%],x∈Reals &&σ>0]
(* define the summand *)
tosum = - D[%,σ,σ]
(* take the sum and evaluate at the injection value*)
Sum[tosum,{x,data}]/.σ->3
(* compute (Subscript[σ, σ]) *)
1/Sqrt[%]

If our samplers did a good job, this should be approximately the width of the posteriors we obtained (because we  assumed uniform priors).  And, in fact, it is! Here we load posterior samples from L05. (The original notebook did not save the samples, so you can either rerun the updated version or just download the sample file from the repo.)

posteriorTrue = Import["L05_bayes/posterior_MCMC.csv"]//Flatten;
Sqrt[Variance[posteriorTrue]]

However, the expression above is not the Fisher Information Matrix (FIM). The FIM is not a statement about a specific data realization but, instead, an expectation value taken over the data generation process:

(Γ^ij) = -⟨(∂/(∂(θ^i)))(∂/(∂(θ^i)))log ℒ(d|θ)⟩

Our data is draw from a normal distribution 𝒩[(Subscript[x, i])|0,σ]. Each term in the sum must be integrated over the PDF

(Γ^(σσ))=  -(∑_(i=1)^N) ∫𝒩[(Subscript[x, i])|0,σ] ((∂^2)/(∂(σ^2)))log 𝒩[(Subscript[x, i])|0,σ]dx

The sum can be replaced by a factor of N,

(Γ^(σσ))=- N ∫𝒩[(Subscript[x, i])|0,σ] ((∂^2)/(∂(σ^2)))log 𝒩[(Subscript[x, i])|0,σ]dx

This FIM only has one element so inverting it simple,

(Subscript[σ, σ]) = ((((Γ^(σσ))))^(-1/2))

We can see that

(Subscript[σ, σ])  ~ (N^(-1/2))

This shows explicitly how increasing  the number of data points will improve the precision  σ

Length[data]Integrate[tosum pdf,{x,-∞,∞},Assumptions->σ>0]
σσ = 1/Sqrt[%] /.σ->3//N

We can approximate the posterior as a Gaussian centered at the true value (Subscript[μ, σ]) = 3 and standard deviation (Subscript[σ, σ])

posteriorFisher = PDF[NormalDistribution[3,σσ],x];

What is causing the shift between the true posterior and the Fisher approximation? How could this be improved?

Binary Parameter Estimation FIM

The binary likelihood is

log ℒ[d|θ] ~ -(1/2)(d-h[θ]|d-h[θ])

and the FIM is

(Γ^ij) = -⟨(∂/(∂(θ^i)))(∂/(∂(θ^j)))log ℒ(d|θ)⟩

		=⟨-(∂/(∂(θ^i)))(d-h[θ]| ((∂h)/(∂(θ^j)))) ⟩

		=⟨(((∂h)/(∂(θ^i)))|((∂h)/(∂(θ^j))))-(d-h[θ]|(((∂^2)h)/(∂(θ^i)∂(θ^j))))⟩

We assume that d = n + h[θ], where n is the detector noise and that n is drawn from a Gaussian distribution with zero mean. 
Therefore,

⟨d-h[θ]⟩ =⟨ n + h[θ]-h[θ]⟩ =0 

	(Γ^ij) = (((∂h)/(∂(θ^i)))|((∂h)/(∂(θ^j))))

When the following criteria are met, (Γ^(-1)) is exactly the covariance matrix of the posterior PDF
- The SNR is large 
- The waveform is a linear function of the model parameters
- The priors are wide and flat on the relevant area of parameter space

These assumptions are equivalent to the posterior PDF being approximately  Gaussian.  See arXiv:gr-qc/0703086 for proofs of these criteria. 

These criteria are never exactly met for any realistic waveform, but we can still use (Γ^(-1)) to approximate the credible intervals we would obtain from a full Bayesian analysis as well as the correlations between pairs of parameters. 


Recall that the SNR is ρ  = (Sqrt[(h|h)]).  
For a given waveform and location in parameter space, how does Γ scale with SNR? 
How will the inferred credible intervals scale with SNR?

mass and spin parameter conversions

The waveform is written in terms of various mass and spin parameters for convenience reasons, but it is useful to have rules to convert expressions into the pairs of parameters we are interested in

m1m2Rule ={m1->m1,m2->m2,m->m1+m2,η->((m1 m2)/(((m1+m2))^2)),ℳ->((((m1 m2))^(3/5))/(((m1+m2))^(1/5))),q->(m2/m1),δ->((m1-m2)/(m1+m2))};
mηRule = {m1->(1/2) m (1+(Sqrt[1-4 η])),m2->(1/2) m (1-(Sqrt[1-4 η])),m->m,η->η,ℳ->m (η^(3/5)),q->((1-(Sqrt[1-4 η]))/(1+(Sqrt[1-4 η]))),δ->(Sqrt[1-4 η])};
ℳqRule = {m1->(((((1+q)/(q^3))))^(1/5)) ℳ,m2->((((q^2) (1+q)))^(1/5)) ℳ,m->(((((1+q))^(6/5)) ℳ)/(q^(3/5))),η->(q/(((1+q))^2)),ℳ->ℳ,q->q,δ->((1-q)/(1+q))};
χ1χ2Rule = {χ1->χ1,χ2->χ2,χ2->χ2,χs->((χ1+χ2)/2),χa->((χ1-χ2)/2)};
χsχaRule = {χ1->χa+χs,χ2->-χa+χs,χs->χs,χa->χa};

parameters  injection rule

Here, we define a function that takes the component mass (Subscript[m, 1])((Subscript[M, ⊙])),m2((Subscript[M, ⊙])), component spins (Subscript[χ, 1]),(Subscript[χ, 2]), luminosity distance (Subscript[d, L]), reference time (Subscript[t, c]), and reference phase (Subscript[ϕ, c]), and returns replacement rules for all of the various parameters we consider in this notebook.

InjectionParams[
m1Input_(*Msol*),m2Input_(*Msol*),
χ1Input_,χ2Input_, 
dLInput_(*Mpc*),tcInput_(*s*),ϕcInput_]:=
{m1->m1Input,
m2 ->m2Input ,
m -> (m  /.m1m2Rule/.{m1->m1Input,m2->m2Input}),
ℳ-> (ℳ  /.m1m2Rule/.{m1->m1Input,m2->m2Input}),
η -> (η /.m1m2Rule/.{m1->m1Input,m2->m2Input}),
q ->( q /.m1m2Rule/.{m1->m1Input,m2->m2Input}),
δ-> (δ/.m1m2Rule/.{m1->m1Input,m2->m2Input}),
χ1->χ1Input,
χ2->χ2Input,
χs->(χ1Input+χ2Input)/2,
χa->(χ1Input-χ2Input)/2,
dL-> dLInput ,
tc -> tcInput ,
ϕc -> ϕcInput}

assumptions = ℳ>0&&q>0&&q<1&&m>0&&η<1/4&&m1>0&&m2>0&&f>0;

Taylor F2 waveform

This is the TaylorF2 waveform, which is the standard PN inspiral waveform (see arXiv:0907.0700)
This is the waveform we computed at Newtonian order in L03.

unit conversions

We want the waveform to take masses in terms of  (Subscript[M, ⊙]) and lengths in terms of Mpc, so we need to include the appropriate factors so that [(OverTilde[h])] = s

GVal =667430×(10^(-25)) (*(km^3) (kg^(-1)) (s^(-2))*);
cVal = 299792458 (10^(-3))(*km (s^(-1))*);

MsolSEC=1988416×(10^24) GVal (cVal^(-3)) (*sec*);
MpcSEC=30856775814913673 (10^3)(cVal^(-1))(*sec*);

post-Newtonian factor

The dimensionless PN velocity

v = (((m MsolSEC π f))^(1/3));

waveform amplitude

The TF2 amplitude contains only the Newtonian term

𝒜0 =(Sqrt[((5π η)/6)]) ((((m MsolSEC))^2)/(dL MpcSEC))(v^(-7/2));

waveform phase

The phase contains terms through 3.5PN (proportional to  (v^7)relative to leading order)

𝒫0 = (3/(128 ))(1/η)(v^(-5));

𝒫[0]=1;
𝒫[1]=0;
𝒫[2]=(20/9)((743/336)+(11/4)η);
𝒫[3]=-16π  +4 (((113/12)-(19/3)η)χs+(113/12)δ χa); 
𝒫[4]=10((3058673/1016064)+(5429/1008)η+(617/144)(η^2)) -10 (η((721/48)(((χs)^2)-((χa)^2))-(247/48)(((χs)^2)-((χa)^2)))+(1-2η)((719/96)(((χs)^2)+((χa)^2))-(233/96)(((χs)^2)+((χa)^2)))+δ((719/48)(χs χa)-(233/48)χa χs)); 
𝒫[5]=π((38645/756)-(65/9)η)(1+3 Log[v])-(((732985/2268)-(24260/81)η-(340/9)(η^2))χs+((732985/2268)+(140/9)η)δ χa )(1+3 Log[v]) ;𝒫[6]=(((11583231236531/4694215680)-((640(π^2))/3)-(6848/21)EulerGamma-(6848/21)Log[4v])+η(-(15737765635/3048192)+((2255(π^2))/12))+(76055/1728)(η^2)-(127825/1296)(η^3)) +π((2270/3)δ χa+((2270/3)-520η)χs);
𝒫[7]=π((77096675/254016)+(378515/1512)η-(74045/756)(η^2))+(-(25150083775/3048192)+(26804935/6048)η-(1985/48)(η^2))δ χa+(-(25150083775/3048192)+(10566655595/762048)η-(1042165/3024)(η^2)+(5345/36)(η^3))χs ;

We will want to be able to generate the waveform to various PN orders

ϕ[order_]:=-(π/4)+ 2 π f tc-2 ϕc +𝒫0 Sum[𝒫[n](v^n),{n,0,order}];

full waveform construction

We are disregarding the dependance on  extrinsic parameters which enter into the amplitude of the waveform . The factor below was chosen to obtain the correct SNR based for extrinsic parameters (Subscript[θ, jn]) = 0.4, ra = 1.375, dec = -1.2108 ,  ψ = 2.659.

ExtrinsicParamFactor  =0.2692959228694876;

Here, we construct the frequency domain strain and its conjugate 
We can choose the PN order of the phase by entering the highest power of v (i.e. 7 for 3.5PN)

FDStrain[order_]:=ExtrinsicParamFactor 𝒜0 (E^(-I ϕ[order]))/.ℳqRule/.χ1χ2Rule;
FDStrainConj[order_]:= ExtrinsicParamFactor 𝒜0 (E^(I ϕ[order]))/.ℳqRule/.χ1χ2Rule;

PSD

We need the PSD to compute inner products

Load the advanced LIGO O4 high power ASD . This is the default ASD Bilby uses for H1 and L1 . 
(you can also find it in the bilby package director: “python/site-packages/bilby/gw/detector/noise_curves/aLIGO_O4_high_asd.txt”, although I was having a hard time getting the latest Mathematica version to import the .txt file)

ASDTable=Quiet[Import[ "L07_FIM_PCA/aLIGO_O4_high_asd.csv"]];

Square it to obtain the PSD, express the PSD and frequencies in terms of km .

frequencies =ASDTable[[;;,1]] ;
psd = (((ASDTable[[;;,2]]))^2)  ;
PSDTable=Transpose[{frequencies,psd}];

Compute PSD interpolating function

PSD = Interpolation[PSDTable][f];

The inner product is

((Subscript[h, 1])|(Subscript[h, 2])) = 4 Re ∫(( (Subscript[(OverTilde[h]), 1])(SuperStar[(Subscript[(OverTilde[h]), 2])]))/(Subscript[S, n]))df

We use the standard minimum frequency 20 Hz and because we are using an inspiral only waveform, we cut off the integration at (Subscript[f, ISCO])

fmin = 20 ;
fISCO[m_]:=4400/m ; (*ISCO, Hz*)

InnerProd[h1Conj_,h2_,fmax_]:=4 Quiet[Re[NIntegrate[(h1Conj h2)/PSD,{f,fmin,fmax}]]]

Fisher matrix

Finally, we define the FIM 
The function takes: 
 - the frequency domain strain (OverTilde[h]) and its conjugate
 - a list of the parameters for which we want to compute the FIM
 - a list of those parameter' s values (in the same order as the previous list)
 - and the maximum frequency value we want to integrate to

FIM[h_,hconj_, Θ_,Θvals_,fmax_]:=
(ΘRule = Table[Θ[[aa]]->ΘVals[[aa]],{aa,1,Length[Θ]}];
(*computes {(∂_i)(Subscript[h, AP])((Subscript[θ, tr])), (∂_j)(SuperStar[(Subscript[h, AP])])((Subscript[θ, tr]))}*)
(*to speed things up, first sub in all param values except the one taking the derivative*)
(*also only computes lower corner, we symmetrize later*)
Γpre1=Table[If[aa>=bb,{D[h(*/.Drop[ΘRule,{aa}]*),Θ[[aa]]],D[hconj(*/.Drop[ΘRule,{bb}]*),Θ[[bb]]]},{0,0}],{aa,1,Length[Θ]},{bb,1,Length[Θ]}]/.ΘRule;
(*computes the inner products*)
Γpre2=Table[Length[Permutations[{aa,bb}]]InnerProd[Γpre1[[aa,bb,1]],Γpre1[[aa,bb,2]],fmax],{aa,1,Length[Θ]},{bb,1,Length[Θ]}];
(*symmetrize*)
Symmetrize[Γpre2]//Normal)

compare to Bayesian parameter estimation

We are going to compare the approximate posterior obtained through a Fisher analysis to the sampling of the true posterior obtained through Bayesian parameter  estimation. 

We will consider a non-spinning system with (Subscript[m, 1]) = 15 (Subscript[M, ⊙ ])and  (Subscript[m, 2]) = 5 (Subscript[M, ⊙ ]), corresponding to ℳ ~ 7.33  (Subscript[M, ⊙ ]) and q = (1/3), and an SNR of approximately 200. We will produce the signal with theTaylorF2 waveform into a single LIGO detector.  

To save time, we will only sample on two parameters, ℳ and q, which is the standard choice for binary parameter estimation. 

The posteriors used for comparison in this notebook were generated with L07_FIM_PC/PE.ipynb and are saved in L07_FIM_PCA/samples_chirp_mass_mass_ratio_q3.csv and “L07_FIM_PCA/samples_chirp_mass_mass_ratio_q1.csv

Define the replacement rule for the injection values

params= InjectionParams[10,10,0,0,40,0,0]

Because we are working with an inspiral only waveform, we will cut off integration at the ISCO

fmax = fISCO[20]

Define the sampled parameters and their injection values

Θ= {ℳ,q};
ΘVals= {ℳ,q}/.params;

Define a rule to set the non sampled parameters to their injection values

nonSampledParams = {χ1,χ2,tc,ϕc,dL};
nonSampledParamsVals =  {χ1,χ2,tc,ϕc,dL}/.params;
nonSampledParamsRule = Table[nonSampledParams[[aa]]->nonSampledParamsVals[[aa]],{aa,1,Length[nonSampledParams]}]

Set up the strain and its conjugate so that they are only functions of q, ℳ, and f

h =Simplify[FDStrain[7]/.nonSampledParamsRule, assumptions];
hconj = Simplify[FDStrainConj[7]/.nonSampledParamsRule, assumptions];

Compute the FIM

fim = FIM[h ,hconj, Θ,ΘVals,fmax]

Invert FIM to obtain the approximate covariance matrix of the posterior

Σ = Inverse[fim]

and compute the standard deviations

{Sqrt[Σ[[1,1]]],Sqrt[Σ[[2,2]]]}

Load the samples of the true posterior obtained through Bayesian analysis and compute the standard deviations

(* samples for q = 1/3*)
posteriorTrue = Import["L07_FIM_PCA/samples_chirp_mass_mass_ratio_q1.csv"];
(* samples for q = 1*)
(*posteriorTrue = Import["L07_FIM_PCA/samples_chirp_mass_mass_ratio_q1.csv"];*)

and compute the standard deviations

((Variance[posteriorTrue])^(1/2))

Produce an approximate posterior by sampling from a multivariate Gaussian distribution with the covariance given by the inverse of the FIM centered at the injection values

posteriorFisher = RandomVariate[MultinormalDistribution[ΘVals,Symmetrize[Σ]],Length[posteriorTrue]];

ListPlot[{posteriorTrue(*,posteriorFisher*)},
PlotStyle->{Orange,Purple},
Frame->True,
FrameLabel->{"ℳ ((Subscript[M, ⊙ ]))","q" },AspectRatio->1,
PlotLegends->{"True Posterior","Fisher Approximation"}]

Generally, we say a FIM has done a good job approximating the posterior when the standard deviations are within an order of magnitude of the true values, so we have done a pretty good job here. What factors may be contributing to the differences between the true posterior and the approximation?

What happens if you set (Subscript[m, 1]) = 10(Subscript[M, ⊙ ])and  (Subscript[m, 2]) = 10 (Subscript[M, ⊙ ]) so that q = 1?  What most likely causes the breakdown?

principal component analysis

Parameter estimation can be tricky when parameters are strongly correlated. We can use the FIM to identify combinations of the model parameters that are less correlated than the original parameters. 

See arXiv:1304.7017 for an overview.

Multivariate Gaussians

Lets return to the multivariate Gaussian distribution

Σ = {{((σ1)^2),ρ σ1 σ2},{ρ σ1 σ2,((σ2)^2)}};
(*Σ = {{Σ11,Σ12},{Σ12,Σ22}};*)
μ = {μ1, μ2};
θ = {θ1,θ2};

pdf = PDF[MultinormalDistribution[μ,Σ],θ]

We will reparametrize it in terms of

(Subscript[φ, i ])= (Subscript[Λ, ij])(Subscript[θ, j])

Where (Subscript[Λ, ij])  is the (j^(th ))component of the (i^th) eigenvector of  (Σ^(-1))

EigenVecs = Eigenvectors[Inverse[Σ]];
EigenVals = Eigenvalues[Inverse[Σ]];

Invert the relationship to get θ in terms of φ

φ = {φ1, φ2};
Φexpr = EigenVecs.θ//Simplify
θtoφRule = Solve[φ == Φexpr,θ]//Flatten//Simplify

And write the PDF in terms of φ

pdfNew = pdf/.θtoφRule

We can determine  (Σ^(-1)) of the distribution on φ

ΣInvNew = -Table[D[Log[pdfNew],φ[[aa]],φ[[bb]]],{aa,1,2},{bb,1,2}]//Simplify;
%//MatrixForm

The off-diagonals are zero! (Subscript[φ, 1]) and (Subscript[φ, 2]) are not correlated!

PCA of the (post-)Newtonian Waveform

When computing the waveform, it is natural to parameterize the masses in terms of total mass m and symmetric mass ration η. However this is not the best parameterization for parameter estimation.

Define the replacement rule for the injection values

params= InjectionParams[12,8,0,0,40,0,0];
fmax = fISCO[m/.params];

Define the sampled parameters and their injection values.
A PCA returns linear combinations of the original parameters.
We would like our final result to be the product of the original parameters so we will work in terms of logs.

Θ= {logm,logη};
ΘVals= {Log[m],Log[η]}/.params;

logRule = {m ->Exp[logm],η ->Exp[logη]};

Define a rule to set the non sampled parameters to their injection values

nonSampledParams = {χ1,χ2,tc,ϕc,dL};
nonSampledParamsVals =  {χ1,χ2,tc,ϕc,dL}/.params;
nonSampledParamsRule = Table[nonSampledParams[[aa]]->nonSampledParamsVals[[aa]],{aa,1,Length[nonSampledParams]}]

Set up the strain and its conjugate so that they are only functions of log[m], log[η], and f 
We will begin with the Newtonian waveform!

h =  Simplify[FDStrain[2]/.mηRule,assumptions]/.logRule/.nonSampledParamsRule
hconj = Simplify[FDStrainConj[2]/.mηRule,assumptions]/.logRule/.nonSampledParamsRule

Compute the FIM

fim = FIM[h ,hconj, Θ,ΘVals,fmax]

Compute the eigen vectors of the FIM and rescale them so we can express them in terms of rational numbers

Eigenvectors[fim];
EigenVecs =%/(%[[1,1]])//Rationalize

Our new parameters are

Φexpr = EigenVecs.Θ

The first parameter should look familiar!

ℳ/.mηRule

We know that the Newtonian waveform only depends on ℳ, so it makes sense it is a principal component

Simplify[FDStrain[0],assumptions]

Invert the relationship and compute the values of log[φ1]  and log[φ2]

Φ = {logφ1, logφ2};
θtoφRule = Solve[Φ==Φexpr,Θ]//Flatten
ΦVals = Φexpr/.(Table[Θ[[aa]]->ΘVals[[aa]],{aa,1,Length[Θ]}])

Express the waveform in terms of log[φ1]  and log[φ2]
As we should expect, for the Newtonian form, log[φ2] drops out of the expression

hNew = Simplify[h/.θtoφRule, logφ1>1&& logφ2>2]
hconjNew = Simplify[hconj/.θtoφRule, logφ1>1&& logφ2>2]

It follows that the new FIM is diagonal (it only has one element)

fimNew= FIM[hNew ,hconjNew,Φ,ΦVals,fmax]

Inverse[fimNew]

What happens when you add PN corrections to the waveform?
Why is the FIM not always diagonal even when the waveform is expressed in terms of the principal components?
What happens when you compute the PCA in a different region of parameter space?

A short coming of PCAs

What happens if we start with the waveform in terms of different mass parameters, such as (Subscript[m, 1]) and (Subscript[m, 2])?

Θ2= {logm1,logm2};
ΘVals2= {Log[m1],Log[m2]}/.params;

logRule2 = {m1 ->Exp[logm1],m2->Exp[logm2]};

h2 =  Simplify[FDStrain[0]/.m1m2Rule,assumptions]/.logRule2/.nonSampledParamsRule;
hconj2 = Simplify[FDStrainConj[0]/.m1m2Rule,assumptions]/.logRule2/.nonSampledParamsRule;

fim2 = FIM[h2 ,hconj2, Θ2,ΘVals2,fmax]

Eigenvectors[fim2];
EigenVecs2 =%/(%[[1,1]])//Rationalize;
Φexpr2 = EigenVecs2.Θ2

A PCA only returns linear combinations of the original parameters. There is no way to express log[ℳ] as a linear combination of  log[ (Subscript[m, 1])] and log[(Subscript[m, 2])] so the PCa cannot return log[ℳ] even though ℳ is clearly the best parameter to sample in

Estimating biases

We can also use FIM to estimate the bias in the maximum likelihood parameters due to noise and mismodeling. 

See arXiv:0707.2982 and arXiv:2410.06254