Context Navigation

← Previous Ticket
Next Ticket →

close Warning: Can't synchronize with repository "(default)" (/var/svn/tolp does not appear to be a Subversion repository.). Look in the Trac log for more information.

Opened 15 years ago

Last modified 14 years ago

#888 closed defect

La función SliceSampler1D falla con dominios inconexos — at Version 1

Reported by:	Víctor de Buen Remiro	Owned by:	Víctor de Buen Remiro
Priority:	highest	Milestone:	Numerical methods
Component:	BSR	Version:	2.0.1
Severity:	blocker	Keywords:
Cc:

Description (last modified by 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(440791623);
//////////////////////////////////////////////////////////////////////////////
Real PutRandomSeed(440791623);
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);

Estas son las salidas que da el programa anterior para cada una de las semillas elegidas:

// 3 de 8 no estacionarias incluida la última: devuelve el valor inicial
TRACE  rndSeed:1836343167
TRACE  SetValue(-0.09944049062890037) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09944049062890037) :-1387.086986834382
TRACE  SetValue(-0.9900840071843289) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.9900840071843289) :-23313.37238721964
TRACE  SetValue(0.009915992815671126) isStationary:0
TRACE  LogDens.Almagro.Z_cond_U(0.009915992815671126) :-1.#INF
TRACE  SetValue(-0.9900840071843289) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.9900840071843289) :-23313.37238721964
TRACE  SetValue(0.009915992815671126) isStationary:0
TRACE  LogDens.Almagro.Z_cond_U(0.009915992815671126) :-1.#INF
TRACE  SetValue(-0.8559963063218914) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.8559963063218914) :-17243.20337303151
TRACE  SetValue(-0.6656894964087795) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.6656894964087795) :-10313.15043543136
TRACE  SetValue(-0.02358345558353336) isStationary:0
TRACE  LogDens.Almagro.Z_cond_U(-0.02358345558353336) :-1.#INF
TRACE x=-0.09944049062890037

//ninguna no estacionaria: devuelve una simulación válida
TRACE  rndSeed:1664310743
TRACE  SetValue(-0.09944049062890037) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09944049062890037) :-1400.100045170432
TRACE  SetValue(-0.09271461170849427) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09271461170849427) :-1400.667034074515
TRACE  SetValue(-0.09271461170849427) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09271461170849427) :-1400.667034074515
TRACE  SetValue(-0.2891919542795427) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.2891919542795427) :-2378.702930021611
TRACE  SetValue(-0.2776974464241599) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.2776974464241599) :-2264.760232978221
TRACE  SetValue(-0.1911206362276727) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1911206362276727) :-1633.055504809866
TRACE  SetValue(-0.1249158202167721) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1249158202167721) :-1419.837272288265
TRACE  SetValue(-0.1183620014027887) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1183620014027887) :-1411.451378464228
TRACE  SetValue(-0.1146000550613002) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1146000550613002) :-1407.673101525984
TRACE  SetValue(-0.1112400614714343) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1112400614714343) :-1404.936913722732
TRACE  SetValue(-0.1101893541977121) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1101893541977121) :-1404.204911943004
TRACE  SetValue(-0.1040738663446942) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.1040738663446942) :-1401.113552350067
TRACE  SetValue(-0.09290987134805526) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09290987134805526) :-1400.616555105531
TRACE  SetValue(-0.09519783142933348) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09519783142933348) :-1400.176636667713
TRACE  SetValue(-0.09717720524250692) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09717720524250692) :-1400.021389962531
TRACE x=-0.09717720524250692

//2 de 8 no estacionarias, pero no la última: devuelve una simulación válida
TRACE  rndSeed:440791623
TRACE  SetValue(-0.09944049062890037) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09944049062890037) :-1414.620951115607
TRACE  SetValue(-0.8534963953190766) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.8534963953190766) :-15809.41847702517
TRACE  SetValue(0.1465036046809234) isStationary:0
TRACE  LogDens.Almagro.Z_cond_U(0.1465036046809234) :-1.#INF
TRACE  SetValue(-0.8534963953190766) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.8534963953190766) :-15809.41847702517
TRACE  SetValue(0.1465036046809234) isStationary:0
TRACE  LogDens.Almagro.Z_cond_U(0.1465036046809234) :-1.#INF
TRACE  SetValue(-0.407688048935803) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.407688048935803) :-3840.162568990362
TRACE  SetValue(-0.2629580672722903) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.2629580672722903) :-2105.638034846308
TRACE  SetValue(-0.09840860067143334) isStationary:1
TRACE  LogDens.Almagro.Z_cond_U(-0.09840860067143334) :-1414.534033071261
TRACE x=-0.09840860067143334

Change History (1)

comment:1 Changed 15 years ago by Víctor de Buen Remiro

Description:	modified (diff)

Note: See TracTickets for help on using tickets.

Download in other formats: