Me gustaría generar muestras aleatorias de una distribución normal bivariante bajo una condición. La primera variable normal es $\varepsilon_1$ , y en segundo lugar normal es $\varepsilon_2$. La condición es $\varepsilon_1>T_1$ donde $T_1$ es una constante, y $ a \varepsilon_1 + b\varepsilon_2 <T_2$ donde $a$, $b$, y $T_1$ son constantes. $\varepsilon_1$ $\varepsilon_2$ son independientes. Por tanto, las condiciones son: la generación de una región en el espacio 2D delimitada por la línea vertical $T_1$ y una línea inclinada. Es allí una manera de hacer esto sin generar muchas muestras aleatorias y tirar los que están fuera de la condición de área? La razón es la probabilidad de que en la región de la condición puede ser muy pequeño, así que a tirar de las muestras no es una opción.
Respuestas
¿Demasiados anuncios?Si usted tuvo otro obligado (como $\epsilon_2 > T3$), se muestra de manera uniforme y, a continuación, los pesos de la muestra utilizando la densidad normal bivariada. Usted tendría cero rechazo. Tal vez en tu aplicación no es demasiado irrazonable imponer un límite?
Probablemente mejor:
Usted encontrar la intersección entre las dos condiciones lineales. Después de generar una r.v. $x_1$ a partir de una exponencial o trunca normal a lo largo de una de las dos condiciones (dicen que a lo largo de $\epsilon_1 = T_1$). Entonces, si el ángulo entre el 2 condiciones lineales es aguda, se dibuja de manera uniforme (una perpendicular a $\epsilon_1 = T_1$) a lo largo de la línea entre el$x_1$$a\epsilon_1 + b\epsilon_2 = T_2$. Si es obtuso, dibujar perpendicularmente a $\epsilon_1 = T_1$ a partir de una normal truncada o exponencial. No hay ningún rechazo involucrados, y no es necesario que el área a ser delimitada, pero se obtiene una muestra ponderada.
He utilizado el muestreo de Gibbs enfoque. De esta manera solo el principio del muestreo de Gibbs es expulsado (periodo de estabilización). Por lo tanto el número de muestras de talle es que no aumenta con el número de muestras.
- Condicional en la observación de $\varepsilon_1$, $\varepsilon_2$ es la muestra de una distribución normal con enlazados $b\varepsilon_2< Th_2 - a\varepsilon_1$.
- Condicional en la observación de $\varepsilon_2$, $Th_1<\varepsilon_1< (Th_2 - b\varepsilon_2)/a$.
A continuación el código establece $a=\sqrt{t1}$, $b=\sqrt{t2-t1}$.
nScens = 1E8;
epsilon1 = randn(nScens, 1);
epsilon2 = randn(nScens, 1);
Th1 = -3;
Th2 = -2.9;
t1 = 700;
t2 = 707;
ind = epsilon1 > Th1 & ( epsilon1*sqrt(t1) + epsilon2*sqrt(t2-t1))/sqrt(t2) < Th2;
sum(ind)
figure(1)
subplot(121)
scatter(epsilon1(ind), epsilon2(ind),'.' )
axis([ -3 -2.5 -5 1])
subplot(122)
smoothhist2D([epsilon1(ind), epsilon2(ind)],5, [100,100],[], 'contour')
axis([ -3 -2.5 -5 1])
% gibbs sampler
nGibbs = 75000;
epsilon1Gibbs = 0;
for i=1:nGibbs
epsilon2Gibbs = norminv( normcdf( (Th2*sqrt(t2) - epsilon1Gibbs*sqrt(t1) )/sqrt(t2-t1) )*rand );
p = ( -normcdf(Th1) + normcdf( (Th2*sqrt(t2) - epsilon2Gibbs*sqrt(t2-t1) )/sqrt(t1) ) )*rand + normcdf(Th1);
epsilon1Gibbs = norminv( p );
epsilonGibbs(i, :) = [epsilon1Gibbs epsilon2Gibbs];
end
indGibbs = 2500:nGibbs;
figure(2)
subplot(121)
scatter(epsilonGibbs(indGibbs,1), epsilonGibbs(indGibbs,2),'.' )
axis([ -3 -2.5 -5 1])
subplot(122)
smoothhist2D( epsilonGibbs(indGibbs,:) ,5, [100,100],[], 'contour')
axis([ -3 -2.5 -5 1])
La fuerza bruta de muestreo:
Muestreo de Gibbs:
Un enfoque sencillo que implicaría una enorme reducción en la tasa de rechazo sería para girar las coordenadas $(\epsilon_1,\epsilon_2)$ decir $(X_1,X_2)$ de manera tal que la línea de $aε_1+bε_2=T_2$ se convierte en vertical ($cX_1=\tau_2$, por ejemplo). Luego se genera a partir de la truncada normal tal que $cX_1<\tau_2$. A continuación, generar un independiente $X_2$ y rechazar aquellos pares que no la otra (girado) condición, y girar la aceptación de los pares de la espalda.
La tasa de rechazo será sustancial (es probable que supere el 50%, por ejemplo), pero probablemente no va a estar en todos los extremos, como sin duda lo sería si no generar desde el extremo de la cola truncada normal, para empezar.