Variance of the Random Variable $|im(f)|$ where $f:[n] \rightarrow [n]$ Here's a question:

Let $f$ be picked randomly from the set of all functions from $[n]$ to $[n]$, where $[n]$ is the set $\{1,2,3,\ldots,n\}$. Give a closed-form expression for the variance of the random variable $|im(f)|$, where bars denote the size of the set, and $im(f)$ is the image of $f$.

I know that for a random variable $X$, variance $V(X) = E(X^2) - (E(X))^2$, where $E(X)$ is the expected value of $X$. Beyond that, however, I'm unsure how to proceed. Can I get some help, please?
 A: Sorry, this is a little bit long, but I don't know how to write it shorter!
Suppose to fix the set $B \subseteq [n]$, with $|B| = m$ ($m$ is the cardinality of the set $B$). Then there are exactly $n^m$ function $f : [n] \rightarrow B$.
Now, we most count the number of subsets $B$ present in $[n]$ with fixed cardinality $|B| = m$. This number is the binomial coefficient of $n$ and $m$. Namely:
$$\mathcal{B}(n,m) = \frac{n!}{m!(n-m)!}$$
At this point, we can say that there are
$$N(n,m) = n^m\mathcal{B}(n,m)$$
functions $f : [n] \rightarrow B$, with $|B| = m$.
Remark 1. Recall the binomial formula $\displaystyle\sum_{m=0}^n a^{n-m}b^m \mathcal{B}(n,m) = (a + b)^n$
Now, we want to know the total number of function $f : [n] \rightarrow B$, where $B$ is a generic subset of $[n]$. We need to sum up all the terms $N(n,m)$. Then:
$$\sum_{m=1}^n n^m \mathcal{B}(n,m) = \sum_{m=0}^n n^m \mathcal{B}(n,m) - 1 = \sum_{m=0}^n 1^{n-m}n^m \mathcal{B}(n,m) - 1 = (1 + n)^n - 1$$
So we pose $N(n) = (1+n)^n - 1$ the total number of possible functions.
Remark 2. for each function $f$, we have that $im(f) = B$.
If we randomly pick a function $f$ among the all possible, then the probability that $|im(f)| = m$ is:
$$Pr(|im(f)| = m) = \frac{N(n,m)}{N(n)}$$ 
Now, can pass to evaluate $E(X)$:
$$E(X) = \sum_{m=1}^n m \frac{N(n,m)}{N(n)} = \frac{1}{N(n)}\sum_{m=1}^n m n^m B(n,m)$$
Do evaluate this, let's pose $g_n(a) = \sum_{m=1}^n a^m \mathcal{B}(n,m) = (1+a)^n - 1$ (binomial formula). If we derivate both terms with respect to $a$, we get: 
$$\frac{d g_n(a)}{da} = \sum_{m=1}^n m a^{m-1} \mathcal{B}(n,m) = n(1+a)^{n-1}$$
Substituting $a$ with $n$, and multiplying by $n$ we get:
$$\sum_{m=1}^n m n^{m} \mathcal{B}(n,m) = n^2(1+n)^{n-1}$$
You can notice that:
$$E(X) = \frac{1}{N(n)}\sum_{m=1}^n m n^{m} \mathcal{B}(n,m) = \frac{n^2(1+n)^{n-1}}{(1+n)^n - 1}$$ 
With similar calculation, we can derive $E[X^2]$. In fact, consider the second derivative of $g_n(a)$:
$$\frac{d^2 g_n(a)}{da^2} = \sum_{m=1}^n m(m-1) a^{m-2} \mathcal{B}(n,m) = n(n-1)(1+a)^{n-2}$$
Also, if we multiply all by $a^2$, we have:
$$a^2\frac{d^2 g_n(a)}{da^2} = \sum_{m=1}^n (m^2-m) a^{m} \mathcal{B}(n,m) = a^2n(n-1)(1+a)^{n-2}$$
In particular we can split the summation, substitute $n$ to $a$ and write the following:
$$\sum_{m=1}^n m^2 n^{m} \mathcal{B}(n,m) = n^3(n-1)(1+n)^{n-2} + \sum_{m=1}^n m n^{m} \mathcal{B}(n,m) = n^3(n-1)(1+n)^{n-2} + E[X]N(n)$$
Finally:
$$E(X^2) = \frac{1}{N(n)}\sum_{m=1}^n m^2 n^m \mathcal{B}(n,m) = \frac{n^3(n-1)(1+n)^{n-2}}{(1+n)^n-1} + E[X]$$ 
Now, you can evaluate the variance of $|im(f)|$.
