Estimated covariance matrix

Pymc3 how to calculate the covariance matrix of the sample and compare it with the real covariance matrix to get the mean square error.

there is a pm.trace_cov function for that.

Thank you for your help. I have another question: how to estimate the correlation matrix?

sd = np.sqrt(np.diag(cov))
corr = np.diag(sd**-1).dot(cov.dot(np.diag(sd**-1)))

Thank you,it works,
now I want to train my model on multiple data sets, and calculate the mean square error of covariance matrix on each data set. Finally, I will take the average after adding these results together.so I wrote the following code

for i in range(M):   # M data sets
    with pm.Model() as model:
        packed_l=pm.LKJCholeskyCov('packed_l',n=2,eta=2,sd_dist=pm.HalfCauchy.dist(2.5))
        L=pm.expand_packed_triangular(2,packed_l)
        Sigma=pm.Deterministic('Sigma',tt.dot(L,L.T))
        mu=pm.Normal('mu',0.,10,shape=2,testval=dataSet[i].mean(axis=0))
        obs=pm.MvNormal('obs',mu=mu,chol=L,observed=dataSet[i])
        trace = pm.sample(random_seed=seed,model=model,njobs=1)
    
    Sigma_post=trace['Sigma'].mean(axis=0)  # Estimated covariance matrix

    MSE0=[]
    MSE0.append(1/(k**2)*np.sqrt(((Sigma_post-Sigma_real)**2).sum())) # mean square error
MSE=np.array(MSE0)
RMSE=np.sqrt((1/M)*MSE.sum())# Root mean square error
print(MSE.shape)
print(RMSE)

It should be MSE.shape=(1,M), but the run results indicate that MSE.shape=(1,),

I am very confused. Where is the mistake? Can you give me some advice?

You have

    MSE0=[]

right above

    MSE0.append(...)

which clears the MSE0 every time. You should put MSE0=[] outside of the for-loop.