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.

enter image description here

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).


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

In R you can use

# Load library

# 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

R Lambert W

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$.)

| cite | improve this answer | |
  • $\begingroup$ Thanks. After spending some time on it I eventually figured out how to manipulate it into a form that Wolfram alpha would solve before timing out (this was a few days back now), but I figured that I would let the bounty run so someone could get the benefit. Your answer is much nicer than anything I'd have posted, so it well deserves the bounty. $\endgroup$ – Glen_b Aug 16 '15 at 7:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.