531f10MCIb

531f10MCIb - h=function(x){ exp(-x)*sqrt(x)/gamma(3/2)} X =...

Info iconThis preview shows pages 1–2. Sign up to view the full content.

View Full Document Right Arrow Icon
############################################################# # STAT 531 MCIb ############################################################# set.seed(111) h=function(x){ 1/(x^2*sqrt(2*pi)*exp(1/(2*x^2)))} par(mfrow=c(2,1)) curve(h,from=0,to=1/4.5,xlab="x",ylab="h(x)",lwd="2") m=10^4 I= (1/4.5)*h(runif(m)/4.5) estint= cumsum(I)/(1:m) esterr= sqrt(cumsum((I-estint)^2))/(1:m) estvar= cumsum((I-estint)^2)/(1:m)^2 plot(estint,xlab="Iterations",ty="l",lwd=2, ylim=mean(I)+20*c(- esterr[m],esterr[m]),ylab="") lines(estint+2*esterr,col="gold",lwd=2) lines(estint-2*esterr,col="gold",lwd=2) estint[m] esterr[m] estvar[m] integrate(h,0,1/4.5) pnorm(-4.5) ############################################################# x =rnorm(10^4) bound= -4.5 mean(x[1:10^4] < bound ) ############################################################# set.seed(111) x=rexp(10^4) mean(dnorm(x+4.5)*exp(x)) #############################################################
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: h=function(x){ exp(-x)*sqrt(x)/gamma(3/2)} X = rexp(10^4,1) + 12.5 I=exp(-12.5)*sqrt(X)/gamma(3/2) estint=cumsum(I)/(1:10^4) esterr=sqrt(cumsum((I-estint)^2))/(1:10^4) estint[10000] esterr[10000] plot(estint,xlab=&quot;Iterations&quot;,ty=&quot;l&quot;,lwd=2, ylim=mean(I)+20*c(-esterr[10^4],esterr[10^4]),ylab=&quot;&quot;) lines(estint+2*esterr,col=&quot;gold&quot;,lwd=2) lines(estint-2*esterr,col=&quot;gold&quot;,lwd=2) abline(a=pchisq(25,3, low=F),b=0,col=&quot;gold3&quot;,lwd=2) integrate(h,12.5,Inf) pchisq(25,3,low=F) ############################################################# set.seed(111) m=10^4 x=rexp(m) + 4.5 I=dnorm(x)/dexp(x-4.5) estint=cumsum(I)/(1:m) esterr=sqrt(cumsum((I-estint)^2))/(1:m) estvar=cumsum((I-estint)^2)/(1:m)^2 plot(est,type=&quot;l&quot;,xlab=&quot;Iterations&quot;,ylab=&quot;&quot;,lwd=2) abline(a=pnorm(-4.5),b=0,col=&quot;gold3&quot;,lwd=2) estint[m] esterr[m] pnorm(-4.5) #############################################################...
View Full Document

Page1 / 2

531f10MCIb - h=function(x){ exp(-x)*sqrt(x)/gamma(3/2)} X =...

This preview shows document pages 1 - 2. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online