Mathematical clarification and Python implementation
In a earlier publish I launched completely different numerical sampling strategies, certainly one of them being significance sampling. In that publish we used this method to permit sampling from complicated distributions, sampling from which might in any other case be infeasible. Nonetheless, significance sampling is steadily used for one more cause, particularly variance discount: that’s, by selecting a suited proposal distribution we are able to scale back the variance of our estimator — which we are going to cowl right here.
Assume we don’t simply need to calculate the expectation E[X] of a random variable X, however as an alternative the expectation of a operate of that variable, f[X]. In a steady setting that is calculated as:
We are able to approximate this expectation utilizing numerical approximation, often known as Monte Carlo strategies, by sampling n random values from the distribution p after which calculate the pattern imply:
The thought behind significance sampling now could be to make use of a easy re-formulation trick and write the expectation as
— giving the expectation of f(x)p(x)/q(x) over the distribution q! And with that, permitting us to calculate the pattern imply by sampling from q:
The variance of the usual Monte Carlo estimator is given by:
The variance for the reformulated significance sampling estimator is:
In order a primary step we positively observe a distinction in variance, that means with excessive likelihood we are able to additionally discover a technique to scale back this. And certainly it’s comparatively simple to see that this variance could be lowered to 0 by selecting q as:
(Insert this time period within the equation above, and movie f(x)p(x) cancelling one another out — leaving Var[E[f(X)]]=0.)
Naturally, we don’t know E[f(X)], as the explanation we’re doing this sampling in spite of everything is to search out the expectation of f.
Nonetheless, we are able to consider E[f(X)] as some normalisation fixed, and no less than take one important perception away from this: we must always assemble q s.t. it has excessive density wherever f(x)p(x) is excessive. And with this, let’s dive right into a sensible instance and apply this studying.
For the sake of demonstration, we would like a sharp operate f, and a likelihood distribution p which don’t overlap too properly. Thus for simplicity allow us to set each to be regular distributions, e.g. f = N(5, 1) and p = N(9, 2):
I hope selecting a traditional distribution for each doesn’t confuse the reader, so let’s re-iterate what we’re attempting to do right here: we need to compute E[f(X)], the place X is a random variable which follows the distribution p — i.e. we need to compute the imply of f below p. Be aware this imply isn’t the imply normally related to a traditional distribution (which is a price on the x-axis, particularly the mode of the distribution), however now we’re after the imply of the y-values below p: on this instance it’s ~0.36 — a a lot lesser identified and used worth.
To approximate this numerically, as acknowledged above we might now pattern values x from the distribution p, and compute the empirical imply of f(x).
Intuitively one can see why sampling from this distribution is a foul thought, hopefully amplified by the earlier part: for many values sampled from p, f will probably be near 0 — however for just a few sampled values f will probably be very giant — thus we acquire a big variance.
Subsequently, following above launched concepts we now suggest a brand new distribution q = N(5.8, 1), which satisfies the derived criterion that its density is excessive in areas the place f(x)p(x) is excessive:
Be aware it’s not trivial to search out this operate, and positively there are rather more troublesome real-word situations. We’ve to attempt to fulfill the criterion in addition to attainable, but in addition care for satisfying the significance sampling requirement of p masking q, and so on. For this instance I really plotted p(x)f(x) after which picked a q which resembled it finest.
Python Implementation
Let’s code this in Python. First, we introduce the required features and distributions, and for comfort use functools.partials to acquire a operate representing a traditional distribution with fastened imply / customary deviation:
MEAN_F, STD_F = 5, 1
MEAN_P, STD_P = 9, 2
MEAN_Q, STD_Q = 5.8, 1def normal_dist(
imply: float, standard_deviation: float, x: np.ndarray
) -> np.ndarray:
return (
1
/ (standard_deviation * np.sqrt(2 * np.pi))
* np.exp(-0.5 * ((x - imply) / standard_deviation) ** 2)
)
f = partial(normal_dist, MEAN_F, STD_F)
p = partial(normal_dist, MEAN_P, STD_P)
q = partial(normal_dist, MEAN_Q, STD_Q)
Then, we generate the plot from above for orientation:
x = np.linspace(0, 15, 100)
plt.plot(x, f(x), "b-", label="f")
plt.plot(x, p(x), "r-", label="p")
plt.plot(x, q(x), "y-", label="q")
plt.legend()
plt.present()
Lastly, we come to the (significance) sampling half. First, we compute the direct Monte Carlo Estimator for E[f(X)]. We generate random samples x from p, and calculate the imply of f(x):
x_p = np.random.regular(loc=MEAN_P, scale=STD_P, measurement=NUM_SAMPLES)
y_p = f(x_p)
Now we apply significance sampling, i.e. pattern from qand proper by way of the significance weights:
x_q = np.random.regular(loc=MEAN_Q, scale=STD_Q, measurement=NUM_SAMPLES)
y_q = f(x_q) * p(x_q) / q(x_q)
Placing all of it collectively:
from functools import partialimport matplotlib.pyplot as plt
import numpy as np
NUM_SAMPLES = 1000000
MEAN_F, STD_F = 5, 1
MEAN_P, STD_P = 9, 2
MEAN_Q, STD_Q = 5.8, 1
def normal_dist(
imply: float, standard_deviation: float, x: np.ndarray
) -> np.ndarray:
return (
1
/ (standard_deviation * np.sqrt(2 * np.pi))
* np.exp(-0.5 * ((x - imply) / standard_deviation) ** 2)
)
f = partial(normal_dist, MEAN_F, STD_F)
p = partial(normal_dist, MEAN_P, STD_P)
q = partial(normal_dist, MEAN_Q, STD_Q)
x = np.linspace(0, 15, 100)
plt.plot(x, f(x), "b-", label="f")
plt.plot(x, p(x), "r-", label="p")
plt.plot(x, q(x), "y-", label="q")
plt.legend()
plt.present()
x_p = np.random.regular(loc=MEAN_P, scale=STD_P, measurement=NUM_SAMPLES)
y_p = f(x_p)
x_q = np.random.regular(loc=MEAN_Q, scale=STD_Q, measurement=NUM_SAMPLES)
y_q = f(x_q) * p(x_q) / q(x_q)
print(
f"Authentic imply / variance: {np.imply(y_p):.6f} / {np.var(y_p):.6f}"
)
print(
f"Significance sampling imply / variance: {np.imply(y_q):.6f} / {np.var(y_q):.6f}"
)
The output will probably be one thing like:
Authentic imply / variance: 0.036139 / 0.007696
Significance sampling imply / variance: 0.036015 / 0.000027
Thus, we nonetheless acquire the right imply, however have lowered variance by ~100x!
Significance sampling is a intelligent reformulation trick, permitting us to compute expectations and different moments by sampling from a distinct proposal distribution. This not solely permits sampling from complicated, in any other case hard-to-sample distributions, but in addition modifications the variance of the ensuing estimator. On this publish we confirmed the best way to make use of this to scale back the variance. Specifically, we proved and confirmed that deciding on a proposal distribution with excessive likelihood in areas the place p(x)f(x) (product of authentic distribution and performance in query) is excessive, yields finest outcomes.
Thanks for studying!
