## Source: chap6.txt in ## http://www.stat.colostate.edu/computationalstatistics/datasets/computationalstatistics.zip # SIR: THE SAMPLING IMPORTANCE SAMPLING ALGORITHM # Example 6.3 on page 165 in Givens and Hoeting, 2nd edition # # Key idea: apply SIR to approximately simulate realizations from a # target distribution. Compare the performance when SIR is used to # simulate realizations from the normal distribution using the slash # distribution given on page 166 and vice versa. This is for # illustrative purposes only - we can generate samples from both of # these distributions easily. dslash=function(x) { # The slash density given on page 166 # input: x=quantile of the distribution # output: the height of the density at x ifelse(x!=0,(1-exp(-x^2/2))/x^2,1/2)/sqrt(2*pi) } #________________________________________________________________ # SIR to simulate draws from the normal distribution using draws from # the slash distribution. Y has a slash distribution if Y=X/U where # X~N(0,1) and U~Unif(0,1) # Step 1: sample candidates from the slash distribution set.seed(3) #set random seed x=rnorm(100000) u=runif(100000) slash.draws=x/u hist(slash.draws,nclass=40) stem(slash.draws) hist(slash.draws,xlim=c(-7,7), nclass=400000, freq=F) curve(dslash, -7, 7, add=T) # Step 2: Calculate the standardized weights in equation (6.12) weights=dnorm(slash.draws)/dslash(slash.draws) std.weights=weights/sum(weights) # Step 3: Resample the candidates with probability equal to the # standardize weights. # Note: you might get a warning message that begins "Walker's alias # method" when you run this. You can ignore it. sir.normal=sample(slash.draws,5000,prob=std.weights,replace=T) # Plot the results (similar to Figure 6.7 on page 166) par(mfrow=c(2,1)) hist(sir.normal,prob=T,yaxs="i",xlab="x",ylab="normal density") # for more options of "yaxs", try "?par" lines(sort(sir.normal),dnorm(sort(sir.normal))) #________________________________________________________________ # SIR to simulate draws from the slash distribution using draws from # the normal distribution # Step 1: sample candidates from normal distribution set.seed(3) #set random seed x=rnorm(100000) # Step 2: Calculate the standardized weights in equation (6.12) weights=dslash(x)/dnorm(x) std.weights=weights/sum(weights) # Step 3: Resample the candidates with probability equal to the # standardize weights sir.slash=sample(x,5000,prob=std.weights,replace=T) # Plot the results (similar to Figure 6.7 on page 166) hist(sir.slash,prob=T,xlim=c(-7,7),yaxs="i",,ylab="slash density") z=seq(-7,7,len=300) lines(z,dslash(z))