Mean value using importance sampling in Monte Carlo simulations

I'm studying Monte Carlo simulations from Mark Newman's book Monte Carlo Methods in Statistical Physics, and I've arrived at the concept of importance sampling. I'm having some trouble understanding the formula of the best estimate for a physical quantity $$Q$$ using importance sampling. The problem is the following:

The expectation value $$\langle Q\rangle$$ of a physical quantity $$Q$$ of a system in thermal equilibrium at temperature $$T$$ is $$\langle Q\rangle = \frac{\sum_{\mu} Q_{\mu} e^{-\beta E_{\mu}}}{\sum_{\mu}e^{-\beta E_{\mu}}},$$ where the sum runs over all states $$\mu$$ of the system. Importance sampling will choose at random a subset $$\{\mu_1, \cdots, \mu_N\}$$ from all states following a given probability distribution $$p_{\mu}$$. The best estimate for $$\langle Q\rangle$$ using this subset is said to be (without further justification)

$$\langle Q\rangle = \frac{\sum_{i=1}^{N} Q_{\mu_i} p_{\mu_i}^{-1} e^{-\beta E_{\mu_i}}}{\sum_{j=1}^{N} p_{\mu_j}^{-1} e^{-\beta E_{\mu_j}}}.$$

I can't understand the last formula. Why is it the best estimate? Where does the $$p_{\mu}^{-1}$$ come from?

Say we want to compute the mean value of $$X$$ given a probability distribution $$P(X)$$. Denote this quantity $$\langle X \rangle_P$$. Then, by definition: $$\begin{eqnarray} \langle X \rangle_P &=& \sum_X X P(X) \end{eqnarray}$$ Suppose we draw $$M$$ samples $$\hat{X}_j$$ ($$i=1,2,...M$$) from $$P$$. Then we can compute an estimate of $$\langle X \rangle_P$$ using an average over the samples... $$\begin{eqnarray} \langle X \rangle_P &=& \frac{1}{M} \sum_{j=1}^M \hat{X}_j \end{eqnarray}$$ Easy enough!

Now suppose we have another probability distribution $$Q(X)$$. We can express $$\langle X \rangle_P$$ in terms of an average $$\langle \cdots \rangle_Q$$ over the distribution $$Q$$ using the following trickery $$\begin{eqnarray} \langle X \rangle_P &=& \sum_X X P(X) \frac{Q(X)}{Q(X)} \\ &=& \sum_X \left[X\frac{P(X)}{Q(X)}\right] Q(X) \\ &=& \langle X \frac{P}{Q}\rangle_Q \end{eqnarray}$$ The basic idea of importance sampling takes advantage of this identity by saying that to compute $$\langle X \rangle_P$$, it may be more convenient to first draw samples from $$Q$$, and then correct the naive average you would obtain using the samples from $$Q$$ by using a weighted average with weights $$P/Q$$. A typical reason to do this would be if $$P$$ has a tail which is difficult to sample from directly; you can instead sample from $$Q$$ which has a more uniform distribution, and apply weights to get the average with respect to $$P$$.

In more concrete terms, suppose we draw samples $$\hat{X}_i$$ ($$i=1,2,...N$$) from $$Q$$. [Very important note: even though I'm using similar notation, $$\hat{X}_i$$, $$i=1,...,N$$, for samples drawn from $$Q$$ that I used for samples drawn from $$P$$ above, $$\hat{X}_j$$, $$i=1,...,M$$, it's crucially important to understand that these two sets of samples are drawn from different distributions have different statistical properties. From here on out in the answer, I will only refer to samples drawn from $$Q$$, denoted by $$\hat{X}_i$$] Then we can estimate $$\langle X \rangle_P$$ as $$\begin{eqnarray} \langle X \rangle_P &=& \frac{1}{N} \sum_{i=1}^N \hat{X}_i \frac{P(\hat{X}_i)}{Q(\hat{X}_i)} \end{eqnarray}$$

Now let's turn to your example. The abstract definitions of $$P$$ and $$Q$$ I gave are related to the distributions in your question via $$\begin{eqnarray} P &=& \frac{e^{-\beta E_\mu}}{\sum_\mu e^{-\beta E_\mu}} \\ Q &=& p_\mu \end{eqnarray}$$ Note: I relabeled the quantity $$Q$$ in your question to $$X$$ to avoid conflicting with the distribution $$Q(X)$$ I am using in my answer.

We can estimate $$\langle X \rangle_P$$ using importance sampling via $$\begin{eqnarray} \langle X \rangle_P &=& \langle X \frac{P}{Q} \rangle _Q \\ &=& \langle X \cdot \frac{e^{-\beta E_\mu}}{\left[\sum_\mu e^{-\beta E_\mu}\right]}\cdot \frac{1}{p_\mu} \rangle _Q \end{eqnarray}$$ If we draw samples $$\hat{X}_i$$ ($$i=1,2,...,N$$) from $$Q$$, then our estimate of $$\langle X \rangle_P$$ is $$\begin{eqnarray} \langle X \rangle_P &=& \frac{1}{N} \frac{ \sum_{i=1}^N \hat{X}_i e^{-\beta E_{\mu_i}} p_{\mu_i}^{-1} }{\sum_{\mu} e^{-\beta E_{\mu}}} \end{eqnarray}$$ This is almost what you want, except we have an overall factor of $$1/N$$ and the sum in the denominator is over states $$\mu$$ and not over samples labeled by $$i$$.

No matter: the trick is that we can also estimate the partition function (the sum in the denominator over $$\mu$$) using importance sampling. On the one hand, of course $$\langle 1 \rangle_P=1$$ (the average of $$1$$ is $$1$$). On the other hand, we have that $$\begin{eqnarray} \langle 1 \rangle_P &=& \langle \frac{P}{Q}\rangle_Q \end{eqnarray}$$ Or, in terms of the samples: $$\begin{eqnarray} \frac{1}{N} \sum_{i=1}^N \frac{P(\hat{X}_i)}{Q(\hat{X}_i)} = 1 \end{eqnarray}$$ Plugging in the explicit expressions for your distributions: $$\begin{eqnarray} \frac{1}{N} \frac{\sum_{i=1}^N e^{-\beta E_{\mu_i}} p_{\mu_i}^{-1}}{\sum_\mu e^{-\beta E_\mu}} = 1 \end{eqnarray}$$ or $$\begin{eqnarray} N \sum_\mu e^{-\beta E_\mu} = \sum_{i=1}^N e^{-\beta E_{\mu_i}} p_{\mu_i}^{-1} \end{eqnarray}$$ Using this trick, we have our final expression for the estimate of $$\langle X \rangle_P$$ written entirely in terms of sums over samples: $$\begin{eqnarray} \langle X \rangle_P &=& \frac{ \sum_{i=1}^N \hat{X}_i e^{-\beta E_{\mu_i}} p_{\mu_i}^{-1} }{\sum_{i=1}^N e^{-\beta E_{\mu_i}} p_{\mu_i}^{-1}} \end{eqnarray}$$ which is the expression you wanted.