I answered 7 question with R and I need (good interpretation for questions 1 to 7), and answer for question 8.
questions 1-7: only interpretation to the graphs (some questions ask for specific interpretation )
I need it done after 6 hours maximum
SRM 620 Advanced Statistical Computing Fall 2018 Homework Assignment #1 Simulation of random variables. Note 1. Always use captions and numbers to identify each graph and table. And they should be referred to in the text. Remember, a figure or table should ideally be self explanatory. That means, if you take the figure out of the context, it should still be meaningful. 2. Do not stress yourself if you cannot change the R codes. You may discuss things with your friends. However, you should prepare the report by yourself. 3. You will be evaluated based on completeness, quality of interpretations, and quality of presentation. 4. Always think that you are writing it for a peer reviewed journal. Submission 1. Answer all the questions. 2. Use Word or any word-processor to typeset your answer. Use 12 point Roman font, 1.5 spaced. Print on one side of the paper. 3. If you are pasting R codes or outputs, make sure to change the font to Courier New (10 or 11 point as appropriate) for those portions of the report. 1 Problems Q.1 (15 points) Triangular distribution. Simulate sums of pairs of independent uniformly distributed random variables. Study the histogram. Use your insight to develop an algorithm for drawing from a triangular distribution. The triangle has it peak at zero and its base runs from 0 to 2. Q.2 (10 points) Sums of squares. Generate sums of squares of 5 independent normal variables. Compare their distribution to the appropriate χ2 distribution. Q.3 (10 points) Kernel density estimation. Generate exponentially distributed random numbers. Use density() to plot the empirical density. Compare to the histogram. Do you notice a problem? Q.4 (15 points) Rayleigh distribution. Simulate pairs of independent standard normal variables and interpret them as x and y of a vector. Compute the length of the vector and study its distribution. Q.5 (15 points) Box-Müller. Take a look at example 4 at the beginning of Section 3.4: “4. If U, V∼Unif(0,1) are independent, then Z1 = √ −2 logUcos(2πV ), Z2 = √ −2 logUsin(2πV ) are independent standard normal variables (see e.g. [238, p. 86]).” It describes the Box-Müller method for generating normally distributed numbers. Implement it in R. Q.6 (10 points) Cauchy distribution. Draw 1000 sets of numbers from the Cauchy distribution. Do this for set size 2, 5, 10 and 20. Compute the median of each set (be smart, use a matrix and apply()). Study the distribution of the medians, for each set size. Q.7 (10 points) Cauchy distribution again. Do the same as above, but do not compute medians of sets, but means. Do you notice anything special? Q.8 (15 points) Acceptance-rejection sampling. In a formula for a density usually a numerical constant occurs, to make it a proper density, i.e. the integral is 1. In 2 the standard normal density, exp(−x2/2)/ √ 2π, it is the factor 1/ √ 2π. Suppose I don’t know or forget this “normalizing constant”, how would that influence on acceptance-rejection sampling? 3 x1<- apply(data.frame(runif(100,0,1),runif(100,0,1)),1,sum)=""><-apply(data.frame(runif(10000,0,1),runif(10000,0,1)),1,sum) hist(x[,1],main="TeX("$n=100$")," xlab="TeX("Sapmple1"))" hist(x[,2],main="TeX("$n=10000$")," xlab="TeX("Sample2"))" xlo="0" xhi="2" x="seq(xlo," xhi,="" by="0.1)" ff=""><- function(x){(2-x)/2}="" f="ff(x)" plot(x,="" f,="" lwd="2," type='l' ,="" col='blue' )="" u="max(f)" n="1000" a="rep(0," n)="" y="runif(n)" *="" (xhi="" -="" xlo)="" +="" xlo="" z="runif(n)" *="" u="" a="1" *="" z="">< ff(y) g = c('-', '+') points(y, z, pch = g[a + 1]) y1=rnorm(10000, mean = 0, sd = 1)^2 y2=rnorm(10000, mean = 1, sd = 2)^2 y3=rnorm(10000, mean = 2, sd = 2)^2 y4=rnorm(10000, mean = 3, sd = 3)^2 y5=rnorm(10000, mean = 4, sd = 2)^2 sumsquared=cbind(y1,y2,y3,y4,y5) sumsquared=cbind(apply(sumsquared,1,sum)) sumsquared=apply(matrix(rnorm(10000, mean = 0, sd = 1)^2*5,10000,5),1,sum) chi=rchisq(10000, 5, ncp = 0) plot(density(sumsquared), xlim = c(-1, 100)) plot(density(chi), xlim = c(-1, 25)) y=rexp(10000) par(mfrow=c(1,2)) plot(density(y), xlim = c(-1, 20)) hist(y) y1=(rnorm(100000, mean = 0, sd = 1) y2=(rnorm(100000, mean = 0, sd = 1) s=data.frame(cbind(rnorm(10000, mean = 0, sd = 1),rnorm(10000, mean = 0, sd = 1))) len=function(x,y){ ss=x^2+y^2 l=sqrt(ss) l } length=(len(s$x1,s$x2)) par(mfrow=c(1,)) plot(density(length), xlim = c(-1, 8)) a=runif(10000) b=runif(10000) z1=sqrt(-2*log(u))*cos(2*pi*v) z2=sqrt(-2*log(u))*sin(2*pi*v) par(mfrow=c(1,2)) plot(density(z1), xlim = c(-6, 6)) plot(density(z1), xlim = c(-6, 6)) hist(apply(matrix(rcauchy(1000*2), nrow=1000,ncol=2),2,median)) hist(apply(matrix(rcauchy(1000*5), nrow=1000,ncol=5),2,median)) hist(apply(matrix(rcauchy(1000*10), nrow=1000,ncol=10),2,median)) hist(apply(matrix(rcauchy(1000*20), nrow=1000,ncol=20),2,median)) hist(apply(matrix(rcauchy(1000*2), nrow=1000,ncol=2),2,mean)) hist(apply(matrix(rcauchy(1000*5), nrow=1000,ncol=5),2,mean)) hist(apply(matrix(rcauchy(1000*10), nrow=1000,ncol=10),2,mean)) hist(apply(matrix(rcauchy(1000*20), nrow=1000,ncol=20),2,mean)) 0 20 40 60 80 100 0. 00 0. 01 0. 02 0. 03 0. 04 density.default(x = sumsquared) n = 10000 bandwidth = 3.26 d en si ty 020 40 60 80100 0 . 0 0 0 . 0 1 0 . 0 2 0 . 0 3 0 . 0 4 density.default(x = sumsquared) n = 10000 bandwidth = 3.26 d e n s i t y 0 5 10 15 20 25 0. 00 0. 05 0. 10 0. 15 density.default(x = chi) n = 10000 bandwidth = 0.425 d en si ty 0510152025 0 . 0 0 0 . 0 5 0 . 1 0 0 . 1 5 density.default(x = chi) n = 10000 bandwidth = 0.425 d e n s i t y 0 5 10 15 20 0. 0 0. 2 0. 4 0. 6 0. 8 density.default(x = y) n = 10000 bandwidth = 0.1184 d en si ty histogram of y y fr eq ue nc y 0 2 4 6 8 10 0 10 00 20 00 30 00 40 00 50 00 60 00 05101520 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 density.default(x = y) n = 10000 bandwidth = 0.1184 d e n s i t y histogram of y y f r e q u e n c y 0246810 0 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 0 2 4 6 8 0. 0 0. 1 0. 2 0. 3 0. 4 0. 5 0. 6 density.default(x = length) n = 10000 bandwidth = 0.0955 d en si ty 0 2 4 6 8 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 density.default(x = length) n = 10000 bandwidth = 0.0955 d e n s i t y -6 -4 -2 0 2 4 6 0. 0 0. 1 0. 2 0. 3 0. 4 density.default(x = z1) n = 10000 bandwidth = 0.1401 d en si ty -6 -4 -2 0 2 4 6 0. 0 0. 1 0. 2 0. 3 0. 4 density.default(x = z1) n = 10000 bandwidth = 0.1401 d en si ty -6-4-20246 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 density.default(x = z1) n = 10000 bandwidth = 0.1401 d e n s i t y -6-4-20246 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 density.default(x = z1) n = 10000 bandwidth = 0.1401 d e n s i t y histogram of apply(matrix(rcauchy(1000 * 2), nrow = 1000, ncol = 2), 2, median) apply(matrix(rcauchy(1000 * 2), nrow = 1000, ncol = 2), 2, median) fr eq ue nc y 0.03 0.04 0.05 0.06 0.07 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 histogram of apply(matrix(rcauchy(1000 * 2), nrow = 1000, ncol = 2), 2, median) apply(matrix(rcauchy(1000 * 2), nrow = 1000, ncol = 2), 2, median) f r e q u e n c y 0.03 0.04 0.05 0.06 0.07 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 histogram of apply(matrix(rcauchy(1000 * 5), nrow = 1000, ncol = 5), 2, median) apply(matrix(rcauchy(1000 * 5), nrow = 1000, ncol = 5), 2, median) fr eq ue nc y -0.10 -0.05 0.00 0.05 0.10 0. 0 0. 5 1. 0 1. 5 2. 0 histogram of apply(matrix(rcauchy(1000 * 5), nrow = 1000, ncol = 5), 2, median) apply(matrix(rcauchy(1000 * 5), nrow = 1000, ncol = 5), 2, median) f r e q u e n c y -0.10 -0.05 0.00 0.05 0.10 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 histogram of apply(matrix(rcauchy(1000 * 10), nrow = 1000, ncol = 10), 2, median) apply(matrix(rcauchy(1000 * 10), nrow = 1000, ncol = 10), 2, median) fr eq ue nc y -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0. 0 0. 5 1. 0 1. 5 2. 0 2. 5 3. 0 histogram of apply(matrix(rcauchy(1000 * 10), nrow = 1000, ncol = 10), 2, median) apply(matrix(rcauchy(1000 * 10), nrow = 1000, ncol = 10), 2, median) f r e q u e n c y -0.020.000.020.040.060.080.10 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 2 . 5 3 . 0 histogram of apply(matrix(rcauchy(1000 * 20), nrow = 1000, ncol = 20), 2, median) apply(matrix(rcauchy(1000 * 20), nrow = 1000, ncol = 20), 2, median) fr eq ue nc y -0.15 -0.10 -0.05 0.00 0.05 0.10 0 2 4 6 8 histogram of apply(matrix(rcauchy(1000 * 20), nrow = 1000, ncol = 20), 2, median) apply(matrix(rcauchy(1000 * 20), nrow = 1000, ncol = 20), 2, ff(y)="" g="c('-'," '+')="" points(y,="" z,="" pch="g[a" +="" 1])="" y1="rnorm(10000," mean="0," sd="1)^2" y2="rnorm(10000," mean="1," sd="2)^2" y3="rnorm(10000," mean="2," sd="2)^2" y4="rnorm(10000," mean="3," sd="3)^2" y5="rnorm(10000," mean="4," sd="2)^2" sumsquared="cbind(y1,y2,y3,y4,y5)" sumsquared="cbind(apply(Sumsquared,1,sum))" sumsquared="apply(matrix(rnorm(10000," mean="0," sd="1)^2*5,10000,5),1,sum)" chi="rchisq(10000," 5,="" ncp="0)" plot(density(sumsquared),="" xlim="c(-1," 100))="" plot(density(chi),="" xlim="c(-1," 25))="" y="rexp(10000)" par(mfrow="c(1,2))" plot(density(y),="" xlim="c(-1," 20))="" hist(y)="" y1="(rnorm(100000," mean="0," sd="1)" y2="(rnorm(100000," mean="0," sd="1)" s="data.frame(cbind(rnorm(10000," mean="0," sd="1),rnorm(10000," mean="0," sd="1)))" len="function(x,y){" ss="x^2+y^2" l="sqrt(ss)" l="" }="" length="(len(S$X1,S$X2))" par(mfrow="c(1,))" plot(density(length),="" xlim="c(-1," 8))="" a="runif(10000)" b="runif(10000)" z1="sqrt(-2*log(U))*cos(2*pi*V)" z2="sqrt(-2*log(U))*sin(2*pi*V)" par(mfrow="c(1,2))" plot(density(z1),="" xlim="c(-6," 6))="" plot(density(z1),="" xlim="c(-6," 6))="" hist(apply(matrix(rcauchy(1000*2),="" nrow="1000,ncol=2),2,median))" hist(apply(matrix(rcauchy(1000*5),="" nrow="1000,ncol=5),2,median))" hist(apply(matrix(rcauchy(1000*10),="" nrow="1000,ncol=10),2,median))" hist(apply(matrix(rcauchy(1000*20),="" nrow="1000,ncol=20),2,median))" hist(apply(matrix(rcauchy(1000*2),="" nrow="1000,ncol=2),2,mean))" hist(apply(matrix(rcauchy(1000*5),="" nrow="1000,ncol=5),2,mean))" hist(apply(matrix(rcauchy(1000*10),="" nrow="1000,ncol=10),2,mean))" hist(apply(matrix(rcauchy(1000*20),="" nrow="1000,ncol=20),2,mean))" 0="" 20="" 40="" 60="" 80="" 100="" 0.="" 00="" 0.="" 01="" 0.="" 02="" 0.="" 03="" 0.="" 04="" density.default(x="Sumsquared)" n="10000" bandwidth="3.26" d="" en="" si="" ty="" 020="" 40="" 60="" 80100="" 0="" .="" 0="" 0="" 0="" .="" 0="" 1="" 0="" .="" 0="" 2="" 0="" .="" 0="" 3="" 0="" .="" 0="" 4="" density.default(x="Sumsquared)" n="10000" bandwidth="3.26" d="" e="" n="" s="" i="" t="" y="" 0="" 5="" 10="" 15="" 20="" 25="" 0.="" 00="" 0.="" 05="" 0.="" 10="" 0.="" 15="" density.default(x="chi)" n="10000" bandwidth="0.425" d="" en="" si="" ty="" 0510152025="" 0="" .="" 0="" 0="" 0="" .="" 0="" 5="" 0="" .="" 1="" 0="" 0="" .="" 1="" 5="" density.default(x="chi)" n="10000" bandwidth="0.425" d="" e="" n="" s="" i="" t="" y="" 0="" 5="" 10="" 15="" 20="" 0.="" 0="" 0.="" 2="" 0.="" 4="" 0.="" 6="" 0.="" 8="" density.default(x="Y)" n="10000" bandwidth="0.1184" d="" en="" si="" ty="" histogram="" of="" y="" y="" fr="" eq="" ue="" nc="" y="" 0="" 2="" 4="" 6="" 8="" 10="" 0="" 10="" 00="" 20="" 00="" 30="" 00="" 40="" 00="" 50="" 00="" 60="" 00="" 05101520="" 0="" .="" 0="" 0="" .="" 2="" 0="" .="" 4="" 0="" .="" 6="" 0="" .="" 8="" density.default(x="Y)" n="10000" bandwidth="0.1184" d="" e="" n="" s="" i="" t="" y="" histogram="" of="" y="" y="" f="" r="" e="" q="" u="" e="" n="" c="" y="" 0246810="" 0="" 1="" 0="" 0="" 0="" 2="" 0="" 0="" 0="" 3="" 0="" 0="" 0="" 4="" 0="" 0="" 0="" 5="" 0="" 0="" 0="" 6="" 0="" 0="" 0="" 0="" 2="" 4="" 6="" 8="" 0.="" 0="" 0.="" 1="" 0.="" 2="" 0.="" 3="" 0.="" 4="" 0.="" 5="" 0.="" 6="" density.default(x="Length)" n="10000" bandwidth="0.0955" d="" en="" si="" ty="" 0="" 2="" 4="" 6="" 8="" 0="" .="" 0="" 0="" .="" 1="" 0="" .="" 2="" 0="" .="" 3="" 0="" .="" 4="" 0="" .="" 5="" 0="" .="" 6="" density.default(x="Length)" n="10000" bandwidth="0.0955" d="" e="" n="" s="" i="" t="" y="" -6="" -4="" -2="" 0="" 2="" 4="" 6="" 0.="" 0="" 0.="" 1="" 0.="" 2="" 0.="" 3="" 0.="" 4="" density.default(x="Z1)" n="10000" bandwidth="0.1401" d="" en="" si="" ty="" -6="" -4="" -2="" 0="" 2="" 4="" 6="" 0.="" 0="" 0.="" 1="" 0.="" 2="" 0.="" 3="" 0.="" 4="" density.default(x="Z1)" n="10000" bandwidth="0.1401" d="" en="" si="" ty="" -6-4-20246="" 0="" .="" 0="" 0="" .="" 1="" 0="" .="" 2="" 0="" .="" 3="" 0="" .="" 4="" density.default(x="Z1)" n="10000" bandwidth="0.1401" d="" e="" n="" s="" i="" t="" y="" -6-4-20246="" 0="" .="" 0="" 0="" .="" 1="" 0="" .="" 2="" 0="" .="" 3="" 0="" .="" 4="" density.default(x="Z1)" n="10000" bandwidth="0.1401" d="" e="" n="" s="" i="" t="" y="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 2),="" nrow="1000," ncol="2)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 2),="" nrow="1000," ncol="2)," 2,="" median)="" fr="" eq="" ue="" nc="" y="" 0.03="" 0.04="" 0.05="" 0.06="" 0.07="" 0.="" 0="" 0.="" 2="" 0.="" 4="" 0.="" 6="" 0.="" 8="" 1.="" 0="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 2),="" nrow="1000," ncol="2)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 2),="" nrow="1000," ncol="2)," 2,="" median)="" f="" r="" e="" q="" u="" e="" n="" c="" y="" 0.03="" 0.04="" 0.05="" 0.06="" 0.07="" 0="" .="" 0="" 0="" .="" 2="" 0="" .="" 4="" 0="" .="" 6="" 0="" .="" 8="" 1="" .="" 0="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 5),="" nrow="1000," ncol="5)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 5),="" nrow="1000," ncol="5)," 2,="" median)="" fr="" eq="" ue="" nc="" y="" -0.10="" -0.05="" 0.00="" 0.05="" 0.10="" 0.="" 0="" 0.="" 5="" 1.="" 0="" 1.="" 5="" 2.="" 0="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 5),="" nrow="1000," ncol="5)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 5),="" nrow="1000," ncol="5)," 2,="" median)="" f="" r="" e="" q="" u="" e="" n="" c="" y="" -0.10="" -0.05="" 0.00="" 0.05="" 0.10="" 0="" .="" 0="" 0="" .="" 5="" 1="" .="" 0="" 1="" .="" 5="" 2="" .="" 0="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 10),="" nrow="1000," ncol="10)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 10),="" nrow="1000," ncol="10)," 2,="" median)="" fr="" eq="" ue="" nc="" y="" -0.02="" 0.00="" 0.02="" 0.04="" 0.06="" 0.08="" 0.10="" 0.="" 0="" 0.="" 5="" 1.="" 0="" 1.="" 5="" 2.="" 0="" 2.="" 5="" 3.="" 0="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 10),="" nrow="1000," ncol="10)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 10),="" nrow="1000," ncol="10)," 2,="" median)="" f="" r="" e="" q="" u="" e="" n="" c="" y="" -0.020.000.020.040.060.080.10="" 0="" .="" 0="" 0="" .="" 5="" 1="" .="" 0="" 1="" .="" 5="" 2="" .="" 0="" 2="" .="" 5="" 3="" .="" 0="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 20),="" nrow="1000," ncol="20)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 20),="" nrow="1000," ncol="20)," 2,="" median)="" fr="" eq="" ue="" nc="" y="" -0.15="" -0.10="" -0.05="" 0.00="" 0.05="" 0.10="" 0="" 2="" 4="" 6="" 8="" histogram="" of="" apply(matrix(rcauchy(1000="" *="" 20),="" nrow="1000," ncol="20)," 2,="" median)="" apply(matrix(rcauchy(1000="" *="" 20),="" nrow="1000," ncol="20),">