Priors on nonlinear parameters in BSR

[Pedro Gea]

Estoy buscando introducir información a priori en los parámetros no lineales.

Soy capaz de ponerles una restricción de dominio, pero no sé (ni siquiera si se puede o qué sentido tiene) ponerle una distribución a priori.

Intenté ponerle un prior normal (como si fuese un parámetro lineal o el de un omitido) definiéndo un bloque "primary" pero obtengo el error:

Nombre de índice duplicado Veh.XXX.Mat__Pib.es.Men__SDelta.7

Miré en el código de BysMcmc y encontré un método get.priorLogDens pero tampoco sé de qué se trata.

[Víctor de Buen Remiro]

La información a priori tiene el mismo sentido en este caso que sobre cualquier
otro parámetro, aunque la implementación, que no siempre está disponible (como el
caso ARMA), es muy diferente en cada bloque de Gibbs. Por eso no tiene sentido
intentar aplicársela como si fuera un parámetro lineal, pues la fórmula de actuación
es completamente distinta.

Afortunadamente sí está programada para los filtros no lineales, y de hecho ibas
muy bien encaminado: sólo hay que crear en el NameBlock que hace de filtro un método

Real get.priorLogDens(Real paramIdx)

que devuelva el logaritmo de la densidad del parámetro número paramIdx.

Hay un ejemplo de prior normal en
_DeltaPrior.tol

  ///////////////////////////////////////////////////////////////////////////
  Real get.priorLogDens(Real paramIdx)
  ///////////////////////////////////////////////////////////////////////////
  {
  //En este caso no se usa el parámetro paramIdx porque sólo hay prior sobre
  //el parámetro _.delta, concretamente un prior normal.
    Log(DensNormal(_.delta, _.delta.info->Nu, _.delta.info->Sigma))
  //Es lo mismo que esto otro, salvo una cosntante irrelevante, pero en TOL 
  //es más rápido la anterior fórmula
  //Log(-0.5*((_.delta-_.delta.info->Nu)(_.delta.info->Sigma))^2)
  };

En el capítulo 6.2.1 del manual de BSR (pag. 48), está la descripción matemática
para el caso normal, auqnue en realidad es extensible a cualquier distribución
escalar.

Podría pensarse incluso en un sistema que ofreciera la distribución conjugada
en los ​casos más típicos

[Pedro Gea]

Las indicaciones anteriores son suficientes para solucionar las dudas planteadas.

Quizá lo único que comentaría, porque quizá no es muy común, es que el nameblock que representa al filtro ha de tener un atributo dinámico (atributo cuyo valor va cambiando) que represente a la matriz de parámetros, ya que esta matriz no se pasa como argumento en la llamada a get.priorLogDens y sin embargo es necesaria para su evaluación. Este atributo (_.delta en el ejemplo) ha de actualizarse en cada llamada al método eval que sí lo recibe como argumento.

[Víctor de Buen Remiro]

El lugar adecuado para esa actualización es

Real set.parameter(Matrix paramValues)

[Pedro Gea]

Sí, yo he creado un método igual que ése, pero ese método no es obligatorio para el filtro, ni es llamado por el objeto @BlockNonLin en ningún lado. Se puede crear o no, en cualquier caso lo que sí hay que hacer es actualizar el atributo dinámico en la llamada al método eval, sea mediante una llamada a ese método set.parameter o no.

[Pedro Gea]

Se encuentra una advertencia con un mensaje incorrecto:

En RealGSLFunction::Eval, valor desconocido en la evaluación de LogDens.ARIMA

cuando el método get.priorLogDens devuelve un omitido. Debería usarse un mensaje general para cualquier tipo de parámetro no lineal (no ha de ser ARIMA). Si no fuera muy complicado quizá podría indicarse el nombre del segmento o algo así.

La advertencia surge cuando se pone un prior muy ajustado:

Real Log(DensNormal(0.5, 0.92, 0.01)) // -> ?

ya que para la función Log el logaritmo de 0 es desconocido y no menos infinito. Aunque si no me equivoco si se devuelve menos infinito (-1/0) la simulación del parámetro se queda atrancada en el valor inicial, así que entiendo que debería devolver un valor negativo grande como -10^10 o algo así.

[Víctor de Buen Remiro]

¿Sería posible hacerme llegar una forma de reproducir esto?
Con las explicaciones que me das sólo me hago una idea perro no soy capaz de localizarlo exactamente para arreglarlo, y menos aún de comprobarlo posteriormente.

Mientras tanto podrías usar una forma de cálculo más adecuada a la aritmética discreta. En general, siempre que se quiere el logaritmo de una expresión que es producto o cociente de varias siempre es mejor hacer el logaritmo de cada factor y luego hacer la suma o resta correspondiente.

Por ejemplo en este caso, si definimos la función

Real LogDensNormal(Real x, Real nu, Real sigma)
{
  -0.5*Log(2*pi)-Log(sigma)-0.5*((x-nu)/sigma)^2
};

entonces la expresión

Real LogDensNormal(0.5, 0.92, 0.01);

da -878.3137683472166, un número bastante más modesto que -10^10

Quizás sería bueno que se dispusiera de las correspondientes funciones LogDens*. Por el momento esta de la normal será incluida en la siguiente release.

[Víctor de Buen Remiro]

(In [3452]) Refs #1150