id	summary	reporter	owner	description	type	status	priority	milestone	component	version	severity	resolution	keywords	cc
62	Problema de la funci� GibbsSampler	jcaballero	Jorge	"Hola a todos. Programando una distribución multivariante con GibsSampler hemos 
apreciado una posible carencia en la función. A saber:
 
  - No permite introducir condiciones iniciales.
  - No permite introducir un contador para generar la función

Nos explicamos mejor con el siguiente ejemplo:

///////////////////////////////////////////////////////////////////////
// NormalMultivariate.TOL
//
// PORPUSE: Generacion por muestreo de una Multinormal utilizando el método
//          de Gibs Sampler
///////////////////////////////////////////////////////////////////////
Set Include(""../sample/fun/funstat.tol"");
///////////////////////////////////////////////////////////////////////
// 1º Introduccion de datos
///////////////////////////////////////////////////////////////////////
Struct multinormal_data 
{
  Matrix mean,
  Matrix vars
};

Real dim  = 2;    // dimensión de la normal
///////////////////////////////////////////////////////////////////////
// 2º Creación del vector de medias y matriz varianzas covarianzas
///////////////////////////////////////////////////////////////////////

// Vector de medias
Set nu_ = For(1,dim,Set(Real k)
{
  [[Round(Rand(0,20))]]
});
Matrix nu = SetMat(nu_);
Set var_ = For(1,dim,Set(Real k)
{
  Set filasup = For(1,k,Real(Real j)
  {
     Round(Rand(1,5))
  });
  Set filainf = For(k+1,dim,Real(Real m)
  {
    0
  });
  filasup << filainf
});
Matrix var = SetMat(Traspose(var_))*SetMat(var_);

Set medias = For(2,dim,Set(Real i)
{
  Matrix mu2 = Sub(nu,i,1,1,1);
  Matrix mu1 = Sub(nu,1,1,i-1,1);
  SetOfMatrix(mu1,mu2)
});

Set mydata = multinormal_data(nu,var);
Set indice = SetOfReal(1);
///////////////////////////////////////////////////////////////////////
// 3º  Fórmulas condicionadas.
///////////////////////////////////////////////////////////////////////

Set NormalConds = BinGroup(""<<"",For(1,dim,Set(Real k)
{
    Real NormalCond(Real X, Matrix par, Real prev, Set data)
    {
      Matrix media  = data->mean;
      Matrix varia  = data->vars;
      
/// AQUÍ ESTÁ EL PROBLEMA ¿CÓMO SE PUEDE HACER QUE indice SE MUEVA
    CON EL BUCLE?
      Set momentos = MomtCondMultNormal(par,indice,media,varia);
      
      Real meanY = MatSet(momentos[1])[1][1];
      Real varY  = MatSet(momentos[2])[1][1];
     -0.5 * Pow((X-meanY)/varY,2)
    };
    [[NormalCond]]
}));

Set FullCondSpec1 = For(1,dim,Set(Real j)
{
  [[NormalConds[j],mydata,-100,100,1]]
});

WriteLn (""-------------------"");
WriteLn ("" Metodo1 :Sample1  "" );
WriteLn (""-------------------"");
Matrix Sample1 = GibbsSampler(FullCondSpec1, 100, 100);

-----------------------------

Obsérvese que al ir generando las diversas normales multivariantes no es factible 
declarárselas una a una, ya que se busca generalidad, para cualquier n. Por eso la 
idea de crear una variable índice que vaya generando las distintas condicionadas.
Pero desafortunadamente la función:

  Real NormalCond(Real X, Matrix par, Real prev, Set data)

del código del ejemplo, no permite introducirle un 5º parámetro que sea el índice 
para controlar las normales(definido en el código como global con valor 1) ya que 
GibsSampler es estricto en el sentido que sólo admite estos 4 parámetros.

Tampoco sirve mover el índice en el bucle pues la función no lo ""recuerda"" por ser 
variable local.

¿Qué opinais? Una posible solución que se me ocurre es permitir que GibsSampler 
admita la posibilidad de un quinto valor real.

Un saludo a todos."	enhancement	closed	normal		Math	head	minor	later