1. The Likelihood Function
We observe $n$ independent data points $y_1, y_2, \dots, y_n$ where each $y_i$ follows a Poisson distribution with rate parameter $\lambda > 0$:
$$y_i \mid \lambda \sim \text{Poisson}(\lambda)$$
The probability mass function for a single observation is:
$$p(y_i \mid \lambda) = \frac{\lambda^{y_i} e^{-\lambda}}{y_i!}$$
Since the $n$ observations are independent, the joint likelihood is the product of individual likelihoods:
$$\mathscr{L}(\lambda \mid \mathbf{y}) = \prod_{i=1}^n \frac{\lambda^{y_i} e^{-\lambda}}{y_i!}$$
$$= \frac{\lambda^{\sum_{i=1}^n y_i}\, e^{-n\lambda}}{\prod_{i=1}^n y_i!}$$
Denoting $\displaystyle\sum_{i=1}^n y_i = S$ for convenience:
$$\boxed{\mathscr{L}(\lambda \mid \mathbf{y}) = \frac{\lambda^{S}\, e^{-n\lambda}}{\prod_{i=1}^n y_i!}}$$
2. The Prior Distribution
Since $\lambda > 0$, a natural and conjugate prior is the Gamma distribution:
$$\lambda \sim \text{Gamma}(\alpha, \beta)$$
The Gamma distribution has the probability density function:
$$p(\lambda) = \frac{\beta^\alpha}{\Gamma(\alpha)}\, \lambda^{\alpha-1} e^{-\beta\lambda}, \quad \lambda > 0$$
where $\alpha > 0$ is the shape hyperparameter and $\beta > 0$ is the rate hyperparameter, and the Gamma function is:
$$\Gamma(\alpha) = \int_0^\infty t^{\alpha-1} e^{-t}\, dt$$
The prior mean is $\dfrac{\alpha}{\beta}$ and the prior strength is interpreted through $\alpha$ and $\beta$ together — $\beta$ can be thought of as the number of prior observations and $\alpha$ as the prior total count.
3. Bayes’ Theorem from First Principles
From the definition of conditional probability, the joint distribution factors two ways:
$$p(\lambda,\, \mathbf{y}) = p(\mathbf{y} \mid \lambda)\cdot p(\lambda)$$
$$p(\lambda,\, \mathbf{y}) = p(\lambda \mid \mathbf{y})\cdot p(\mathbf{y})$$
Setting these equal and dividing both sides by $p(\mathbf{y})$:
$$\boxed{p(\lambda \mid \mathbf{y}) = \frac{p(\mathbf{y} \mid \lambda)\cdot p(\lambda)}{p(\mathbf{y})}}$$
where the marginal likelihood is:
$$p(\mathbf{y}) = \int_0^\infty p(\mathbf{y} \mid \lambda)\, p(\lambda)\, d\lambda$$
4. The Marginal Likelihood — Full Derivation
Substituting the likelihood and prior:
$$p(\mathbf{y}) = \int_0^\infty \frac{\lambda^{S}\, e^{-n\lambda}}{\prod_{i=1}^n y_i!} \cdot \frac{\beta^\alpha}{\Gamma(\alpha)}\, \lambda^{\alpha-1} e^{-\beta\lambda}\, d\lambda$$
Pulling out all constants with respect to $\lambda$:
$$p(\mathbf{y}) = \frac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!} \int_0^\infty \lambda^{S} e^{-n\lambda} \cdot \lambda^{\alpha-1} e^{-\beta\lambda}\, d\lambda$$
Combining powers of $\lambda$ and the exponential terms:
$$p(\mathbf{y}) = \frac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!} \int_0^\infty \lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}\, d\lambda$$
Now we recognise this integral as a Gamma integral. Recall:
$$\int_0^\infty \lambda^{a-1} e^{-b\lambda}\, d\lambda = \frac{\Gamma(a)}{b^a}$$
Matching with $a = \alpha + S$ and $b = \beta + n$:
$$\int_0^\infty \lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}\, d\lambda = \frac{\Gamma(\alpha+S)}{(\beta+n)^{\alpha+S}}$$
Substituting back:
$$\boxed{p(\mathbf{y}) = \frac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!} \cdot \frac{\Gamma(\alpha + S)}{(\beta+n)^{\alpha+S}}}$$
5. Substituting into Bayes’ Theorem
Numerator:
$$p(\mathbf{y} \mid \lambda)\cdot p(\lambda) = \frac{\lambda^{S} e^{-n\lambda}}{\prod_{i=1}^n y_i!} \cdot \frac{\beta^\alpha}{\Gamma(\alpha)}\,\lambda^{\alpha-1} e^{-\beta\lambda}$$
$$= \frac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!}\cdot \lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}$$
Denominator:
$$p(\mathbf{y}) = \frac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!} \cdot \frac{\Gamma(\alpha+S)}{(\beta+n)^{\alpha+S}}$$
Dividing numerator by denominator:
$$p(\lambda \mid \mathbf{y}) = \frac{\dfrac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!}\cdot \lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}}{\dfrac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!} \cdot \dfrac{\Gamma(\alpha+S)}{(\beta+n)^{\alpha+S}}}$$
The term $\dfrac{\beta^\alpha}{\Gamma(\alpha)\prod_{i=1}^n y_i!}$ cancels from numerator and denominator:
$$p(\lambda \mid \mathbf{y}) = \frac{\lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}}{\dfrac{\Gamma(\alpha+S)}{(\beta+n)^{\alpha+S}}}$$
Bringing the denominator up:
$$p(\lambda \mid \mathbf{y}) = \frac{(\beta+n)^{\alpha+S}}{\Gamma(\alpha+S)}\cdot \lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}$$
6. Identifying the Posterior
Recalling the Gamma density $\text{Gamma}(a,\, b) \propto \lambda^{a-1} e^{-b\lambda}$ with normalising constant $\dfrac{b^a}{\Gamma(a)}$, we match:
$$a = \alpha + S = \alpha + \sum_{i=1}^n y_i, \qquad b = \beta + n$$
Therefore:
$$\boxed{p(\lambda \mid \mathbf{y}) = \frac{(\beta+n)^{\alpha+S}}{\Gamma(\alpha+S)}\,\lambda^{(\alpha+S)-1}\, e^{-(\beta+n)\lambda}}$$
$$\therefore \quad \lambda \mid \mathbf{y} \sim \text{Gamma}\left(\alpha + \sum_{i=1}^n y_i,\ \beta + n\right)$$
7. Interpretation of Updated Hyperparameters
Shape: The prior shape is $\alpha$ and the posterior shape is $\alpha + S$, where $S = \sum_{i=1}^n y_i$ is the total observed count. The data adds the observed counts directly onto the prior shape.
Rate: The prior rate is $\beta$ and the posterior rate is $\beta + n$. Each new observation adds one unit to the prior rate, reflecting the accumulation of $n$ new observations.
Posterior Mean: The prior mean is $\frac{\alpha}{\beta}$. After observing the data, the posterior mean becomes $\frac{(\alpha + S)}{(\beta + n)}$, which is the updated estimate of the Poisson rate $\lambda$.