Distribution of one Poisson variable relative to another Assume A and B are two independent Poisson random numbers with means a and b.
What is the conditional distribution of A given that A>B, A=B and A<B respectively, i.e. out of all the samples, what is the distribution of A when its larger, equal or smaller than B?
Note that I am not asking for the distribution of the largest or smallest of the two. I am specifically interested in the first variable. 
I can use Skellam to get probabilities for the occurrence of each of these 3 scenarios, but this does not help me getting the distribution of A itself. I fear they can not be Poisson since Prob (A=0)=0 when A>B.
 A: Let's tackle the equality problem first, then the other two:
$$P(A=i|A=B) = \frac{P(A=i)P(B=i)}{\sum_j P(A=j)P(B=j)}$$
$$= \frac{\exp(-a)a^i/i! \cdot \exp(-b)b^i/i! }
  {\sum_j \exp(-a)a^j/j! \cdot \exp(-b)b^j/j!}$$
$$= \frac{ (ab)^i /(i!)^2 }{ \sum_j (ab)^j/(j!)^2}\,,\: i=0,1,2,...$$
Which can be further simplified but I think the principle is clear. Even this form is quite suitable for numerical calculation (the values decrease quite rapidly to the right of the mode, so good accuracy can usually be had simply by truncating the series at a sensible place. If the parameters are very large you might need to be a bit more sophisticated about it. 
I simulated 10000 values each from $A$ and $B$ with $a=2$ and $b=1$ to compare with that calculation. It turned out that 21.4% of the simulations gave A=B, and the proportion of A at each value when that was the case was as shown in the histogram:

The blue points represent the corresponding theoretical calculation above. As you can see the values are very close together.

$$P(A=i|A<B) = \frac{P(A=i)\sum_{k=i+1}^\infty P(B=k)}{\sum_j (P(A=j)\sum_{k=j+1}^\infty P(B=k))}$$
$$ = \frac{P(A=i)S_B(i)}{\sum_j P(A=j)S_B(j)}$$
Where $S_B$ is the survivor function for variable $B$. It would be possible to cancel factors of $e^{-(a+b)}$ as in the previous case but it's convenient to leave them there and do the calculations directly in this form.
Here's that calculation done in R code for that cut off somewhere reasonable for the parameters involved (2 and 1 again):
i <- 0:20
p <- dpois(i,2)*ppois(i,1,lower.tail=FALSE)
p <- p/sum(p)

And here's a comparison with simulation:


Now the $A>B$ case:
$$P(A=i|A>B) = \frac{P(A=i)\sum_{k=0}^{i-1} P(B=k)}{\sum_j (P(A=j)\sum_{k=0}^{j-1} P(B=k))}$$
$$ = \frac{P(A=i)F_B(i-1)}{\sum_j P(A=j)F_B(j-1)}\,,\: i=1,2,...$$
And again we can simulate and compare simulated proportions with the above algebraic calculations:

