Firstly, I want to give you a maximum of informations and precisions about my issue. If I can't manage to get the expected results, I will launch a bounty, maybe some experts or symply people who have been already faced to a similar problem will be able to help me
I have 2 covariance matrices known $\text{Cov}_X$ and $\text{Cov}_Y$ that I want to cross-correlate.
1) For this, I have performed a diagonalisation of each Fisher matrix $F_X$ (FISH_eigen_sp_flat in code) and $F_Y$ (FISH_eigen_xc_flat in code) associated of Covariance matrices $\text{Cov}_X$ and $\text{Cov}_Y$.
So, I have 2 different linear combinations that are uncorraleted, i.e just related by eigen values ($\dfrac{1}{\sigma_i^2}$) as respect of their combination.
Then, I get the diagonal matrices $D_X$ and $D_Y$. I can't build a "global" Fisher matrix directly by summing the 2 diagonal matrices since the linear combination of random variables is different between the 2 Fisher matrices.I have eigen vectors represented by $X$ (with $D_X$ diagonal) and $Y$ matrices (with $D_Y$ diagonal matrix)
That's why I think that I could perform a "global" combination of eigen vectors where I can respect the MLE estimator as eigen value :
$$\dfrac{1}{\sigma_{\hat{\tau}}^{2}}=\dfrac{1}{\sigma_1^2}+\dfrac{1}{\sigma_2^2}\quad(1)$$
because $\sigma_{\hat{\tau}}$ corresponds to the best estimator from MLE method (The only inconvenient in this approach for me is that I respect it only on a linear combination of random variables, and not just one, if someone could tell me if am wrong or not about this)
So, I thought a convenient linear combination that could allow to achieve it would be under the form :
$$M = X + Y + Z$$
with $$Z =aX + bY$$
where $a$ and $b$ have to be determined and such that :
$$\text{Var}(X+Y+Z)=\text{Var}(X)+\text{Var}(Y)$$
Therefore, I am looking for a way to find $Z$ (not null) that could verify :
$$\text{Var}(X+Y+Z)=\text{Var}(X+Y)+\text{Var}(Z)+2\,\text{Cov}(X+Y,Z)$$
$$= \text{Var}(X)+\text{Var}(Y)+2\,\text{Cov}(X,Y)+\text{Var}(Z)+2\,\text{Cov}(X,Z)+2\,\text{Cov}(Y,Z)=\text{Var}(X)+\text{Var}(Y)$$
So, we can remove on both sides the terms $\text{Var}(X)+\text{Var}(Y)$ to get the equation :
$$2\,\text{Cov}(X,Y)+\text{Var}(Z)+2\,\text{Cov}(X,Z)+2\,\text{Cov}(Y,Z)=0$$
It follows with solution under the form $Z=a\,X+b\,Y$ :
$$2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+b^2\,\text{Var}(Y)+2ab\,\text{Cov}(X,Y)+2\,\text{Cov}(X,aX+bY)+2\,\text{Cov}(Y,aX+bY)$$
$$=2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+b^2\,\text{Var}(Y)+2ab\,\text{Cov}(X,Y)+2\,\text{Cov}(X,aX)+2\,\text{Cov}(X,bY)$$
$$=2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+b^2\,\text{Var}(Y)+2ab\,\text{Cov}(X,Y)+2a\,\text{Var}(X)+2b\,\text{Cov}(X,Y)+2a\,\text{Cov}(X,Y)+2b\,\text{Var}(Y)=0\quad(2)$$
Finally, If I fix a value for $a$, I have to solve this equation like a second order equation with $b$ as unknown parameter (below b=x=unknown) :
$$\alpha\,x^2+\beta x + \delta$$
with :
$\alpha=\text{Var}(Y)$
$\beta=2\,\text{Cov}(X,Y)+2a\,\text{Cov}(X,Y)+2\,\text{Var}(Y)$
$\delta=2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+2a\,\text{Var}(X)+2a\,\text{Cov}(X,Y)$
To avoid complex solutions of this second order equation, I have to respect :
$$\beta^2-4\,\alpha\delta > 0$$
2) We could for example take $a=-1$ and try to find the b unknown like this (more precision, I think that b unknow parameter is surely a matrix or a vector : if someone could confirm it since I take the scalar value $a=1$, this would fine to tell it). This choice of $a=-1$ may seem to be arbitral, it is not only the unique solution for this issue but I have not any error concerning my issue.
Therefore, by taking $a=-1$, I have implemented a function to compute varX, varY and CovXY in Python :
# Solution Z = aX + bY ; Taking a = -1 => solve unknown b represented by b_coef :
b_coef = np.zeros((7,7))
VarX = np.zeros((7,7))
VarY = np.zeros((7,7))
CovXY = np.zeros((7,7))
# Compute covariance between vectors of matrix
def compute_Cov(A,B):
C = np.zeros((7,7))
for i in range(7):
C[0:7,i]= np.mean(A[0:7,i]*B[0:7,i]) - np.mean(A[0:7,i])*np.mean(B[0:7,i])
return C
I can this way compute the variance of X (compute_Cov(X,X)), Y (compute_Cov(Y,Y) and covariance(X,Y) (compute_Cov(X,Y)).
Diagonalisation) By diagonalising the 2 Covariance matrices (inverse of Fisher matrices FISH_sp_flat and FISH_xc_flat) representing the variances on each parameters that I want to compute the standard deviations, I think I build a linear combination of these parameters that is uncorrelated :
This way, I can sum respectively for each combination the eigen values like this :
# Get eigen values (marginalized error since we handle Covariance error)
# and eigen vectors giving matrix "P" with convention : F = P D P^-1
eigen_sp_flat, eigenv_sp_flat = np.linalg.eig(np.linalg.inv(F_X))
eigen_xc_flat, eigenv_xc_flat = np.linalg.eig(np.linalg.inv(F_Y))
# Sum of FISH_eigen sp and xc
FISH_eigen_sp_flat = np.linalg.inv(np.diag(eigen_sp_flat))
FISH_eigen_xc_flat = np.linalg.inv(np.diag(eigen_xc_flat))
# MLE method (Maximum Likelihood Estimator) : Sum of FISH_eigen sp and xc : 1/sig^2 = 1/sig1^2 + 1/sig^2
FISH_eigen_sum = FISH_eigen_sp_flat + FISH_eigen_xc_flat
3) Once linear combination of parameters and Var(X), Var(Y), Cov(X,Y) computed, I build my final covariance matrix like this :
# Sum of passing matrix X+Y+Z with X = eigenv_sp_flat, Y = eigenv_xc_flat and Z to infer
eigenv_final = eigenv_sp_flat + eigenv_xc_flat + a*eigenv_sp_flat + b_coef*eigenv_xc_flat
# Covariance and Fisher matrix final
COV_final = np.dot(eigenv_final, np.linalg.inv(FISH_eigen_sum))
COV_final = np.dot(COV_final, np.linalg.inv(eigenv_final))
FISH_final = np.linalg.inv(COV_final)
The issue is that I get a better FoM (1389) (Figure of Merit = 1/area(w0,wa) which is the inverse area of 1 C.L contour for the joint parameters (w0,wa)) compared to a simple sum betwen the 2 starting Fisher matrix (1235): I would expect a more significant improvement and some anomalies appears also in the fnal constraints (for example, for a random variable, it doesn't respect the equation(1) ( $\dfrac{1}{\sigma_{\hat{\tau}}^{2}}=\dfrac{1}{\sigma_1^2}+\dfrac{1}{\sigma_2^2}\quad(1)$), i.e I have no gain for $\sigma_{\hat{\tau}}^{2}$ compared to the smallest individual $\sigma$ ($\sigma_1$ or $\sigma_2$).
Anyone could confirm me if my approach is correct, especially the calculation of a new basis that cheks V(X+Y+Z) = Var(X)+Var(Y) which seems for me essential ?
I hope that I have been enough clear in my explanations : the issue is simple and a little tricky in the same time.
EDIT 1: I realized that the condition of orthogonality ($F.T^T\neq I_d$) on the final building of eigen vectors is not respected in my code.
So, this characteristic of orthogonality is essential if I want to check the condition $\text{Var}(X+Y+Z)=\text{Var}(X)+\text{Var}(Y)$ where $Z =ax + bY$.
Therefore, I think this condition of orthogonality adds a complementary condition on the choice (and so the computation) of coefficient $a$ and $b$ into $Z =ax + bY$.
That's why I have asked another question about this criterion of orthogonality on : Building an orthogonal basis from the sum of different orthogonal basis.
Hope this will help you
Anyone could see how to translate this condition to determine a single value for $a$ and for $b$ ?
Don't hesitate to ask further information.
Any help/fix/suggestion is welcome. Sorry if this is a little bit long to read.
EDIT 2 : I have slightly modified the function compute_Cov : Does it make sense to compute variance $\text{Var}(X)$, $\text{Var}(Y)$, and $\text{CoVar}(X,Y)$ with $X$, $Y$ covariance matrices like this ? :
# Compute covariance between vectors of matrix
def compute_Cov(A, B):
C = np.zeros(7)
for i in range(7):
C[i]= np.cov(A[0:7,i]*B[0:7,i])
return C
VarX = compute_Cov(X,X)
VarX = compute_Cov(Y,Y)
CoVarX = compute_Cov(X,Y)
I have serious doubts about what I do since the beginning, anyone could see clearer ?
You are trying to find $u, v, \text{and }w$ such that:
$$ u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y = 0 \\ u \geq 0 \\ -1 \leq v, w \leq +1, $$
where $u = \sigma_Z$, $v = \rho_{X,Z}$ and $w = \rho_{Y,Z}$. Once you know $u, v, \text{and }w$, it is straightforward to find $a \text{ and } b$ in $Z = aX + bY$.
This needs to be solved numerically. You might have to preclude edge cases ($\rho_{XY} = \pm 1$).
An example: $\sigma_1 = \sigma_X = 0.1, \sigma_1 = \sigma_Y = 0.2$ and $\rho_{XY}$ varies from -0.9 to 0.9. $a$ and $b$ from minimising $(u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y)^2$ subject to the constraints look like this. The optimal objective values are of the order of $10^{-15}$, so practically $0$.