Course: STA6934- Monte Carlo Statistical MethodsInstructor: Professor CasellaAssignment 1Problem 1.1mix <- function(x, e, u1, u2, sigma1, sigma2)(e*dnorm(x, mean=u1, sd=sigma1)+ (1-e)*dnorm(x, mean=u2, sd=sigma2))min <- function(x, u1, u2, sigma1, sigma2)((1-pnorm((x-u1)/sigma1))/sigma2*dnorm((x-u2)/sigma2)+ (1-pnorm((x-u2)/sigma2))/sigma1*dnorm((x-u1/sigma1)))censplot <- function(e,u1,u2,sigma1,sigma2){lowpoint=pmin(u1,u2)-3*pmax(sigma1,sigma2)uppoint=pmax(u1,u2)+3*pmax(sigma1,sigma2)xplot <- seq(from=lowpoint, to=uppoint, length=1000)mixxplot <- mix(xplot, e, u1, u2, sigma2, sigma2)minxplot <- min(xplot, u1, u2, sigma1, sigma2)plot(xplot, mixxplot, xlim=c(lowpoint, uppoint), ylim=c(0,0.8),type="l", lty=1, ylab="density", col="blue")lines(xplot, minxplot, lty=2, col="red")legend(lowpoint, 0.8, c("mixed", "minimum"), lty=c(1,2), col=c("blue", "red"))mtext(bquote(paste("u=", .(e), ",Normal(", .(u1),",", .(sigma1),"),Normal(", .(u2),",", .(sigma2), ")")))}#library(lattice)#trellis.device(pdf, file="HW1p1", height=20, width=17)par(mfrow=c(3,2))censplot(0.3,1,1,1,1)censplot(0.3,-1,1,1,1)censplot(0.3,1,1,2,1)censplot(0.3,-3,1,3,1)censplot(0.5,3,1,2,1)censplot(0.5,-3,1,1,3)1Problem 1.4 In order to find an explicit form of the integralZ∞ωαβxα−1e−βxαdx,we use the change of variable y = xα. We have dy = αxα−1dx and the integral becomesZ∞ωαβxα−1e−βxαdx =Z∞ωαβe−βydy = e−βωα.Problem 1.7 The density f of the vector Ynisf(yn, µ, σ) =1σ√2πnexp −12nXi=1yi− µσ2!, ∀yn∈ Rn, ∀(µ, σ2) ∈ R ×R∗+This function is strictly positive and the first and second order partial derivatives with respect to µ andσ exist and are positive. The same hypotheses are satisfied for the log-likelihood functionlog(L(µ, σ, yn)) = −n log√2π −n log σ −12nXi=1yi− µσ2thus we can find the ML estimator of µ and σ2. The gradient of the log-likelihood is∇log (L) =(∂ log(L(µ,σ,yn))∂µ∂ log(L(µ,σ,yn))∂σ=(1σ2Pni=1(yi− µ)−nσ+Pni=1(yi−µ)2σ3if we equate the gradient to the null vector, ∇log (L) = 0 and solve the resulting system in µ and σ, wefindˆµ =1nnXi=1yi= ¯y ,ˆσ2=1nnXi=1(yi− ¯y)2= s2.Problem 1.13For X ∼ We(α, β, γ), where α > 0 is the shape parameter, β > 0 is the scale parameter, and γ isthe translation parameter, the density is given as:f(x; α, β, γ) =αβ(x − γβ)α−1e−(x−γβ)α, for x ≥ γ2For X1, . . . , Xnare iid as We(α, β, γ), the likelihood function is given as:L(α, β, γ|x1, . . . , xn) =nYi=1f(xi|α, β, γ)= (αβ)n·nYi=1(xi− rβ)α−1· exp{−X(xi− rβ)α}l(α, β, γ|x1, . . . , xn) = logL(·)= n[logα − logβ] + (α − 1) ·X[log(xi− r) − logβ] −X(xi− rβ)α= nlogα − nαlogβ + (α − 1) ·Xlog(xi− r) −X(xi− rβ)α∂∂αl(·) =nα− nlogβ +Xlog(xi− r) −X[logxi− rβ· (xi− rβ)α]∂∂βl(·) = −nαβ+Xαβ−(α+1)· (xi− r)α∂∂γl(·) = (α − 1) ·X(−1xi− r) +Xαβα(xi− r)α−11. γ = 100, α = 3n=19, in this case∂∂βl(·) = −19 ∗ 3β+X3(xi− 100)3· β−4and we getˆβ = 125.6846Also we could use the nlm function in R to solve this problem. nlm function carries out a min-imization of the function using a Newton-type algorithm. Thus to find the maximum likelihoodestimator for the parameter is equivalent to get the value minimizing the negative log-likelihood.The R code and the output is given as below:>logweib1 <- function (beta){ alpha=3; size=19;x <- c(143, 164, 188, 188, 190, 192, 206, 209, 213, 216,220, 227, 230, 234, 246, 265, 304, 216, 244)-(size*log(alpha)-size*alpha*log(beta)+(alpha-1)*sum(log(x-100))-sum(((x-100)/beta)ˆalpha))}>nlm(logweib1,100)3$minimum[1] 95.13106$estimate[1] 125.6845As we can see from the above, the two beta values are quite close, hence the fitted model isW eibull(3, 125.6845, 100). We will use this method for the following steps.2. γ = 100, α unknown;>logweib2 <- function (p){ size=19;x <- c(143, 164, 188, 188, 190, 192, 206, 209, 213, 216,220, 227, 230, 234, 246, 265, 304, 216, 244)-(size*log(p[1])-size*p[1]*log(p[2])+(p[1]-1)*sum(log(x-100))-sum(((x-100)/p[2])ˆp[1]))}>nlm(logweib2, c(3, 125.6845)) # use the result from (a) as the initial value$minimum[1] 94.75059$estimate[1] 3.504232 128.104290Hence the fitted model is W eibull(3.504232, 128.104290, 100).3. γ and α both unknown;>logweib3 <- function (p){ size=19;x <- c(143, 164, 188, 188, 190, 192, 206, 209, 213, 216,220, 227, 230, 234, 246, 265, 304, 216, 244)-(size*log(p[1])-size*p[1]*log(p[2])+(p[1]-1)*sum(log(x-p[3]))4-sum(((x-p[3])/p[2])ˆp[1]))}>nlm(logweib3,c(3.504232, 128.104290, 100))# use the result from (b) as the initial value$minimum[1] 94.59973$estimate[1] 2.849366 105.345531 121.425922Hence the fitted model is W eibull(2.849366, 105.345531, 121.425922).Problem 1.22(a). Since L(δ, h(θ)) ≥ 0 by using Fubini’s theorem, we getr(π, δ) =ZΘZXL(δ, h(θ))f(x|θ)π(θ)dxdθ=ZXZΘL(δ, h(θ))f(x|θ)π(θ)dθdx=ZXZΘL(δ, h(θ))m(x)π(θ|x)dθdx=ZXϕ(π, δ|x)m(x)dx ,where m is the marginal distribution of X and ϕ(π, δ|x) is the posterior average cost.The estimator that minimizes the integrated risk r is therefore, for each x, the one that minmizesthe posterior average cost and it is given byδπ(x) = arg minδϕ(π, δ|x) .(b). The average posterior loss is given by :ϕ(π, δ|x) = Eπ[L(δ, θ)|x]= Eπ||h(θ) − δ||2|x = Eπ||h(θ)||2|x + δ2− 2 < δ, Eπ[h(θ)|x] >A simple derivation shows that the minimum is attained forδπ(x) = Eπ[h(θ)|x] .5(c). Take m to be the posterior median and consider the auxiliary function of θ, s(θ), defined ass(θ) =(−1 if h(θ) < m+1 if h(θ) > mThen s satisfies the proprietyEπ[s(θ)|x] = −Zm−∞π(θ|x)dθ +Z∞mπ(θ|x)dθ= −P(h(θ) < m|x) + P(h(θ) > m|x) = 0For δ < m, we have L(δ, θ) − L(m, θ) = |h(θ) − δ| − |h(θ) − m| from which it follows thatL(δ, θ) − L(m, θ) =δ −m = (m − δ)s(θ) if δ > h(θ)m − δ = m − δ if m < δ2h(θ) − (δ + m) > (m − δ)s(θ) if δ < h(θ) < mIt turns out that L(δ, θ) − L(m, θ) > (m − δ)s(θ) which implies thatEπ[L(δ, θ) − L(m, θ)|x] > (m − δ)Eπ[s(θ)|x] = 0 .This still holds, using similar argument when δ > m, so the minimum of Eπ[L(δ, θ)|x] is reachedat δ = m.Problem 1.23(a). When X|σ ∼ N(0, σ2),1σ2∼ Ga(1, 2), the posterior distribution isπσ−2|X∝ f(x|σ)π(σ−2)∝1σe−(x2/2+2)σ2=