1
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ 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). $\endgroup$ – JimB Oct 31 '20 at 5:11
  • $\begingroup$ @JimB, I think you are right on the density. $\endgroup$ – Leslie Chiu Oct 31 '20 at 16:42
3
$\begingroup$

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[12345];
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[12345];
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 *)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.