Exact Likelihood Evaluation by Series

Why the Normalizing Series Matters

For a Tweedie model with power index between one and two, a response can be zero or positive. For a positive observation, its density is not determined by the unit deviance alone: it also contains a normalizing factor that depends on the response, mean, dispersion, and power index. Exact likelihood evaluation therefore requires evaluating an infinite series. The series has positive terms, which makes both its derivation and reliable numerical evaluation particularly transparent.

Write $\mu>0$ for the mean, $\phi>0$ for the dispersion, and $p$ for the power index. Throughout, assume $1<p<2$. The compound Poisson–Gamma representation gives a Poisson count of Gamma increments; conditional on a positive total, summing over all possible positive counts gives the density series below.

Definition: Set $\alpha=(2-p)/(p-1)$, $\lambda=\mu^{2-p}/\{\phi(2-p)\}$, and $\gamma=\phi(p-1)\mu^{p-1}$. For $y>0$, the exact Tweedie density is $$ f(y;\mu,\phi,p)=\frac{\exp(-\lambda-y/\gamma)}{y} \sum_{j=1}^{\infty} \frac{\lambda^j(y/\gamma)^{j\alpha}}{j!\,\Gamma(j\alpha)}. $$ Here $\Gamma$ is the Gamma function. The zero outcome is not evaluated with this density: it has probability $P(Y=0)=\exp(-\lambda)$.

The index $j$ counts the number of Gamma increments. In particular, the terms in the series are nonnegative, and no cancellation between positive and negative contributions is involved. That fact is useful when constructing a truncation bound. It does not, however, mean that direct computation of the terms is always safe: at some parameter values the terms can be extremely small or large.

Deriving the Positive-Response Series

In the compound Poisson–Gamma representation, let $N\sim\operatorname{Pois}(\lambda)$, and let the increments have Gamma shape $\alpha$ and scale $\gamma$, independently of $N$. Conditional on $N=j\geq1$, the sum has a Gamma distribution with shape $j\alpha$ and scale $\gamma$. Its density at $y>0$ is $\{y^{j\alpha-1}\exp(-y/\gamma)\}/\{\Gamma(j\alpha)\gamma^{j\alpha}\}$. Multiplying this conditional density by $P(N=j)=\exp(-\lambda)\lambda^j/j!$, then summing over $j\geq1$, gives

\begin{aligned} f(y;\mu,\phi,p) &=\sum_{j=1}^{\infty}e^{-\lambda}\frac{\lambda^j}{j!} \frac{y^{j\alpha-1}e^{-y/\gamma}}{\Gamma(j\alpha)\gamma^{j\alpha}}\\ &=\frac{e^{-\lambda-y/\gamma}}{y} \sum_{j=1}^{\infty} \frac{\lambda^j(y/\gamma)^{j\alpha}}{j!\Gamma(j\alpha)}. \end{aligned}

This is a mixture-density calculation: for $y>0$, the component densities are weighted by their Poisson probabilities. The event $N=0$ contributes the point mass at zero instead, which explains why it is absent from the positive-response sum.

Term Ratios and a Certified Truncation Bound

For computation, isolate the series term $t_j=\lambda^j(y/\gamma)^{j\alpha}/\{j!\Gamma(j\alpha)\}$, and put $x=\lambda(y/\gamma)^\alpha$. Then $t_j=x^j/\{j!\Gamma(j\alpha)\}$. The ratio of successive terms is

R_j=\frac{t_{j+1}}{t_j} =\frac{x}{j+1}\frac{\Gamma(j\alpha)}{\Gamma((j+1)\alpha)}.

These ratios give a practical stopping rule, but a stopping rule based only on a small last term is not by itself a guarantee about the omitted tail. The following result supplies such a guarantee.

Theorem: For fixed $x>0$ and $\alpha>0$, the series $\sum_{j=1}^{\infty}t_j$ converges. If $R_{m+1}<1$, then the omitted tail after retaining terms $1,\ldots,m$ satisfies $$ 0<\sum_{j=m+1}^{\infty}t_j \leq \frac{t_{m+1}}{1-R_{m+1}}. $$

Proof. The function $\log\Gamma(u)$ is strictly convex for $u>0$. Thus the increment $\log\Gamma(u+\alpha)-\log\Gamma(u)$ increases with $u$. Consequently, $\Gamma((j+1)\alpha)/\Gamma(j\alpha)$ increases with $j$, so its reciprocal decreases. Since $1/(j+1)$ also decreases, $R_j$ decreases with $j$.

Stirling’s formula gives $\Gamma(u+\alpha)/\Gamma(u)\sim u^\alpha$ as $u$ tends to infinity. Applying this with $u=j\alpha$ shows that $R_j\sim x/\{(j+1)(j\alpha)^\alpha\}$, which tends to zero. The ratio test therefore proves convergence. If $R_{m+1}<1$, then for every $k\geq0$, monotonicity gives $t_{m+2+k}\leq t_{m+1}R_{m+1}^{k+1}$. Summing these geometric upper bounds yields \[ \sum_{j=m+1}^{\infty}t_j \leq t_{m+1}\sum_{k=0}^{\infty}R_{m+1}^{k} =\frac{t_{m+1}}{1-R_{m+1}}. \] All terms are positive, so the tail is also strictly positive. $\square$

To translate this into a density bound, multiply the series-tail bound by $e^{-\lambda-y/\gamma}/y$. To obtain a relative bound on the series, divide the bound by the computed partial sum $S_m=\sum_{j=1}^{m}t_j$. For the log-density, if the relative tail bound is $b$, then the log-density error from using $S_m$ is at most $\log(1+b)\leq b$.

Worked Calculations

Worked Example: Six Terms and a Certified Error Bound

Take $p=1.5$, $\mu=1$, $\phi=1$, and $y=1$. Then $\alpha=1$, $\lambda=2$, $\gamma=0.5$, and $x=\lambda(y/\gamma)^\alpha=4$. The series terms are $t_j=4^j/\{j!(j-1)!\}$. In order, the first six are

$t_1=4$, $t_2=8$, $t_3=64/12=16/3$, $t_4=256/144=16/9$, $t_5=1024/2880=16/45$, and $t_6=4096/86400=0.0474074074$. Therefore $S_6=19.5140740741$.

The next term is $t_7=16384/3628800\approx0.004514991$. Here $R_j=4/\{j(j+1)\}$, so the ratios decrease and $R_7=1/14$. The omitted tail is consequently no larger than $t_7/(1-1/14)\approx0.0048623$. Since the density prefactor is $e^{-4}/1$, the six-term density approximation is $e^{-4}S_6\approx0.3574$, with an absolute truncation error less than $e^{-4}(0.0048623)\approx0.0000891$. The relative error in the series is less than $0.0048623/19.5140741$, or about $0.000250$.

Worked Example: A Noninteger Power in the Series Terms

Take $p=4/3$, $\mu=1$, and $\phi=1.5$. These values give $\alpha=2$, $\lambda=1$, and $\gamma=0.5$. At $y=0.5$, $y/\gamma=1$, so $x=1$ and $t_j=1/\{j!\Gamma(2j)\}$. The first four terms are

$t_1=1$, $t_2=1/12\approx0.0833333$, $t_3=1/720\approx0.00138889$, and $t_4=1/120960\approx0.00000826720$. Thus $S_4\approx1.08473049$. In this case $R_j=1/\{(j+1)(2j)(2j+1)\}$. Therefore $t_5=t_4/360\approx2.296\times10^{-8}$, and $R_5=1/660$. The tail after four terms is bounded by $t_5/(1-1/660)$, which is less than $2.30\times10^{-8}$.

The prefactor is $e^{-1-1}/0.5=2e^{-2}$. Hence the four-term density approximation is $2e^{-2}S_4\approx0.2936$, and its absolute truncation error is less than $2e^{-2}(2.30\times10^{-8})$, or $3.2\times10^{-9}$. The bound is small because the term ratios are already much less than one.

Log-Sum Evaluation and Numerical Checks

Directly forming $t_j$ can overflow or underflow, particularly when factorials and Gamma functions are large. Instead calculate the log terms $\ell_j=j\log x-\log\Gamma(j+1)-\log\Gamma(j\alpha)$. For a finite set of terms, let $M=\max_{1\leq j\leq m}\ell_j$. Then

\log S_m=M+\log\left(\sum_{j=1}^{m}\exp(\ell_j-M)\right).

Each exponent in the inner sum is at most zero, which avoids overflowing exponentials. The approximate log density is $-\lambda-y/\gamma-\log y+\log S_m$. Evaluate $\log\Gamma$ with a reliable log-Gamma routine rather than first evaluating $\Gamma$ and then taking its logarithm.

Worked Example: Computing the Log Density from Log Terms

Continue with $p=1.5$, $\mu=1$, $\phi=1$, and $y=1$, for which $x=4$. Retaining six terms gives $S_6=19.5140740741$. The largest term is $t_2=8$, so use $M=\log 8$. The shifted terms are $t_j/8$; their sum is

$0.5+1+0.666666667+0.222222222+0.044444444+0.005925926=2.439259259$.

Thus $\log S_6=\log 8+\log(2.439259259)\approx2.97114$, and the six-term log-density approximation is $-4+\log S_6\approx-1.02886$. The relative tail bound from the first example is less than $0.000250$; hence the difference between this approximate log density and the exact log density is at most $\log(1+0.000250)<0.000250$. The bound concerns series truncation, not floating-point rounding, so high-precision arithmetic may still be needed when very tight accuracy is required.

A reliable implementation should check accuracy rather than treating a plausible decimal output as evidence of convergence. Compute terms by the recurrence $t_{j+1}=t_jR_j$, and compare selected terms with values obtained from $\ell_j$ and log-Gamma calculations. Increase $m$ until the geometric tail bound is below the desired absolute or relative tolerance. Repeating the calculation at a larger $m$ is a useful diagnostic, but the explicit tail bound is what certifies the omitted-series error.

Common Pitfalls

Dropping part of a term when factoring out the density prefactor. Factoring out $e^{-\lambda-y/\gamma}/y$ leaves $t_j=\lambda^j(y/\gamma)^{j\alpha}/\{j!\Gamma(j\alpha)\}$. For the first example, this leaves $4^j/\{j!(j-1)!\}$, not $1/\{j!(j-1)!\}$.
Using the last term as if it were the whole remaining error. The omitted tail contains every later term. Use $t_{m+1}/(1-R_{m+1})$ when $R_{m+1}<1$, not merely $t_{m+1}$.
Using a positive-response density at zero. Zero is a point mass with probability $e^{-\lambda}$. It is not obtained by substituting $y=0$ into the positive-response density formula.
Confusing truncation accuracy with total numerical accuracy. A valid tail bound controls omitted terms, but rounding, inaccurate Gamma evaluations, or overflow can still affect the computed sum.

Key takeaway: For $1<p<2$, exact positive-response likelihood evaluation reduces to a positive Gamma-function series. Compute its terms on the log scale when needed, and certify truncation with a geometric tail bound based on successive-term ratios.

Check Your Understanding

Use the series, ratio, and error-bound methods developed above.

For $p=1.5$, $\mu=1$, $\phi=1$, and $y=1$, write down $x$ and $t_3$.
Why does the zero outcome contribute a probability mass rather than a value from the positive-response density?
Suppose the terms decrease from $t_{m+1}$ onward and $R_{m+1}=0.2$. What upper bound does the theorem give for the tail after $m$ terms?
Explain how subtracting the largest log term before exponentiating helps stabilize a finite series calculation.
What does a geometric truncation bound certify, and what sources of numerical error does it not automatically certify?