5
$\begingroup$

I'm trying to understand the SIR algorithm. In order to do so it seemed best to look at some existing r code (see below), taken from http://hedibert.org/wp-content/uploads/2013/12/example-iii.R.txt.

# SIR: Uniform proposal
N      = 10000
m      = 2000
draw1  = runif(N,-10,10)
w      = dnorm(draw1)/0.05
ind    = sample(1:N,size=m,replace=TRUE,prob=w)
draws1 = draw1[ind]
length(unique(draws1))/m

I don't really understand any of the steps taken after computing w. Why not just take the mean of w to approximate the integral? Why resample? Also, what is the last step about: Why calculate how many unique draws you have?

$\endgroup$
  • $\begingroup$ This particular proposal distribution doesn't even dominate the target. Technically, this example is incorrect. $\endgroup$ – Taylor Dec 8 '17 at 22:52
9
$\begingroup$

SIR uses two ideas. The first idea is importance sampling. The main idea is that you draw from one probability distribution (in your case, it's the uniform), in order to get information about another. You do this by drawing from one distribution, then weighting the samples. Generally you are trying to get information about a tricky distribution, but since this is a tutorial, you're trying to get information about a normal distribution with samples from a uniform.

Say you sample $X_1, \ldots, X_N \overset{iid}{\sim} \text{Uniform}(-10,10)$. In your code these are called draw1. Then you can weight these with the unnormalized weights $w_i = p_{\text{normal}}(X_i)/p_{\text{uniform}}(X_i)$. You can use these particles/weighted-samples to approximate things like expectations:

$$ E_{\text{normal}}[h(X)] \approx \sum_{i=1}^N \tilde{w}_i h(x_i), $$ where $\tilde{w}_i = w_i/\sum_j w_j$ are the normalized weights.

However these samples (without the weights) are not distributed according to the normal distribution. If you want samples that are (asymptotically) normally distributed, you need to resample from your weighted samples. Samples with higher weights are more likely to be picked. But at the end, all resampled things will have equal weight, as you are sampling with replacement.

So say you draw indexes $$ I_1, \ldots, I_m \overset{iid}{\sim} \text{Multinomial}(1, \tilde{w}_1, \ldots, \tilde{w}_N). $$ Your code calls these ind. The samples they select are distributed approximately normally. You can verify this empirically with a command like hist(draws1); it will look like a bell curve. Mathematically you would write these samples as $$ X_{I_1}, X_{I_2}, \ldots, X_{I_N} $$ (instead of $X_1, \ldots, X_N$.)

Now, you may have drawn duplicates here. Even though you might have two $3.6$s, they all are treated equal. Every resampled bit has equal weight. It should be $1/m$ for each of them.

Lastly, as I mentioned in a comment above, this example is technically incorrect. One of the requirements of a proposal/instrumental/importance distribution is that it should be able to draw samples evwrywhere in the support of your target distribution. In other words, the proposal should “dominate” your target. This proposal does not satisfy this criterion, though, because even if you choose a uniform that covers most of the support of the normal target (perhaps it is centered at he normal mode and is very wide), it still will leave some tail area uncovered. It might be true that this is negligible in practice, or that it is not a requirement for approximating certain specific expectations, but it’s worth mentioning.

$\endgroup$
  • $\begingroup$ So why would i want to have samples that are distributed according to the normal distribution? I wouldn't need this to approximate the integral, right? $\endgroup$ – gusdadjdk123 Aug 17 '16 at 13:25
  • 1
    $\begingroup$ @gusdadjdk123 not really this is just a demo. In real life, if you want samples you'd probably use rnorm and if you wanted integrals you'd do it by hand or look it up on wiki. But this would come in handy for distributions that might not have a name or nice properties/wiki pages. Or you could use it to verify pencil and paper work, maybe. $\endgroup$ – Taylor Aug 17 '16 at 13:31
  • $\begingroup$ Why do we divide dnorm(draw1) by 0.05? $\endgroup$ – MHall Feb 19 '20 at 2:08

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.