random hyperharmonic series The Harmonic Series is defined as:
$\sum_{n} \frac{1}{n}$ where $n=1,2,3,4....$. 
This series is known to be divergent.
A generalization of this series can be made by raising each term to $p$:
$\sum_{n} \frac{1}{n^p}$  which is also known as the hyperharmonic series and is known to be convergent when  $p>1$.
On the other hand, for $p=1$, if the signs of the terms are alternating the sum:
$\sum_{n} \frac{(-1)^n}{n}$ is convergent and approaches $\ln{2}$.
A natural extension would be to introduce randomness in the sign of each terms. 
$\sum_{n} \frac{\epsilon_{n}}{n}$ where $\epsilon_{n}$ is defined by the probability of its outcome: $P(\epsilon_{j} = 1)=P(\epsilon_{j} = -1)=1/2$. This is called the random harmonic series in Schmuland (http://www.stat.ualberta.ca/people/schmu/preprints/rhs.pdf).
My question is: What would happen if we generalize this to the case of the random HYPERharmonic series? What would be the distribution of the result of the summation?
 A: The case $p=2$ is treated briefly in the final section of the cited
Schmuland paper, which gives a picture of the distribution.  The observations
in the paper's first few sections adapt to arbitrary $p>1$: the distribution
is not expected to have a simple formula, but the moment-generating function
has a product formula (exhibited below),
which we can use to compute each power moment ${\mathbb E}(X_p^{2k})$
of this random variable $X_p = \sum_n \epsilon_n/n^p$
as a polynomial in $\zeta(2p), \zeta(4p), \zeta(6p), \ldots, \zeta(2kp)$.
[The odd-order moments vanish by symmetry.]
This is because the distribution of $X_p$ is the convolution of
an infinite series of distributions, the $n$-th of which is supported on
$\pm 1 /n^p$ each with probability $1/2$; therefore
$$
{\mathbb E}(\exp(tX_p))
= \prod_{n=1}^\infty \left(\frac12 e^{t/n^p} + \frac12 e^{-t/n^p} \right)
= \prod_{n=1}^\infty \cosh(t/n^p).
$$
We can recover the power moments from the expansion of
${\mathbb E}(\exp(tX_p))$ in a Taylor series about $t=0$,
writing
$$
\begin{eqnarray}
\log {\mathbb E}(\exp(tX_p))
&=& \sum_{n=1}^\infty \log(\cosh(t/n^p)) \cr
&=& \sum_{n=1}^\infty
\frac12 \left(\frac{t}{n^p_{\phantom1}}\right)^2
- \frac1{12} \left(\frac{t}{n^p_{\phantom1}}\right)^4
+ \frac1{45} \left(\frac{t}{n^p_{\phantom1}}\right)^6
- \frac{17}{2520} \left(\frac{t}{n^p_{\phantom1}}\right)^8
+ - \cdots \cr
&=&
\frac{\zeta(2p)}{2} t^2
- \frac{\zeta(4p)}{12} t^4
+ \frac{\zeta(6p)}{45} t^6
- \frac{17\zeta(6p)}{2520} t^8
+ - \cdots
\end{eqnarray}
$$
and exponentiating.
While the random variables $X_p$ may be no more than objects of curiosity,
similar constructions arise naturally in analytic number theory, such as
the value at a fixed $s>1$ of the Dirichlet $L$-series associated to a random
real character $\chi$.  In this example, the terms $\chi(n)/n^s$ in the
Dirichlet series are correlated, but the Euler product
$L(s,\chi) = \prod_l (1 - \chi(l)/l^s)^{-1}$ has independent factors, so
$\log L(s,\chi)$ is an infinite convolution of the same kind.  Likewise
if the $\epsilon_n$ were random complex numbers of unit length:
the analogue could be $\log L(s,\chi)$ for a random Dirichlet charater $\chi$
that need not be real, or $\log \zeta(\sigma + it)$ for fixed $\sigma>1$
and random real $t$.  For these sums of complex-valued random variables,
the factors in the moment generating function get more complicated than
hyperbolic cosines, but are still tractable.
