PUBDEV-8899: Tweedie Variance Power Maximum Likelihood estimation

[tomasfryda]

https://h2oai.atlassian.net/browse/PUBDEV-8899

NOTE: If you see a comma in math expression and it not within a polygamma function ($\psi(a,b)$) then it's just happen to be misrendered \, (small space; 3/18 of a quad).

Implementation

Optimization

1. Bracket a region where $p$ maximizes likelihood ($\mathcal{L}$)

   • if response contains zeros ( $\implies p \in (1,2)$ ) use golden ratio search to find a small neighborhood where $p$ maximizes likelihood
   • if response does not contain zeros, look at $p = 3$
     • if $\mathcal{L}(p=2) > \mathcal{L(p=3)}$ search $p$ between 1 and 3 (note both $p=2$ and $p=3$ are cheap to compute - gamma and inverse gaussian distributions respectively)
     • else set $p_i = 2p_{i-1}$ until $\mathcal{L}(p_{i-1}) > \mathcal{L(p_i)}$ then set the search region to $[p_{i-2}, p_i]$

2. Once the bracket is small enough ($\Delta p = 0.1$) try to use Newton's method which should converge much faster

   • if the Newton's method steps out of the bracket 3 times, we will use only golden ratio search
   • if $p \in [1.95, 2.1]$ we won't attempt to use Newton's method since the gradient and hessian can get very noisy there

This optimization procedure is further optimized to speed up the convergence a bit:

• If iterCnt > 1: evaluate p from previous iteration, and p-0.01, p+0.01, if p has highest likelihood search only in [p-0.01, p+0.01]. This can eliminates multiple evaluations if p >> 2 and if it doesn't fail it eliminates ~5 evaluations for p in (1, 3] each iteration. It tends to fail only in first couple iterations.

Likelihood, gradient and hessian calculation

Likelihood is calculated by the best method available depending on the values of $y, \mu, \varphi$ using the $\xi\text{-based}$ heuristic proposed in section 8 in [1]. I also tried to use a heuristic from R's tweedie::dtweedie which improved results for $p=1.2$ but overall the results worsened and weren't as stable - starting Tweedie variance significantly influenced (±0.05) the estimated Tweedie var. power and at the same time it utilized the Inversion method more often for cases that took long time to estimate which increased the runtime ~10x.

I also tried several different combinations of those two heuristics but it seems that using the one from [1] is the best but it has troubles with $p &lt; 1.3$ and sometimes when $p$ is close to 2 it gets noisy (2, 2.1) ends up often in $-\infty$ so for this part I force the likelihood during the optimization to be estimated using the inversion method that provides more stable results and use the golden ratio search. Problem with this is that there are slight differences between the two methods that can add up to a difference that would influence the optimization procedure, for that reason the inversion method is forced to be used only if the bracket is small enough (lower bound > 1.95 && upper bound < 2.1).

Gradient and hessian is calculated only by the series method. This results in more robust likelihood calculation that can be used to correct the optimization procedure when using the Newton's method.

Details

What methods are used to estimate likelihood?

Depending on $p$, $y$, $\phi$ different methods are used for log likelihood estimation:
Let
$$\xi = \frac{\phi}{y^{2-p}}$$

If $p = 2$ use log likelihood of gamma distribution:

$$\log(p) = \begin{cases} -\infty \quad y = 0 \\ \frac{1}{\phi} log(\frac{1}{\phi\mu}) - \mathrm{logGamma}(\frac{1}{\phi}) + \log(y) (\frac{1}{\phi} - 1) + (-\frac{1}{\phi\mu} y)\quad y> 0\end{cases}$$

If $p = 3$ use inverse Gaussian distribution:

$$\log(p) = \begin{cases} -\infty \quad y = 0 \\ \frac{1}{2}\left(-\log(\phi\mu) - \log(2\pi) - 3 \log(\frac{y}{\mu}) - \frac{(\frac{y}{\mu} - 1)^2} { \phi\mu \frac{y}{\mu}}\right) - \log(\mu)\quad y> 0\end{cases}$$

If $p &lt; 2$ and $\xi \leq 0.01$ use Fourier inversion method[1] briefly described below.

If $p &gt; 2$ and $\xi \geq 1$ use Fourier inversion method.

else use the Series method[4] (also described in more detail below). If the series method fails (outputs NaN) try to use Fourier inversion method.

If Series method and Fourier inversion method failed or if the Fourier inversion method was chosen based on the $\xi$ criterium and it failed, estimate the log likelihood using the Saddlepoint approximation (also described below).

During the optimization we sometimes want to force usage of the Fourier inversion method, this happens when the whole bracket is within [1.95, 2.1].

Series Method [4]

NOTE: $\psi(x)$ is digamma function and $\psi(1,x)$ is a trigamma function

1<p<2

Log likelihood for Y = 0

$\log \mathrm{P}(y=0) = -\frac{w\mu^{2-p}}{\phi(2-p)}$

first derivative with respect to $p$:
$\displaystyle \frac{\partial \log \mathrm{P}(y=0)}{\partial p} = -\frac{\frac{\mu^{-p + 2} w \log\left(\mu\right)}{p - 2} + \frac{\mu^{-p + 2} w}{{\left(p - 2\right)}^{2}}}{\phi}$

second derivative with respect to $p$:
$\displaystyle \frac{\partial \log \mathrm{P}(y=0)}{\partial^2 p} = \frac{\frac{\mu^{-p + 2} w \log\left(\mu\right)^{2}}{p - 2} + \frac{2 , \mu^{-p + 2} w \log\left(\mu\right)}{{\left(p - 2\right)}^{2}} + \frac{2 , \mu^{-p + 2} w}{{\left(p - 2\right)}^{3}}}{\phi}$

Log likelihood for Y > 0

$$W_j = \frac{{\left(p - 1\right)}^{\frac{j p}{p - 1} - \frac{2 , j}{p - 1}} w^{\frac{j}{p - 1}} y^{-\frac{j p}{p - 1} + \frac{2 , j}{p - 1}}}{{\left(-p + 2\right)}^{j} \phi^{\frac{j}{p - 1}} \Gamma\left(j + 1\right) \Gamma\left(-\frac{j p - 2 , j}{p - 1}\right)}$$

$$\log\mathrm{P}(y) = -\log\left(y\right) + \log\left({\sum_{j=1}^{+\infty} W_j}\right) +\frac{w}{\phi} (y\theta - \kappa(\theta))$$

First derivative with respect to $p$:
$$\frac{\partial W_j}{\partial p} = \frac{{\left(p {\left(\log\left(\frac{{\left(p - 1\right)} \phi}{w y}\right) - 2\right)} + {\left(p - 2\right)} \psi\left(-\frac{j {\left(p - 2\right)}}{p - 1}\right) + 2 , \log\left(\frac{w y}{{\left(p - 1\right)} \phi}\right) + 3\right)} \left({\left(p - 1\right)}^{\frac{p - 2}{p - 1}} w^{\left(\frac{1}{p - 1}\right)}\right)^{j} j}{{\left(p^{3} - 4 , p^{2} + 5 , p - 2\right)} \left(-{\left(p - 2\right)} \phi^{\left(\frac{1}{p - 1}\right)} y^{\frac{p - 2}{p - 1}}\right)^{j} \Gamma\left(j + 1\right) \Gamma\left(-\frac{j {\left(p - 2\right)}}{p - 1}\right)}$$

$$\frac{\partial \log\mathrm{P}(y)}{\partial p} = \dfrac{{\sum}{j=1}^{+\infty} \cfrac{\partial W_j}{\partial p}}{{\sum}{j=1}^{+\infty} W_j } + \frac{{\left(\frac{\mu^{-p + 1} y \log\left(\mu\right)}{p - 1} - \frac{\mu^{-p + 2} \log\left(\mu\right)}{p - 2} + \frac{\mu^{-p + 1} y}{{\left(p - 1\right)}^{2}} - \frac{\mu^{-p + 2}}{{\left(p - 2\right)}^{2}}\right)} w}{\phi}$$

Second derivative with respect to $p$: (breaks math rendering in github)

\frac{\partial \log\mathrm{P}(y)}{\partial^2 p} =-\frac{{\left({\left(p {\left(\psi\left(-\frac{j p - 2 \, j}{p - 1}\right) - 2\right)} + p \log\left(\frac{{\left(p - 1\right)} \phi}{w y}\right) - 2 \, \psi\left(-\frac{j p - 2 \, j}{p - 1}\right) + 3\right)} j {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2 \, p}{p - 1}} {\left(p - 2\right)} \log\left(y\right) - {\left(2 \, p^{4} - 11 \, p^{3} + 22 \, p^{2} - {\left(p^{4} - 6 \, p^{3} + 13 \, p^{2} - 12 \, p + 4\right)} \psi\left(-\frac{j p - 2 \, j}{p - 1}\right) - 19 \, p + 6\right)} j {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2}{p - 1}} \log\left(w\right) - {\left({\left(p^{2} - 4 \, p + 4\right)} j \psi\left(-\frac{j p - 2 \, j}{p - 1}\right)^{2} - 2 \, {\left(j \log\left(p - 1\right) - 2 \, j + 11\right)} p^{2} + 5 \, p^{3} - {\left(2 \, p^{2} - {\left(p^{2} - 4 \, p + 4\right)} \psi\left(-\frac{j p - 2 \, j}{p - 1}\right) - 7 \, p + 6\right)} j \log\left(\phi\right) - {\left(p^{2} - 4 \, p + 4\right)} j \psi\left(1, -\frac{j p - 2 \, j}{p - 1}\right) + {\left(7 \, j \log\left(p - 1\right) - 12 \, j + 32\right)} p - 6 \, j \log\left(p - 1\right) + {\left({\left(j \log\left(p - 1\right) - 4 \, j + 10\right)} p^{2} - 2 \, p^{3} - 2 \, {\left(2 \, j \log\left(p - 1\right) - 7 \, j + 8\right)} p + 4 \, j \log\left(p - 1\right) - 12 \, j + 8\right)} \psi\left(-\frac{j p - 2 \, j}{p - 1}\right) + 9 \, j - 15\right)} {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2 \, p}{p - 1}} + {\left({\left(p^{3} - 4 \, p^{2} + 5 \, p - 2\right)} j {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2}{p - 1}} p \log\left(w\right) - {\left(j {\left(p - 2\right)} p \log\left(\phi\right) + j {\left(p - 2\right)} p \psi\left(-\frac{j p - 2 \, j}{p - 1}\right) + {\left(j \log\left(p - 1\right) - 2 \, j + 4\right)} p^{2} - 2 \, p^{3} - j p {\left(2 \, \log\left(p - 1\right) - 3\right)} - 2\right)} {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2 \, p}{p - 1}}\right)} \log\left(\frac{{\left(p - 1\right)} \phi}{w y}\right) + 2 \, {\left({\left(p^{3} - 4 \, p^{2} + 5 \, p - 2\right)} j {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2}{p - 1}} \log\left(w\right) + j {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2 \, p}{p - 1}} {\left(p - 2\right)} \log\left(y\right) - {\left(j {\left(p - 2\right)} \log\left(\phi\right) + j {\left(p - 2\right)} \psi\left(-\frac{j p - 2 \, j}{p - 1}\right) + {\left(j \log\left(p - 1\right) - 2 \, j + 8\right)} p - 3 \, p^{2} - 2 \, j \log\left(p - 1\right) + 3 \, j - 5\right)} {\left(p - 1\right)}^{\frac{j p}{p - 1} + \frac{2 \, p}{p - 1}}\right)} \log\left(\frac{w y}{{\left(p - 1\right)} \phi}\right)\right)} \left(-1\right)^{j} w^{\frac{j}{p - 1}} y^{\frac{2 \, j}{p - 1}}}{{\left(p^{6} - 8 \, p^{5} + 26 \, p^{4} - 44 \, p^{3} + 41 \, p^{2} - 20 \, p + 4\right)} {\left(j - 1\right)} {\left(p - 1\right)}^{\frac{2 \, j}{p - 1} + \frac{2 \, p}{p - 1}} {\left(p - 2\right)}^{j} \phi^{\frac{j}{p - 1}} y^{\frac{j p}{p - 1}} \left(j - 2\right)! \left(-\frac{{\left(j + 1\right)} p - 2 \, j - 1}{p - 1}\right)!}-\frac{{\left(\frac{\mu^{-p + 1} y \log\left(\mu\right)^{2}}{p - 1} - \frac{\mu^{-p + 2} \log\left(\mu\right)^{2}}{p - 2} + \frac{2 \, \mu^{-p + 1} y \log\left(\mu\right)}{{\left(p - 1\right)}^{2}} - \frac{2 \, \mu^{-p + 2} \log\left(\mu\right)}{{\left(p - 2\right)}^{2}} + \frac{2 \, \mu^{-p + 1} y}{{\left(p - 1\right)}^{3}} - \frac{2 \, \mu^{-p + 2}}{{\left(p - 2\right)}^{3}}\right)} w}{\phi}

p > 2

Undefined for y = 0.

For y > 0:

Log likelihood:
$$V_k = \frac{\left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 , \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)}$$
$$\log\mathrm{P}(y) = - \log(\pi y) + \log\left(\sum_{k=1}^{+\infty}V_k\right)+\frac{w}{\phi} (y\theta - \kappa(\theta))$$

First derivative with respect to $p$:

$$\frac{\partial V_k}{\partial p} = \frac{{\left(\frac{{\left(\frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}} + \frac{k \log\left(p - 1\right)}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 , \pi k}{p - 1}\right)}{\phi^{\frac{k}{p - 1}}} - \frac{\pi k \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 , \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)}{{\left(p - 1\right)}^{2} \phi^{\frac{k}{p - 1}}} - \frac{k \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 , \pi k}{p - 1}\right)}{{\left(p - 2\right)} \phi^{\frac{k}{p - 1}}} + \frac{k {\left(\log\left(\phi\right) - \log\left(w\right) - \log\left(y\right) + \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 , \pi k}{p - 1}\right)}{{\left(p - 1\right)}^{2} \phi^{\frac{k}{p - 1}}}\right)} \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}}}{\left({\left(p - 2\right)} w^{\frac{p - 2}{p - 1} - 1} y^{\frac{p - 2}{p - 1}}\right)^{k} \Gamma\left(k + 1\right)}$$

$$\frac{\partial \log\mathrm{P}(y)}{\partial p} = \dfrac{{\sum}{k=1}^{+\infty} \cfrac{\partial V_k}{\partial p}}{{\sum}{k=1}^{+\infty} V_k} +\frac{{\left(\frac{\mu^{-p + 1} y \log\left(\mu\right)}{p - 1} - \frac{\mu^{-p + 2} \log\left(\mu\right)}{p - 2} + \frac{\mu^{-p + 1} y}{{\left(p - 1\right)}^{2}} - \frac{\mu^{-p + 2}}{{\left(p - 2\right)}^{2}}\right)} w}{\phi}$$

Second derivative with respect to $p$: (breaks math rendering in github)

\frac{\partial \log\mathrm{P}(y)}{\partial^2 p} = \sum_{k=1}^{+\infty}\frac{\left(-1\right)^{k} k^{2} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right)^{2} \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k^{2} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right) \log\left(w\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{\left(-1\right)^{k} k^{2} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right)^{2} \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right) \log\left(y\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right) \log\left(y\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right)^{2} \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{4 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{\left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right)^{2} \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k^{2} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{k - 1} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, {\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k^{2} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{k - 1} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, {\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{p - 1} - \frac{p - 2}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{k - 1} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, {\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{\left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)^{2} \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{{\left(p - 1\right)}^{2}} - \frac{p - 2}{{\left(p - 1\right)}^{3}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(\phi\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{1}{{\left(p - 1\right)}^{2}} - \frac{p - 2}{{\left(p - 1\right)}^{3}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(w\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{k}{{\left(p - 1\right)}^{2}} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{3}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \log\left(y\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{2 \, k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{k - 1} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{\left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)}^{2} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \psi\left(1, \frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, {\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)} {\left(p - 2\right)}^{k - 1} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, {\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, {\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \left(-1\right)^{k} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{k - 1} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{\left(-1\right)^{k} {\left(k - 1\right)} k {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{k - 2} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{2 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{{\left({\left(\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)} \log\left(p - 1\right) + \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}\right)}^{2} \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{\left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(\frac{\pi k}{p - 1} - \frac{\pi k p}{{\left(p - 1\right)}^{2}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{2}}\right)}^{2} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{k}{{\left(p - 1\right)}^{2}} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{3}}\right)} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \psi\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} {\left(\frac{\pi k}{{\left(p - 1\right)}^{2}} - \frac{\pi k p}{{\left(p - 1\right)}^{3}} + \frac{2 \, \pi k}{{\left(p - 1\right)}^{3}}\right)} \cos\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right) \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} + \frac{2 \, \left(-1\right)^{k} k^{2} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} {\left(p - 2\right)}^{2 \, k - 2} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{3 \, k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} - \frac{{\left(2 \, {\left(\frac{k}{{\left(p - 1\right)}^{2}} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{3}}\right)} \log\left(p - 1\right) - \frac{\frac{k}{p - 1} - \frac{k {\left(p - 2\right)}}{{\left(p - 1\right)}^{2}}}{p - 1} - \frac{k}{{\left(p - 1\right)}^{2}} + \frac{2 \, k {\left(p - 2\right)}}{{\left(p - 1\right)}^{3}}\right)} \left(-1\right)^{k} {\left(p - 1\right)}^{\frac{k {\left(p - 2\right)}}{p - 1}} \phi^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} \Gamma\left(\frac{k {\left(p - 2\right)}}{p - 1} + 1\right) \sin\left(-\frac{\pi k p}{p - 1} + \frac{2 \, \pi k}{p - 1}\right)}{{\left(p - 2\right)}^{k} w^{k {\left(\frac{p - 2}{p - 1} - 1\right)}} y^{\frac{k {\left(p - 2\right)}}{p - 1}} \Gamma\left(k + 1\right)} -\frac{{\left(\frac{\mu^{-p + 1} y \log\left(\mu\right)^{2}}{p - 1} - \frac{\mu^{-p + 2} \log\left(\mu\right)^{2}}{p - 2} + \frac{2 \, \mu^{-p + 1} y \log\left(\mu\right)}{{\left(p - 1\right)}^{2}} - \frac{2 \, \mu^{-p + 2} \log\left(\mu\right)}{{\left(p - 2\right)}^{2}} + \frac{2 \, \mu^{-p + 1} y}{{\left(p - 1\right)}^{3}} - \frac{2 \, \mu^{-p + 2}}{{\left(p - 2\right)}^{3}}\right)} w}{\phi}

Fourier Inversion method

Implemented as method 3 described in [1]:

(3) $$f(y; \mu, \phi) = b(y, \phi)\exp\left\{-\frac{1}{2\phi} d(y, \mu) \right\}$$
where $d(y, \mu)$ is deviance.

(6) $$f(y;\mu, \phi) = cf(cy; c\mu, c^{2-p}\phi)$$

Method 3: Rescale y to 1 and evaluate b. Compute I = a(1, φ/y2−p) = f (1; 1, φ/y2−p) by Fourier inversion, use b(y, φ) = I/y, and substitute into (3). This method uses (6) with c = 1/y.

Let

$$\alpha = \frac{2-p}{1-p}$$

$$\zeta = \tan^{-1}\left\{ (1-p)t\phi\right\}$$

$$\mathcal{R}k(t) = \frac{1}{\phi(2-p)} \frac{\cos(\zeta\alpha)}{(\cos\zeta)^\alpha} - \frac{1}{\phi(2-p)}$$

$$\mathcal{I}k(t) = \frac{1}{\phi(2-p)} \frac{\sin(\zeta\alpha)}{(\cos\zeta)^\alpha} - yt$$

1<p<2

Two conditional densities (y=0, y>0) merged together

For y=0:

$$\pi_0 = P(y=0|\mu=1) = \exp\left\{\frac{-1}{\phi(2-p)}\right\}$$

For y > 0:

$$a^{*}(y, \phi) = \frac{1}{\pi(1-\pi_0)}\int_0^\infty \left\{ \exp \mathcal{R}k(t)\cos\mathcal{I}k(t) - \pi_0 \cos(ty)\right\}\mathrm{d}t$$

Resulting formula:

$$a(y, \phi) = a^{*}(y, \phi)(1-\pi_0)$$

p > 2

$$a(y,\phi) = \frac{1}{\pi} \int_0^{\infty} \exp\mathcal{R}k(t)\cos\mathcal{I}k(t)\mathrm{d}t $$

Saddlepoint approximation method [3]

(Used only if previous methods produce NaN)

$$p(y;\mu,\phi) = (2\pi\phi y ^p)^{-1/2}\exp\left\{ \frac{-d(y,\mu)}{2\phi} \right\} \left\{1+\mathcal{O}(\phi)\right\} $$

Log likelihood:
$$\log(p) = -0.5 (\log(2 \pi \phi) + p \log(y)) + (-\mathrm{d}(y, \mu) / (2 \phi))$$

Usage

For all cases I assume the dispersion method is set to ML

• with fix_tweedie_variance_power = True and fix_dispersion_parameter = True|False I don't optimize either but the likelihood is calculated properly so it can be used
  for external optimization of both p and phi(using the init_dispersion_parameter), e.g., Bayesian optimisation with Gaussian process with squared exponential kernel could work (having nice smooth estimates could mitigate the noisyness of the loglikelihood; but I did not try this approach) Using the old calculation as the new one caused NonSPDMatrixError in pyunit_pubdev_8685_tweedie_dispersion_factor_exceed2.

• with fix_tweedie_variance_power = False and fix_dispersion_parameter = True - it will optimize the Tweedie variance power as described in the Optimization section of this PR. As dispersion (phi, $\varphi$) is used value set in init_dispersion_parameter

• with fix_tweedie_variance_power = False and fix_dispersion_parameter = False - will in future optimize both Tweedie variance power (p) and dispersion parameter ($\varphi$)

TODO

• test on more datasets
• test with link power estimation doesn't work so I removed the related code

Current Limitations

When response contains zeros, Tweedie Var. Power (p) will be estimated only between 1 and 2

When skipping likelihood calculation for zeros with $p&gt;=2$, we can see an artifact created by this around 2:

If we use the same dataset but filter out the rows with $\text{response} <= 0 $ we could estimate some $p$ but the value is different then using just the proper likelihood where the log likelihood of response 0 or lower turns to $-\infty$ for $p &gt;= 2$:

Properly calculated log likelihood:

Not estimating Poisson

As can be seen from the previous plot $p = 1$ which corresponds to poisson (when $\varphi = 1$) Poisson looks very likely but there are 2 problems:

• can likelihood calculated from pmf be compared to likelihood calculated to pdf?
• neighborhood around $p=1$ is not very stable for likelihood calculation

Problems

• when p gets close to 1, the region gets more noisy. Generally $p &gt; 1.3$ doesn't seem to be affected.
• duration - for some combinations of $y, p, \varphi$ the likelihood is hard to compute and it takes some time. This could be partially improved by using an approximation scheme using Chebyshev polynomials[3] (I believe this is also implemented in R's tweedie::dtweedie.interp) but I wasn't able to try it yet.
• (potentially) not supporting canonical link function - [2; bottom of the page 124] describes the canonical link function as $-\frac{1}{p-1} \mu ^{-(p-1)}$ and in H2O we don't divide by the $-(p-1)$.
• for some cases the optimization can fail due to NonSPDMatrixException

References

[1] DUNN, Peter and SMYTH, Gordon, 2008. Evaluation of Tweedie exponential dispersion model densities by Fourier inversion.
[2] OHLSSON, Esbjörn and JOHANSSON, Björn, 2006. Exact Credibility and Tweedie Models. ASTIN Bulletin: The Journal of the IAA. Online. May 2006. Vol. 36, no. 1, p. 121–133. link
[3] DUNN, Peter and SMYTH, Gordon, 2001. Tweedie Family Densities: Methods of Evaluation.
[4] DUNN, Peter K. and SMYTH, Gordon K., 2005. Series evaluation of Tweedie exponential dispersion model densities.

Attachments

Sage notebooks.zip

[tomasfryda]

I'm not sure about this behavior. I think it might be better to end up with some model than none, but the user has to be warned that the model didn't finish training properly. If we decide to keep it, we should probably improve the message.

[wendycwong]

When you nave Non SPD matrix, it cannot be avoided. It means the dataset is bad the the gram matrix cannot be inverted because it is not SPD (semi-definite positive definite).

[wendycwong]

Sure, this is fine.

[wendycwong]

I tried your code and it works really well for most cases. It has problem for power = 2.1 and power = 2.6. I don't think setting the bounds of power to between 1 and 2 is a good idea. I know the reason you choose to do. I have attached several datasets that the algorithm failed to arrive at a good estimate. You can access them here:

s3://h2o-public-test-data/smalldata/glm_test/tweedie_p2p1_disp3_5Cols_10kRows_est1p89.csv
s3://h2o-public-test-data/smalldata/glm_test/tweedie_p2p6_disp2_5Cols_10krows_est1p94.csv

These datasets do contains y <= 0 and the algo somehow locks it within the 1 and 2 range and will not consider the outside range at all.

I think a good remedy here is to put into the documentation that if the user believes that the variance power is > 2, they should remove all rows that have resp <= 0. However this does not seem to work for s3://h2o-public-test-data/smalldata/glm_test/tweedie_p2p6_disp2_5Cols_10krows_est1p94.csv. I will slack you my code.

[tomasfryda]

Thank you @wendycwong for the thorough testing.

I have tried both datasets with response above 0 training_data[training_data[Y] > 0, :] and it seems to work(1), however, when I don't remove the zeros it doesn't converge to the desired values(2).

I can remove the logic that sets the range to (1,2) when there are zeros (or lower numbers in the response), however, it will not change anything. For $p \geq 2$ zero in the response will end up with log likelihood = -infinity. So it will be locked to $p \in (1, 2)$ anyway.

I have tried ignoring zeros in the response but that increases the likelihood of $p \geq 2$ so even if the true value would be e.g. 1.8 it could be estimated to 2.5 because this way we would be comparing likelihoods on two different datasets.

I wonder if the difference can be just from the definition - if you want to have non-zero chance of having a zero in the response you have to allocate for it some non-zero probability so there is less probability for the pdf since:
$p \in (1,2)$ :
$\pi_0 + \int_{x\in (0, +\infty)} P(x|\mu,\phi) dx = 1\qquad \pi_0 = \exp \left\{-\frac{\mu^{2-p}}{\phi(2-p)}\right\}$
$p \geq 2$ :
$\int_{x\in (0, +\infty)} P(x|\mu, \phi) dx = 1$

And this tiny difference is added with every single row so it might get big enough to influence the results.

(1)

estimated variance power: 2.144064040064248, true variance power: 2.1
estimated variance power: 2.67437207518084, true variance power: 2.6

(2)

estimated variance power: 1.8851764629473313, true variance power: 2.1
estimated variance power: 1.9477466113563873, true variance power: 2.6

I looked more in detail on the case with true p value = 2.6:
I ran tweedie::tweedie.profile in R for both whole dataset and the dataset without zeros:

Note: vertical lines denote the optimum found out by H2O, all the points are, however, calculated only using R's package tweedie. Y axis is log likelihood, X axis is the tweedie variance power.

With zeros:

Without zeros:

Combined - filled points are from the dataset with zeros:

[wendycwong]

Can you separate out tests that take longer time and move them to ***_large.py? Also, the dataset they use should be found in bigdata/laptop/glm_test as well?
This will help reduce our test timeout.