What is the variance of a Polya Gamma distribution? I have a simple application that needs the variance of a Polya Gamma distribution (I know the mean since I found it here- http://arxiv.org/abs/1205.0310). 
This paper says that there is a closed form solution for the variance, but I am not really a mathematician and cannot calculate it.
In short, does anyone know the closed form solution for the variance of the Polya Gamma distribution.
I would be very grateful.
 A: In equation (6), the paper obtains the Laplace transform of these distributions as
$$\phi(t) = \prod_{k=1}^\infty \left(1 + \frac{t}{d_k}
\right)^{-b};\ d_k = 2\left(k-\frac{1}{2}\right)^2 \pi^2 + c^2/2$$
where $b\gt 0$ and $c\in \mathbb R$ are the parameters.  Taking logarithms yields 
$$\psi(t) = \frac{d}{dt}\phi(t) = \sum_{k=1}^\infty -b \log\left(1 + \frac{t}{d_k}\right).$$
This is a cumulant generating function (for imaginary values of $t$, at any rate) whose Taylor series around $t=0$ begins
$$\psi(t) = -\mu_1^\prime t + \frac{1}{2!} \left( \mu_2^\prime - \mu_1^{\prime \,2} \right)t^2 + \cdots$$
with $\mu_j^\prime$ representing the raw moment of order $j$: thus, the negative of the coefficient of $t$ is the mean and twice ($=2!$) the coefficient of $t^2$ is the variance.  The summation formula for $\psi$ can be expanded term-by-term and collected in common powers of $t$ to produce
$$\psi(t) = \sum_{k=1}^\infty -b \left(\frac{t}{d_k} + \frac{t^2}{2 d_k^2} + \cdots\right) = -b\sum_{k=1}^\infty \frac{1}{d_k} t - b\sum_{k=1}^\infty \frac{1}{2d_k^2} t^2 + \cdots.$$
Such sums, whose terms are the reciprocals of quadratic (and higher) functions of the integral index $k$, are straightforward to evaluate using the Weierstrass Factorization Theorem and yield
$$\mu_1^\prime = \frac{b }{2 c}\tanh \left(c/2\right);\ \mu_2^\prime - \mu_1^{\prime\,2} = \frac{b }{4 c^3}(\sinh (c) - c) \text{sech}^2\left(c/2\right).$$
The former agrees with the mean reported in the paper (adding some confidence to the overall correctness of this approach) while the latter answers the question: it is a closed form expression for the variance.  (These series can be continued in higher powers of $t$ to develop closed formulas for any cumulants, from which higher moments can be extracted.)
Partial contour plots show that the signs of the results (at least) are correct.

Nothing is shown along the $b$ axis (where $c=0$) because these formulas are not defined for $c=0$.  However, the plots make it clear that the formulas can be extended to continuous functions along that axis by taking the limits as $c\to 0$.  They give $b/4$ for the mean and $b/24$ for the variance.
