conditional sampling of bivariate normals 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.
 A: 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. 
A: 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:

A: 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.
