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

#849 closed defect

Muestrear una normal con restriciones de igualdad — at Version 1

Reported by: jlaybar Owned by: Víctor de Buen Remiro
Priority: normal Milestone: Mantainance
Component: Math Version: 2.0.1
Severity: blocker Keywords:
Cc:

Description (last modified by Víctor de Buen Remiro)

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);

Change History (1)

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

Component: BSRMath
Description: modified (diff)
Milestone: Manteinance
Status: newaccepted
Version: 2.0.1

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.

Note: See TracTickets for help on using tickets.