Muestrear una normal con restriciones de igualdad

[jlaybar]

Hola Victor estoy intentando muestrear dos variables aleatorias cuya suma sea igual a uno x1+x2=1 y x1>0, x2>0

Si embargo con la función RandConstrainedMNormal no lo consigo

Matrix mu = Col(0.2,0.8);
Matrix Cov =((1,0),(0,1));
Matrix B = (
 (1,1),
(-1,-1),
(-1,0),
(0,-1)
);
Matrix b = Col(1,1,0,0);
Matrix sample = RandConstrainedMNormal(mu, Cov, B, b);

[Víctor de Buen Remiro]

Aunque no he tenido tiempo imagino que el problema es que la distribución coin restricciones de igualdad está degenerada por definición y el método empleado no está preparado para ello.
No sé si se podrá modificar para que lo admita pero no lo creo.
En estos casos lo mejor es aplciar ingeniería inversa, es decir, reducir las dimensiones al caso del rango completo y generar en la nueva base para luego, mediante las ecuaciones de reducción invertir el cambio de base y obtener las simulaciones en el espacio original.

En este caso el cambio de variable es

Y = X1+X2 ~ N(nu_y = 0.2+0.8=1, var_y=1+1=2);
Y>=0

1º) Se genera una realización y de Y con una Normal(1,2) truncada en [0,+Inf)

2º) Se genera una realización x1 de X1 condicionado a Y=y con una Normal(0.2,1) truncada en [0,y]

3º) Se calcula la realización x2 de X2 que es determinista condicionada a X1=x1 e Y=y x2 = y - x1

Espero que te sirva. Para un caso general puede ser bastante complicado invertir las condiones de desigualdad pero iré pensando en cómo resolverlo.

[Víctor de Buen Remiro]

Este ticket está generalizado en el #861