The variance function of the classical occupancy distribution is given in various references$^\dagger$ as:
$$V(n,m) = m \Bigg[ \Big( 1 - \frac{1}{m} \Big)^n + (m-1) \Big( 1-\frac{2}{m} \Big)^n - m \Big( 1-\frac{1}{m} \Big)^{2n} \Bigg] \quad \quad \quad \text{for all } n,m \in \mathbb{N}.$$
This is a variance so it must be non-negative over all argument values. However, when I compute this function in R I get some negative values for low values of $n$ and large values of $m$. I suspect this is an underflow problem, but I first want to verify that the variance function is indeed correct (i.e., check that it is a non-negative function).
Question 1: Show that this function is non-negative over all argument values.
Question 2: Assuming that this function is indeed non-negative, what is the best way to compute this in R so that I avoid negative values due to underflow?
$^\dagger$ See e.g., Johnson, N.L. and Kotz, S. (1969) Discrete Distributions. John Wiley and Sons: New York, p. 252.
For any fixed $0<n<1$, and large $m$ you will actually have $V(n,m)<0$; this is not an issue with your computational software. Indeed, for any fixed $n>0$, by computing the Taylor series for $V(n,m)$ at $m=\infty$, we see that as $m\to\infty$, we have $V(n,m)={n(n-1)\over m}+O(m^{-2})$. Hence, for any fixed $0<n<1$ and $m$ large enough compared to $n$, we have $V(n,m)<0$.
Unfortunately, I don't see how to prove that $V(n,m)\geq 0$ for $n,m\in\mathbb{N}$ at the moment.