How to approximate this expression that involve trigamma functions: $\psi ^{(1)}(n+1)-\psi ^{(1)}(n+m+2)$

Question

I have come up with this expression for the standard deviation for a certain distribution:

$$\psi ^{(1)}(n+1)-\psi ^{(1)}(n+m+2)$$

where $n$ and $m$ are integers and $\psi^{(1)}(x)$ is a trigamma function.

For the software engineering project I am involved with calculation of the trigamma function causes technical problems. How this expression can be approximated by something easier to calculate for large $m$ and $n$?

One expression I arrived at is

$$ \frac{(n+1)^3 (m+1)}{(n+m+2)^4 (n+m+3)} $$

The way I've arrived there is quite convoluted, and the expression most definitely may not be optimal.

Perhaps there is a better approximation?

Answer 1

Don't know if it would be useful, but:

The trigamma function satisfies the following recurrence $$\psi^{(1)}(x+1)=\psi^{(1)}(x)-\frac{1}{x^2}.$$ Using this for the second term in your expression, $$\psi^{(1)}(n+1)-\lbrack\psi^{(1)}(n+m+1)-\frac{1}{(n+m+1)^2}\rbrack$$ Repeat this until we get $\psi^{(1)}(n+1)$: $$\psi^{(1)}(n+1)-[\psi^{(1)}(n+1)-\sum_{k=1}^{m+1}\frac{1}{(n+k)^2}]$$ $$=\sum_{k=1}^{m+1}\frac{1}{(n+k)^2}.$$

From here I'm sure there would be a lot of ways to compute or approximate, but I don't know them or which are the best. If you do, go for it. I will do a little more but it isn't necessarily a good method: $$\sum_{k=1}^{m+1}\frac{1}{(n+k)^2} = \sum_{k=1}^{n+m+1}\frac{1}{k^2} - \sum_{k=1}^{n}\frac{1}{k^2}$$ $$=\sum_{k=1}^{\infty}\frac{1}{k^2} - \sum_{k=n+m+2}^{\infty}\frac{1}{k^2}- \sum_{k=1}^{n}\frac{1}{k^2}$$ $$=\frac{\pi^2}{6} - a_{n+m+1} - s_{n}$$ where $s_k$ is the $k$-th partial sum and $a_k$ is the $k$-th remainder. You can find estimates of both, for example, $$\frac{1}{n+m+2}<a_{n+m+1}<\frac{1}{n+m+1}, \space\space s_{n}\leq 2-\frac{1}{n}$$but I don't know if those will do. Hope this helps.

Answer 2

By Taylor, when $p$ is large $$\psi ^{(1)}(p)=\frac{1}{p}+\frac{1}{2 p^2}+\frac{1}{6 p^3}-\frac{1}{30 p^5}+\frac{1}{42 p^7}+O\left(\frac{1}{p^9}\right)$$ from wich the successive approximations

$$A_2=\frac{m+1}{2 (n+1)^2 (m+n+2)^2}P_2$$ with $$P_2=2 n^2+(2 m+8) n+(3 m+7)$$ and $$A_3=\frac{m+1}{6 (n+1)^3 (m+n+2)^3}P_3$$ with $$P_3=6 n^4+6(2 m+7) n^3+3 \left(2 m^2+19 m+36\right) n^2+$$ $$3 \left(5 m^2+29 m+40\right) n+\left(10 m^2+44 m+49\right)$$

Trying for $m=123$ and $n=456$ $$A_2=\color{red}{0.00046792}62719924$$ $$A_3=\color{red}{0.0004679271684}101$$

Without looking for such kind of expression, let $a=(n+1)$ and $b=(n+m+2)$ and use $$\psi ^{(1)}(a)-\psi ^{(1)}(b)=\left(\frac{1}{a}-\frac{1}{b}\right)+\frac 12\left(\frac{1}{a^2}-\frac{1}{b^2}\right)+\frac 16\left(\frac{1}{a^3}-\frac{1}{b^3}\right)+\cdots$$

Edit

If you want to make it more compact, switching from Taylor series to Padé approximant, you could use $$\psi ^{(1)}(p)\sim \frac{3 p (30 p+1)}{2 \left(45 p^3-21 p^2+3 p+2\right)}$$ whose error is $\sim\frac{1}{45 p^6}$.

For the worked case, the error would be $1.93\times 10^{-18}$.

Even better $$\psi ^{(1)}(p)\sim \frac{5 p \left(1302 p^2-573 p+697\right)}{6 \left(1085 p^4-1020 p^3+910 p^2-285 p+27\right)}$$ whose error is $\sim \frac{1}{74 p^8}$.