Conditional distribution of a discrete random variable given a continuous one I can´t solve a conditional distribution problem using the simple Bayes´ Theorem for PDFs.
The thing is, I know that X|Y ~ Poisson(Y) and that Y ~ Gamma(α,β). I´m being asked to find the PDF for Y|X but I´m complete clueless about how to do so. Can anyone point me in the right direction?
 A: (hint hint)
Try starting out by writing out the Bayes formula with the full PDFs of $p(X|Y)$ and $p(Y)$ (we can ignore $p(X)$: can you tell me why?). Multiple $p(X|Y)$ and $p(Y)$ together and see what terms you can combine. It should start to look familiar.
Here, I'll get you started. Let's start out simple (NB: this isn't the full likelihood below, it's really just the likelihood of a single observation of $x$. What does the likelihood for multiple observations of $x$ look like?):
$
\begin{align}
BayesNumerator &= p(X|Y)  p(Y) \\
&= \left ( \frac{y^x }{x!}e^{-y} \right )
   \left ( \frac{\beta^\alpha}{\Gamma\left( \alpha \right)}y^{\alpha-1} e^{-\beta y}     \right ) \\
&= \left (? \right ) y^{(?)}e^{(?)}
\end{align}$ 
Can you fill in the blanks?
UPDATE:
Importing from the google doc you linked in the comments (for posterity), you've been able to reduce the numerator of the posterior to 
$$p(X|Y)p(Y) = e^{-y(\beta+1)} y^{x+\alpha-1} \frac{\beta^\alpha}{\Gamma\left(\alpha\right)}\frac{1}{x!}$$
And you've correctly determined that the posterior is only a function of $y$, so the $\alpha$, $\beta$, and $x$ terms are actually constants in the posterior formulation. 
Some additional thoughts:


*

*What terms in the posterior are constant with respect to $y$?

*Let's take a look back at the denominator. Earlier I prompted you that we could ignore the denominator. There's a trick there, but since you don't know it, I'm going to assume that deriving it is part of your assignment. Instead of taking you all the way, here's a motivating example:


I suggested to you earlier that you had a conjugate prior, and you've correctly guessed that your posterior will be gamma but you're still confused about the parameters that it will take. Instead of a gamma-poisson model, let's consider a beta-bernoulli model.
$p(W|\theta) \sim Bernoulli(\theta)$
$p(\theta) \sim Beta(\alpha, \beta)$
Again, assuming a single observation from $W$ and constraining our attention to the numerator of $p(\theta|W) = \frac{p(W|\theta)p(\theta)}{p(W)}$ :
\begin{equation}
\begin{split}
p(W|\theta)p(\theta) &= \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}*\theta^w(1-\theta)^{1-w} \\
&=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} \\
\end{split}
\end{equation}
Now, this looks sort of Beta-ish, but we're not quite there yet. The parameters in the ratio on the left seem to disagree with the parameters in the multiplication on the right. This is because we've ignored a "normalizing constant," the denominator in the bayesian formula, the marginal distribution of $p(W)$. We often call this a "normalizing constant" because it's constant with respect to the variable of interest (in this case $\theta$, in your problem $Y$) and it's "normalizing" in the sense that it ensures that the pdf integrates to 1. 
Remember that. By definition, a continuous PDF integrates to 1, and a discrete PDF sums to 1 (basically the same thing). 
With that in mind, let's revisit the denominator. 
We wrote the denominator as $p(W)$, but an equivalent formulation is $p(W) = \int_\theta p(W|\theta)p(\theta)d\theta$. This is also known as "marginalizing out" the parameter theta. We're essentially evaluating the likelihood for all possible values of theta and summing the results. Now, we know $p(W)$ is a valid pdf, therefore, we know that this term integrates to 1. This is tremendously useful information. Let me show you why. 
From earlier, I showed that
$p(W|\theta)p(\theta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1}$
$
\begin{align}
\therefore \int_\theta p(W|\theta)p(\theta) d\theta &= \int_\theta \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} d\theta \\
&=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \int_\theta  \theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} d\theta
\end{align}
$
I pulled that ratio out because it's not in terms of $\theta$, the parameter we're integrating with respect to (and not coincidentally the parameter that the posterior is in terms of). 
Allright, time for the mind blowing part: if you squint very carefully, that integrand looks very similar to the Beta pdf. If we refer to the portion of a pdf that is in terms of the "input" variable as the "kernel" of the distribution, than the integrand is the kernel of $Beta(\alpha+w, \beta+1-w)$. If we can get the appropriate normalizing constant into the integrant, then the integrand completely disappears because it integrates to 1. Bam. Magic. Let's do that. 
The normalizing constant we need for a $Beta(\alpha+w, \beta+1-w)$ distribution is:
$$\frac{\Gamma(\alpha+w)\Gamma(\beta+1-w)}{\Gamma(\alpha+w+\beta+1-w)} = \frac{\Gamma(\alpha+w)\Gamma(\beta+1-w)}{\Gamma(\alpha+\beta+1)}$$ 
We can multiply anything we want by "1", so if we multiply the integrand above with the ratio of the normalizing constant with itself, we can push it into the integration:
$\int_\theta p(W|\theta)p(\theta) d\theta = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \int_\theta  \theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} d\theta$
$= \left[\frac{\Gamma(\alpha+w)\Gamma(\beta+1-w)}{\Gamma(\alpha+\beta+1)} 
\frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+w)\Gamma(\beta+1-w)}
\right] \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \int_\theta  \theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} d\theta$
$=\left[
\frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+w)\Gamma(\beta+1-w)}
\right] \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \int_\theta  \frac{\Gamma(\alpha+w)\Gamma(\beta+1-w)}{\Gamma(\alpha+\beta+1)} \theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} d\theta$
$=\frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+w)\Gamma(\beta+1-w)}
\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} (1)$
$=\frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+w)\Gamma(\beta+1-w)}
\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}$
Ta da! The denominator is constant with respect to $\theta$, and we just evaluated it exactly. So let's put the numerator and denominator together and see what happens:
$p(W|\theta)p(\theta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1}$
$p(W) = \frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+w)\Gamma(\beta+1-w)}
\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}$
$
\begin{align}
\therefore p(\theta|W) &= \frac{p(W|\theta)p(\theta)}{p(W)} \\
&= \left [ \frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+w)\Gamma(\beta+1-w)}
\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \right ]^{-1} \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1} \\
&= \frac{\Gamma(\alpha+w)\Gamma(\beta+1-w)}{\Gamma(\alpha+\beta+1)} \theta^{\alpha+w-1}(1-\theta)^{\beta+1-w-1}
\end{align}
$
Which is exactly the $Beta(\alpha+w, \beta+1-w)$ pdf. Because the numerator is the integrand of the denominator, the denominator gives us exactly the term we need to cancel out the dangling constant in the numerator that is in terms of the prior parameters and replace it with the constant in terms of the posterior parameters we need to get the pdf we recognized in the kernel. Neat trick, right? Welcome to the wonderful world of conjugate priors. 
This looks complicated, but it's really not that bad. You're almost there. You're so close, I can taste it, and I hope you can taste it too (maybe I'm a little drunk, but it's making me productive). You just need to apply the same procedure I applied here to your problem. You've already reduced the numerator: throw that term you've got into the denominator and try and apply an analogous trick I showed you to get it to integrate to 1. You said it looks similar to a gamma distribution right? That's because you have a gamma kernel. Now all you need is the proper normalizing constant. 
PS: I've been doing all of this assuming we only have a single observation. If your sampling distribution is distributed according to $f(x)$, then the full likelihood is $\prod _{i=1}^n f(x_i)$. In your case, this will give you a sufficient statistic in your likelihood that is a summation. Once you solve this for a single observation, you shouldn't have any trouble generalizing to this slightly more complex likelihood.
Ready set go!
