# How to do importance sampling to estimate variance?

For a random variable X with density $$q(x)=e^{(1/3)∗|x|^3}$$, I can estimate its variance by numerical integration method.

Then, how can I do importance sampling using a standard normal proposal (compare the tails of q and the standard Normal distribution). That is, using the importance sampling method, I should sample $$Y ∼ g = N(0, 1)$$, get a new estimate for the variance of X (Based on a sample of size N = 1000). Thanks in advance.

• But $q(x)$ as written is not a probability density function. Do you mean $q(x)=\exp \left(-\frac{| x| ^3}{3}\right)$ (which is proportional to a pdf) or $q(x)=\frac{3^{2/3} \exp \left(-\frac{| x| ^3}{3}\right)}{2 \Gamma \left(\frac{1}{3}\right)}$ (which is a pdf).
– JimB
Oct 31 '20 at 5:11
• @JimB, I think you are right on the density. Oct 31 '20 at 16:42

Once $$q(x)$$ is normalized, then the mean of $$X$$ is 0. Therefore, the variance is the expectation of $$X^2$$. To estimate the variance using self-normalizing importance sampling one can perform the following steps:

(* Define the function proportional to the pdf of the nominal distribution *)
q[x_] := Exp[-Abs[x]^3/3]

(* pdf of proposal distribution (also called the importance distribution *)
p[x_] := PDF[NormalDistribution[], x]

(* Get a random sample from the proposal distribution *)
SeedRandom;
y = RandomVariate[NormalDistribution[], 1000];

(* Calculate the weights *)
w = q[#]/p[#] & /@ y;

(* Calculate the estimate of the variance *)
variance = (y^2).w/Total[w]
(* 0.767868 *)


If the proportionality constant for the pdf of the nominal distribution was known, then the following steps would obtain an estimate of the variance using importance sampling:

(* Define the function proportional to the pdf of the nominal distribution *)
q[x_] := (3^(2/3)/(2 Gamma[1/3])) Exp[-Abs[x]^3/3]

(* pdf of proposal distribution (also called the importance distribution *)
p[x_] := PDF[NormalDistribution[], x]

(* Get a random sample from the proposal distribution *)
SeedRandom;
y = RandomVariate[NormalDistribution[], 1000];

(* Calculate the weights *)
w = q[#]/p[#] & /@ y;

(* Calculate the estimate of the variance *)
variance = Mean[w y^2]
(* 0.768766 *)


The "true" variance is given by

Integrate[x^2 (3^(2/3)/(2 Gamma[1/3])) Exp[-Abs[x]^3/3], {x, -∞, ∞}]
(* 3^(2/3)/Gamma[1/3] *)
(* Approximately 0.776458 *)