Problema de la funci� GibbsSampler

[jcaballero]

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.

[Jorge]

this is an enhancement and will be implemented.

[Jorge]

marked as LATER because the new function MetropolisHastings solve this problem.
May be GibbsSampler could be declared DEPRECATED.