Hi,

I am trying to sample spherical harmonics.

Initially I thought of the model coefficients as independent and Gaussian distributed. I was successful in sampling to degree 25 (675 dimensions).

I have now been trying to think of the coefficients as dependent by using some simulations (from someone else) to define a covariance matrix for the coefficients. My plan was then to use a singular value decomposition of the covariance matrix to project to and from eigen-space. In that way I could sample independent coefficients and project them back to the original space.

But when trying my new prior I cannot get above degree 4 (24 dimensions), when using similar data and sampling inputs, before reaching Rhat < 0.9 and Rhat > 1.1.

Can anyone see if there is anything wrong with my model?

I have included the models below and attached data as rdump.

Covariance model:

```
data {
int<lower=0> L; // Length of data-array
int<lower=0> n; // SH-degree
int<lower=0> mMax; // Amount of model parameters
vector[L] d; // Declare size of data-array
matrix[L, mMax] G; // Declare size of design matrix
vector[mMax] orgMean; // Mean of original model data
vector[mMax] orgStd; // Standard deviation of original model data
matrix[mMax, mMax] T; // Transfer matrix
}
parameters {
vector[mMax] q; // Declare model array (eigen-space)
}
transformed parameters {
vector[mMax] m;
m = (T*q+orgMean).*orgStd; // Project back to original space
// q in projected to standardized orginal space
// Then it is scaled using orgMean and orgStd
}
model {
vector[L] res; // Declare residual array
// Prior distribution
for (i in 1:mMax) {
q[i] ~ normal(0, 1); // Assume normal distribution
}
// Likelihood
res = d-G*m;
res ~ normal(0,1);
}
```

Original model assuming independent coefficients:

```
data {
int<lower=0> L; // Length of data-array
int<lower=0> n; // SH-degree
int<lower=0> mMax; // Amount of model parameters
vector[L] d; // Declare size of data-array
matrix[L, mMax] G; // Declare size of design matrix
vector[mMax] limits;
}
parameters {
vector[mMax] m; // Declare model array
}
model {
vector[L] res; // Declare residual array
// Prior distribution
for (i in 1:mMax) {
m[i] ~ normal(0, limits[i]/2);
}
// Likelihood
res <- d-G*m;
res ~ normal(0,1);
}
```