Fallo en la estimación de BSR

[imendez]

Hola, tengo un problema con BSR, a ver si me podéis ayudar a encontrar el fallo.
Tengo una estimación con 240 parámetros lineales, 37 hiperparámetros, 3 parámetros ARMA y 28 omitidos:

[BysMcmc::DefineBlock] Block (1) BSR.MainLinBlock.Modelo_CEN02 of size 277
[BysMcmc::DefineBlock] Block (2) BSR.SigmaBlock.Modelo_CEN02 of size 1
[BSR.ArimaBlock.Modelo_CEN02]  Initializing block
  [BSR.ArimaBlock.Modelo_CEN02]  Original data length = 3726
  [BSR.ArimaBlock.Modelo_CEN02]  Differenced data length = 3719
  [BSR.ArimaBlock.Modelo_CEN02]  ARMA parameters = 3
[BysMcmc::DefineBlock] Block (3) BSR.ArimaBlock.Modelo_CEN02 of size 3
[BysMcmc::DefineBlock] Block (4) BSR.InputMissingBlock.Modelo_CEN02 of size 0
[BysMcmc::DefineBlock] Block (5) BSR.OutputMissingBlock.Modelo_CEN02 of size 28

Me da el siguiente error:

CHOLMOD warning: matrix not positive definite
Warning: [372] [CHOLMOD 1] at line c:\users\jsperez\suitesparse\cholmod\supernodal\t_cholmod_super_numeric.c:614: matrix not positive definite
ERROR: [1] No es posible aplicar CholeskiFactor a una matriz virtual no definida positiva  Cholmod.R.Sparse(277x277)

Por el error creo que se debe a que hay una combinación lineal, pero no hay ningún Warning previo que indique eso. Es más, si estimo sin jerarquías todo funciona aparentemente bien, aunque da el siguiente Warning que no sé si tiene algo que ver:

Simulating iterations   201-  300 [done  18.27%] [time/sim: 0.135625 s] [remaining    122 s]
Warning: [372] [BSR.MainLinBlock.Modelo_CEN02] Cholesky solving had rounding error in interval [5.784e-012,1.432454155292362e-011 for 22 times of 300
 Simulating iterations   301-  400 [done  27.36%] [time/sim: 0.134219 s] [remaining    107 s]

Si a alguien se le ocurre alguna pista de qué puedo estar haciendo me vendrá muy bien, porque estoy bastante perdido.
Adjunto dos archivos por si os ayudan: el ".bsr" para poder reproducir la estimación, y el ".oza" con la matriz fallida (Failed factorization matrix).

Los parámetros de conifguración de BSR utilizados son:

NameBlock bsr.config =
[[
//MCMC dimensions
  Real mcmc.burnin       =  100;//100;//2000
  Real mcmc.sampleLength =  1000;//500;//10000
  Real mcmc.cacheLength  =  100;
//Basic master configuration
  Real basic.cholesky.epsilon          = 1.E-13;
  Real basic.cholesky.warningFreq      = 100;
  Real basic.truncMNormal.gibbsNumIter = 5;
//Report configuration
  Real report.raftery.diag.q   = 0.025;
  Real report.raftery.diag.r   = 0.007;   
  Real report.raftery.diag.s   = 0.950;   
  Real report.raftery.diag.eps = 0.001; 
  Real report.acf.lag          = 20,
  Real report.histogram.parts  = 100;
  Real report.kerDens.points   = 0;
  Real report.kerDens.numIter  = 2;
  Real report.kerDens.epsilon  = 0.001;
//If it is true the automatic generated input missing variables will use 
//id_node as prefix of identifier.
  Real DBApi.useNodeInMissingIdentifier = True;
//Use Metropolis-Hastings instead og Gibbs for linear blocks after
//this number of iterations
  Real bsr.linBlk.MH.useAfterIter   = 10;
//Use Gibbs method after these rejections of Metropolis-Hastings
  Real bsr.linBlk.MH.maxRejected    = 2;
//Use Gibbs method after these consecutive iterations of Metropolis-Hastings
  Real bsr.linBlk.MH.maxConsecutive = 10; 
//Generic flags
  Real do.resume               = False;
  Real do.report               = True;
  Real do.eval                 = True;
  Real do.linear.effects       = True;
  Real do.save.DBApi.Estim.Oza = False;
//Metodo de aproximacion del arima
  Text bsr.arimaFilter   = BysMcmc::Options::Arima.Filter::FastCholSea;
  Real bsr.arimaSkipIter = 10//mcmc.burnin + mcmc.sampleLength + 1
]];

Gracias de antemano.

[Víctor de Buen Remiro]

Lo primero que se debe hacer en estos casos es determinar de qué tipo de singularidad se trata. En los casos como este en los que no hay demasiadas variables es posible usar la SVD para ver si hay autovalores nulos o casi nulos y cuántos hay

//Recupera la matriz que ha dado el problema
Matrix X = {VMat2Mat(Include(
"C:/Documents and Settings/vdebuen.FROBENIUS/Datos de programa/"
"tol/tmp/_BSR.MainLinBlock.Modelo_CEN02___error_factorize.oza")[1])};
//Realiza la descomposición de valor singular
Set svd = SVD(X);
//Toma la diagonal principal que contiene a los autovalores
Matrix D = Tra(SubDiag(svd[2],0));

Esto da en este caso nada más y nada menos que 16 autovalores casi nulos, por debajo de {{ 6.4624638401e-015 }}} que es prácticamente cero en la aritmética binaria usual, y llega hasta el mínimo de 5.85564053314e-113 . No se tratapues de un problema numérico transitorio o algo evitable algorítmicamente, sino de un diseño totalmente inviable.

Ahora contruimos la matriz del núcleo de X extrayendo las 16 últimas columnas de la matriz de autovectores

  VMatrix kernel = Mat2VMat(SubCol(svd[3],Range(261,277,1)),False,0.5,1.E-14);

Aunque nosiempre hay tanta suerte, en este caso la matriz es bastante dispersa y tiene los coeficientes distribuidos en cajas por lo que es facil deducir las variables colineales. Por ejemplo, las tres siguientes variables lo son

[241] Hyp.VeranoSemanasCasos.VentaEmu.CAMCAPCEN02.Ver_IndSem_00_NotSol,
[254] Hyp.VeranoSemanasCasos.VentaEmu.CAMCAPCEN02.Ver_IndSem_16_NotSol
[256] Hyp.VeranoSemanasCasos.VentaEmu.CAMCAPCEN02.Ver_IndSem_17_NotSol

Los grupos de variables colineales son

[[241,254,256]],
[[261,262,263]],
[[264,265,266]],
[[267,268]],
[[272,273]],
[[274,275]]

La matriz de datos originales parece estar bien formada, es decir no es colineal, pero la matriz que le lllgaal simulador del bloque lineal tiene todo ceros en esas variables. Al mirar en la definición del modelo se observa que todas ellas corresponden a ecuaciones de nodos latentes con varianza cero, lo cual es absurdo

Lat.VeranoSemanasCasos.VentaEmu.CAMCAPCEN02.Ver_IndSem_00_Sol::Noise[7] ~ Normal(0,0);

Esto es algo que el parser debería detectar para avisar de que hay un error.

[Víctor de Buen Remiro]

(In [1779]) Fixes #798

[Víctor de Buen Remiro]

(In [1780]) Fixes #798