Closed form function relating $\mu$ to the natural parameter for the logarithmic series distribution?

Question

While answering another question here, I mentioned the logarithmic series distribution as a possible model for species per genus.

In the course of looking at the pmf while answering that I realized that it was exponential family (a fact I hadn't previously been aware of).

The parameterization of the pmf at the Wikipedia page (which up to change of symbols is the only form I've seen it in) is of the form

$$f(y) = c(p) . p^y/y\,,\qquad 0<p<1\,;\: y=1,2,3,...$$

where $c(p)=\frac{1}{\log(\frac{1}{1-p})}=-\frac{1}{\log(1-p)}$ is the normalizing constant.

By examining it in "standard form" for the exponential family we see the natural parameter, $\theta=\log(p)$ (and also that the sufficient statistic $T(x)=x$, so the MLE for a constant $p$ should be a function of the sample mean, and in that case we'd expect MLE should be equivalent to method of moments, which turns out to be the case).

[So it looks like if we had a glm with a model $\theta=\log p=X\beta$, then $X'y$ should be sufficient for $\beta$.]

It's easy enough to derive $\mu$ as a function of $p$ (the series is trivial to sum), or (looking at it in terms of the exponential family) to take $b'(\theta)$, and either way write $\mu$ in terms of $\theta$:

$$\mu = c(p)\cdot \frac{p}{1-p} = -\frac{1}{\log(1-e^\theta)}\cdot \frac{e^\theta}{1-e^\theta}\,.$$

It's even quite feasible to compute the variance function $V(\mu)$.

However it's not immediately clear to me how to invert this mean function to write $\theta$ (or equivalently, $p$) as a function of $\mu$ (which would be convenient to have, since it's the link function). It may be tiredness or it might actually be difficult, but I couldn't see how to get anywhere with it... so to the question:

Can $\theta(\mu)$ be written in closed form?

($p(\mu)$ will do just as well)

I don't wish to be overly pedantic about what's included in closed form; it's useful to be able to write these things in terms of well-known functions with known properties, which may be relatively easily evaluated by calling some function. Expressions in terms of more-or-less "standard" functions that are commonly available in packages, like Gamma functions or even say Lambert-W functions (or even the various forms of hypergeometric functions I guess) could count -- I'm more interested in what useful things we can say about $\theta(\mu)$ and how easily we can evaluate it (beyond the obvious step of taking the above equation for $\mu$-as-a-function-of-$\theta$ and just solving for $\theta$ as needed).

Answer 1

You have the answer in your question: Lambert-W.

In R you can use

# Load library
  library(LambertW)

# Define function to obtain p from mu
  p <- function(mu) {1 - exp((1/mu) + W(-1/(mu*exp(1/mu)),-1))}

# Show (visually) that the function provides the correct values
  mu <- 1 + 10*c(0:100)/100
  plot(p(mu), mu, xlab="p", ylab="mu", las=1)
  pp <- c(1:99)/100
  lines(pp, -pp/((1-pp)*log(1-pp)))

to obtain

In Mathematica you can use

data = Table[{1 - Exp[(1/mu) + 
       ProductLog[-1, -1/(mu Exp[1/mu])]], mu}, {mu, 1.01, 5, 0.01}];
cp = -1/Log[1 - p];
ListPlot[{data, Table[{p, cp p/(1 - p)}, {p, 0.1, 0.93, 0.01}]}, 
 PlotStyle -> {{PointSize[0.03]}, {PointSize[0.01]}}]  

(Mathematica was used to determine the general form of the relationship for determining $p$ from $\mu$.)