I’m trying to implement the Bayesian hierarchical correlated t-test from Corani et. al. (see: https://link.springer.com/content/pdf/10.1007%2Fs10994-017-5641-9.pdf) in which they build a generative model describing the generation of an array of performance measures for two classifiers across multiple datasets.

In short, the model looks as follows:

Given an array of observations of shape (n_datasets, n_obs), where each column holds measurement of performance for a given train/test split, and each row refers to the observations for a given dataset, the generative model is as follows:

Assume 3 datasets and 5 observations of performance, i.e. `arr.shape=(3, 5)`

.

For each dataset (i.e. each row) we have;

```
x ~ MVN(mu, cov)
mu ~ t(mu_0, sigma_0, nu)
sigma ~ unif(0, 100),
```

Where `x`

is the vector of observations (i.e. 5 dimensional) for the given dataset which are jointly multivariate-normal distributed with the same mean, `mu`

, and the same variance, `sigma^{2}`

, and the same correlation, `p`

.

More specifically the mean is a vector of size 5 whose elements are all equal to `mu`

and the covariance matrix is defined to have diagonal elements `sigma^{2}`

and out-of-diagonal elements `p * sigma^{2}`

, with the correlation, `p`

, fixed a priori.

My question is: how do I define such a covariance matrix, i.e. one based on a random variable sampled from a uniform distribution and the given correlation? I would like to be able to sample the sigma and then create the covariance matrix which will then be fed to the multivariate likelihood. Would a good way to do this be:

cov = math.dot(np.eye(n_obs), math.sqr(sigma))) +

math.dot(np.ones((n_obs, n_obs)) - np.eye(n_obs), p * math.sqr(sigma))

This is my first question on this forum so apologies if I’m not formulating the question in the correct manner, happy to amend if required.

Thanks!