La función SliceSampler1D falla con dominios inconexos

[Víctor de Buen Remiro]

Tal y como se puede ver en el siguiente código TOL, la función SliceSampler1D falla a veces con dominios inconexos, es decir, cuando en parte del intervalo la densidad vale 0, o lo que es lo mismo su logaritmo es -Inf.

Como no sucede de forma determinista hay que usar una semilla de geración distinta para reproducir cada caso posible.

Por lo que he podido ver, devuelve el último valor generado, por lo que si éste tiene densidad 0 el valor devuelto está fuera del dominio.

El ejemplo reproduce la simulación del coeficiente de grado 1 de un modelos AR puro de grado 6, tal y como lo hace internamente BSR.

//////////////////////////////////////////////////////////////////////////////
// Semilla de generación aleatoria fijada para reproducir el comportamiento
// exacto
// 
// 3 de 8 no estacionarias incluida la última: devuelve el valor inicial
// Real PutRandomSeed(1836343167);
// 
// ninguna no estacionaria: devuelve una simulación válida
// Real PutRandomSeed(1664310743);
// 
// 2 de 8 no estacionarias, pero no la última: devuelve una simulación válida
// Real PutRandomSeed(1664310743);
//////////////////////////////////////////////////////////////////////////////
Real PutRandomSeed(0);
Real rndSeed = GetRandomSeed(0);

WriteLn("TRACE  rndSeed:"<<rndSeed);

//Se define un modelo AR con un polinomio que tiene un dominio
//inconexo de estacionariedad para el coeficiente de grado 1
Polyn ar = 1-0.09944049062890037*B+0.4208976407813203*B^2-0.07979027150930651*B^3-0.524509273364506*B^4+0.1794384784972413*B^5+0.01042162544590616*B^6;
Polyn ma = 1;

//Se copia el polinomio ma para evitar efectos secundarios en la evaluación
Polyn p = Copy(ar);
Real d = 1;

//Esta función evalúa la estacionariedad de un polinomio pol con respecto al 
//coeficiente de grado d en el intervalo [-d,d] con longitud de paso h
Matrix get.stationary.map(Polyn pol, Real d, Real h)
{
  SetMat(EvalSet(Range(-d,d,h), Set(Real c)
  {
    Polyn p = Copy(pol);
    Real c_ = If(c, c, 1.E-10);
    Real PutCoef(p,d,c_);
    [[c,IsStationary(p)]]
  }))
};
Matrix map = get.stationary.map(p,d,0.001);

//Se debe comenzar por un polinomio que sea estacionario
Real isStationary = IsStationary(ar);

//Se generan aleatoriamente los residuos
VMatrix res = Gaussian(1000,1,0,1);

//Para ilustrar mejor como el dominio de estacioanriedad no es conexo
//construimos series de ruidos en tres valores

Real c1 = -0.5; VMatrix z1 = { Real PutCoef(p,d,c1); DifEq(ma/p,res) };
Real c2 = +0.3; VMatrix z2 = { Real PutCoef(p,d,c2); DifEq(ma/p,res) };
Real c3 = +0.7; VMatrix z3 = { Real PutCoef(p,d,c3); DifEq(ma/p,res) };

//Para las pruebas nos sirve cualquiera estacionario: tomamos el original
//Retomamos el valor inicial
VMatrix z = DifEq(ma/ar,res);
Polyn p := Copy(ar);

//Se genera el evaluador de la verosimilitud ARMA
NameBlock almagro = ARMAProcess::Eval.Almagro(ar,ma,z,1);
//Se calculan los valores iniciales condicionados
VMatrix u0 = almagro::Draw.U_cond_Z(0);
VMatrix z0 = almagro::Get.Z0(u0);
VMatrix a0 = almagro::Get.A0(u0);

//Se toma el coeficiente de grado 1
Real x0 = Coef(ar,d);
Real L = -1;    //Mínimo valor para el que se da estacionariedad
Real U = 0.807; //Máximo valor para el que se da estacionariedad

//Función de actualización del polinomio AR que devuelve cierto si es 
//estacionario con el nuevo coeficiente
Real SetValue(Real x0)
{
  //Se asegura de que no pone coeficiete 0 porque destruye la forma del 
  //polinomio como vector de monomios. Es algo que habría que arreglar.
  Real x = If(1+x==1, 0.000001, x);
  Real PutCoef(p, d, x);
  Real isStationary = IsStationary(p);
  WriteLn("TRACE  SetValue("<<x+") isStationary:"<<isStationary);
  isStationary
};
//Función de evaluación del logaritmo de la densidad
Real LogDens.Almagro.Z_cond_U(Real x)
{
  Real llh =  If(!SetValue(x), -1/0, {
    Copy(almagro::LogLH.Z_cond_U(p,ma,z0,a0))
  }) ;
  WriteLn("TRACE  LogDens.Almagro.Z_cond_U("<<x+") :"<<llh);
  llh
};
//Llamada al método de simulación escalar slice
Real x = MatDat(SliceSampler1D(LogDens.Almagro.Z_cond_U, L, U, x0),1,1);

WriteLn("TRACE x="<<x);