HW 6 Solution ============= (a) > g = function(x,y) (3*x^2 + x*y + 5*y^3)/(1+x^2 + y^2)^3 > g.inner = function(y) { ny = length(y) out = err = numeric(ny) for(i in 1:ny) { tmp =integrate(g, -2,2, y=y[i]) out[i] = tmp$val err[i] = tmp$abs.error } attr(out,"max.error") = max(err) out } > curve(g.inner(x),-3,3) ### This works which shows that g.inner vectorizes > attr(g.inner(seq(-3,3,by=.005)),"max.error") [1] 0.0001583349 > integrate(g.inner,-3,3) 1.672408 with absolute error < 4.7e-05 ### But this estimates only the error due to outer integral. ### Overall error is bounded by 1.583e-4+4.73e-5 = 2.056e-4 gvals = g(runif(1e6,-2,2),runif(1e6,-3,3)) > cat("MC integral 99% CI=",4*6*mean(gvals)," + or - ", 24*qnorm(.995)*sd(gvals)/sqrt(1e6),"\n") MC integral 99% CI= 1.671374 + or - 0.01477804 (b) Find the breakpoints between analytically defined pieces first: f(x) = C* max( exp(-abs(2*x)), dnorm(x), dbeta((x+2)/4,2,1) ) Quadratic formula says first expression smaller than second iff |x| between > w = 2 + c(-1,1)*sqrt(4-2*log(sqrt(2*pi))) [1] 0.5295841 3.4704159 Also dnorm(x) < dbeta((x+2)/4,2,1) for 0 uniroot(function(x) dnorm(x)- dbeta((x+2)/4,2,1),c(-2,0)) $root [1] -1.85799 $f.root [1] 1.191693e-07 $iter [1] 4 $init.it [1] NA $estim.prec [1] 6.103516e-05 > y = .Last.value$root > y [1] -1.85799 > z = uniroot(function(x) dbeta((x+2)/4,2,1) - exp(2*abs(x)),c(-2,0))$root > z [1] -1.960318 > curve(exp(-2*abs(x)),-2,2) curve(dnorm(x),-2,2, add=T, col="red") curve( dbeta((x+2)/4,2,1), add=T, col="blue") ### From the plot, or from the comparisons indicated above, we can see ### that ### that the max density f defined above can be alternatively expressed as: f(x) = C*dnorm(x) for -2 < x < -1.85799, and = C*dbeta((x+2)/4 = C*(x+2)/2 for -1.85799 < x < 2 = 0 otherwise ## The constant is then const = 1/(pnorm(-1.85799)-pnorm(-2) + 4 - 0.25*(2-1.85799)^2) = 0.249763 ## Using the method explained in class, we then simulate as follows p1 = const*(pnorm(-1.85799)-pnorm(-2)) = 0.002206679 p2 = 1-p1 q1 = pnorm(-1.85799)-pnorm(-2) = 0.0088350 q2 = 4 - 0.25*(2-1.85799)^2 Lvec = rbinom(1e5, 1, p1) Xvec = numeric(1e5) Xvec[Lvec==1] = qnorm( pnorm(-2) + q1*runif(sum(Lvec==1)) ) Xvec[Lvec==0] = -2 + 2*sqrt(q2*runif(sum(Lvec==0)) + 0.25*(2-1.85799)^2) > f = const*pmax( exp(-abs(2*x)), dnorm(x), dbeta((x+2)/4,2,1) ) > hist(Xvec,nclass=50,prob=T) curve(f(x), add=T, col="green") #### OK !!! #------------------------------------------------------------------------ ### In my own notes, I originally stated the problem with f defined as ### the "min" of the three densities. In that case, the density would be ### proportional to: dbeta((x+2)/4,2,1) on (-2,z) exp(-2*abs(x)) on (z,-w[1]) U (w[1],2) dnorm(x) on (-w[1],w[1]) To find constant, note const = 1/(4*pbeta((z+2)/4,2,1)+0.5*(exp(-2*w[1])-exp(2*z))+ 0.5*(exp(-2*w[1])-exp(-4)) + pnorm(w[1])-pnorm(-w[1])) > const [1] 1.366745 ## masses of the three subintervals are qvec = c(4*pbeta((z+2)/4,2,1), 0.5*(exp(-2*w[1])-exp(2*z)), pnorm(w[1])-pnorm(-w[1]), 0.5*(exp(-2*w[1])-exp(-4))) pvec = const*qvec > pvec [1] 0.0005380417 0.2234051484 0.5516178259 0.2244389839 > f = function(x) ifelse(x curve(f(x), add=T, col="green") #### OK, have the right density function ### Then simulation gives: > Lv = sample(1:4, 1e5, rep=T, prob=pvec) > table(Lv) Lv 1 2 3 4 53 22418 55262 22267 Xv = numeric(1e5) Xv[Lv==1] = qbeta(runif(sum(Lv==1))*qvec[1]/4,2,1)*4-2 Xv[Lv==2] = 0.5*log(2*qvec[2]*runif(sum(Lv==2))+exp(2*z)) Xv[Lv==3] = qnorm(runif(sum(Lv==3))*qvec[3]+pnorm(-w[1])) Xv[Lv==4] = -0.5*log(exp(-2*w[1])-2*qvec[4]*runif(sum(Lv==4))) > hist(Xv,nclass=50,prob=T) > curve(f(x)*const, add=T, col="green") #### OK !!!