singular matrix

[MaggieLieu]

Hi,

I have a sporadic problem that occasionally I will run into the following error:

Error in qr.solve(qx, Y.tilde[1:n + (j - 1) * n, 1]): singular matrix 'a' in solve
Traceback:

1. Gibbs.regression(x.in = x_obs, y.in = y_obs, M = cov_mat, Nsamples = nmc)
2. qr.solve(qx, Y.tilde[1:n + (j - 1) * n, 1])
3. stop("singular matrix 'a' in solve")

Seems to occur more often when I'm using the hierarchical implementation, and it can go away as long as I choose a really small sample size.

This is what I am trying to do:
nx = 100
x = runif(nx, 13,15)
alpha = 3.5
beta = 1.5
y = beta*x + alpha
dx = 1.4
dy = 2
dint = 1 #intrinsic scatter in y
x_obs = x + rnorm(nx, 0, dx)
y_obs = y + rnorm(nx, 0, dy) + rnorm(nx, 0, dint)
cov_mat = array(0, dim=c(2,2,nx))
cov_mat[1,1,] <- dx^2
cov_mat[2,2,] <- dy^2
nmc = 5000 #number of samples
post = Gibbs.regression(x.in = x_obs, y.in = y_obs, M = cov_mat, Nsamples = nmc) #fit line

[abmantz]

I think this can only occur in the same circumstances and for the same reason as you might get a singular matrix in ordinary least squares, i.e. a huge degeneracy between slope and intercept parameters. Certainly the only place that qr.solve appears in the source code is when sampling the slope and intercept parameters, which is (conditionally) an OLS calculation.

Centering the covariates by subtracting 14 from x_obs appears to mostly eliminate the problem. There is still a strong degeneracy between the intercept and slope, presumably because the constraining power of the data set is so small, and some chains eventually did crash as a result, but mostly they stayed stable. This is probably the most one can expect when the priors are uniform and the measurement error in x is wider than the actual range.

I couldn't guess why using a small sample size would help (in the uncentered case), unless it's simply that having fewer data keeps the parameter correlation from going numerically unstable by inflating the overall uncertainties.

[MaggieLieu]

Yes, pivoting the x and y values does seem to help on the standard model, but still many of the errors when using the Dirichlet process which can only really be used with small errors, probably will go away with tighter priors though.

[abmantz]

Can you clarify whether the "many" errors are the same one, or are these new errors?

It seems to me that the same failure mode is likely to exist if the x measurement errors are large compared with their range. The Dirichlet processes should then be consistent with a single cluster, meaning that internally all true x values will be identical and a slope will be completely unconstrained by the data.