I would like to generate random samples from a bivariate normal distribution under a condition. First normal variable is $\varepsilon_1$ , and second normal is $\varepsilon_2$. The condition is $\varepsilon_1>T_1$ where $T_1$ is a constant, and $ a \varepsilon_1 + b\varepsilon_2 <T_2$ where $a$, $b$, and $T_1$ are constants. $\varepsilon_1$ and $\varepsilon_2$ are independent. Thus the conditions are generating a region in 2D space bounded by the vertical line $T_1$ and a tilted line. Is there a way to do this without generating many random samples and throwing the ones outside the condition area? The reason is the probability in the region of condition can be quite small, thus throwing away samples is not an option.
$\begingroup$
$\endgroup$
8
-
1$\begingroup$ When you say throwing away values is not an option, do you mean that there's some reason you can't throw any away, or that you can't throw as many away as a naive approach would suggest? $\endgroup$– Glen_bCommented Jan 2, 2014 at 22:53
-
$\begingroup$ Depending on $T_1$ and $T_2$ acceptable region can be far away from the center (0,0) of 2D space of $\varepsilon_1$ and $\varepsilon_2$ standard normals ($N(0,1)$). Thus naive approach of just sampling then throwing away might require in some cases more than a billion samples to get acceptable 1000 pairs of $\varepsilon_1$ and $\varepsilon_2$. A throw away method will work if we are still sampling around the acceptable region, but we need to throw away a small amount of samples for some reason, and use the rest for Monte Carlo simulation. What do you have in mind, Gibbs sampling?? $\endgroup$– adamCommented Jan 3, 2014 at 8:48
-
$\begingroup$ Not Gibbs, no; much simpler rejection algorithms should work. Can you give a typical T1, T2 a and b? $\endgroup$– Glen_bCommented Jan 3, 2014 at 11:42
-
$\begingroup$ For example $T_1=-3$, $T_2=-2.9$, $a=0.995$, $b=0.0995$ $\endgroup$– adamCommented Jan 3, 2014 at 13:41
-
$\begingroup$ And the normals are standard ($\sigma=1$)? $\endgroup$– Glen_bCommented Jan 3, 2014 at 15:24
|
Show 3 more comments
3 Answers
$\begingroup$
$\endgroup$
2
If you had another bound (such as $\epsilon_2 > T3$), you could sample uniformly and then weights the sample using the bivariate normal density. You would have zero rejection. Maybe in your application it is not too unreasonable to impose such a bound?
Probably better:
You find the intersection between the two linear conditions. Then you generate a r.v. $x_1$ from an exponential or a truncated normal along one of the two conditions (say along $\epsilon_1 = T_1$). Then, if the angle between the 2 linear conditions is acute, you draw uniformly (an perpendicularly to $\epsilon_1 = T_1$) along the line between $x_1$ and $a\epsilon_1 + b\epsilon_2 = T_2$. If it is obtuse, you draw perpendicularly to $\epsilon_1 = T_1$ from a truncated normal or exponential. There is no rejection involved, and you don't need the area to be bounded, but you get a weighted sample.
answered Jan 3, 2014 at 17:10
-
$\begingroup$ Let me make sure if I understand it correctly. In your suggestion of uniform sampling along the vertical line to $T1$ for each sample we are going to find the weight by multiplying normal probability densities for sample $\varepsilon_1$ and $\varepsilon_2$. Thus each sample will have a different weight. I think it is a good idea, but I would prefer to have equally weighted samples, since I look at the envelope of the distribution. It becomes a little complicated to back out the envelope from the weighted samples. $\endgroup$– adamCommented Jan 6, 2014 at 12:21
-
$\begingroup$ Yes, that was my idea. Generally when I have a weighted sample I just re-sample with probability proportional to the weights. In that way you obtain an equally weighted sample out of a weighted one. $\endgroup$ Commented Jan 6, 2014 at 14:54
$\begingroup$
$\endgroup$
1
I have used the Gibbs sampling approach. This way only the beginning of the Gibbs sampling is thrown out (stabilization period). Thus number of waisted samples is not increasing with the number of required samples.
- Conditional on observing $\varepsilon_1$, $\varepsilon_2$ is sampling from normal distribution with bound $b\varepsilon_2< Th_2 - a\varepsilon_1$.
- Conditional on observing $\varepsilon_2$, $Th_1<\varepsilon_1< (Th_2 - b\varepsilon_2)/a$.
Below code sets $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])
Brute force sampling:
Gibbs sampling:
-
$\begingroup$ Nice, how much did you gain by using Gibbs in terms of rejection rate? $\endgroup$ Commented Jan 6, 2014 at 14:49
$\begingroup$
$\endgroup$
One simple approach that would involve a huge reduction in the rejection rate would be to rotate the coordinates $(\epsilon_1,\epsilon_2)$ to say $(X_1,X_2)$ such that the line $aε_1+bε_2=T_2$ becomes vertical ($cX_1=\tau_2$, say). Then generate from the truncated normal such that $cX_1<\tau_2$. Then generate an independent $X_2$ and reject those pairs which fail the other (rotated) condition, and rotate the accepted pairs back.
The rejection rate will be substantial (it will likely exceed 50%, for example), but probably won't be at all extreme, as it certainly would be if you didn't generate from the extreme-tail truncated normal to begin with.
