Statistical Computing HW6

(1) C=1, when \(V_1^2 + V_2^2 \le 1\)

Under the case ‘C = 1’, there are two subcases, including \(S = 1, (V_1^2 + V_2^2 + V_3^2 \le 1)\) given \(C = 1\), and \(S = 0, (V_1^2 + V_2^2 + V_3^2 > 1)\) given \(C = 1\).

The first subcase ‘S = 1 given C = 1’ has no special meaning because if \(V_1^2 + V_2^2 + V_3^2 \le 1\) holds, \(V_1^2 + V_2^2 \le 1\) is always true. In other words, \(\{S = 1\} \subset \{C = 1\}\), which leads the probability \(\mathbb P(S = 1,C=1)\) to \(\mathbb P(S=1)\). Then, \[\begin{aligned} \mathbb P(S=1|C=1) &= \mathbb P(S = 1,C=1)/\mathbb P(C=1) \\ & = \mathbb P(S=1)/\mathbb P(C=1) \\ &= (\pi/6)/(\pi/4) = 2/3. \end{aligned}\]

The second subcase ‘S = 0 given C = 1’ implies the points that are outside of the sphere (with radius 1) and inside the cylinder (with radius 1 and height 2). Then, the volume of such space is \(2\pi*1^2(\text{volumn of cylinder}) - \frac{4}{3}\pi*1^3(\text{volumn of shpere}) = \frac{2}{3}\pi\), and \(\mathbb P(S=0,C=1) = \frac{2}{3}\pi/8 = \pi/12\). Thus, \[\mathbb P(S=0|C=1) = \mathbb P(S=0,C=1)/\mathbb P(C=1) = (\pi/12)/(\pi/4) = 1/3.\]

(2) C = 0: when \(V_1^2 + V_2^2 > 1\)

There are two subcases, S = 1 given C = 0, and S = 0 given C = 0.

The first subcase, S = 1 given C = 0, has no special meaning, because when \(V_1^2 + V_2^2 > 1\), there are no such case that leads \(V_1^2 + V_2^2 +V_3^2\)to be less than or equal to 1. In other words, \(\{S=1\} \cap \{C=0\} = \emptyset\), which leads to \(\mathbb P(S=1,C=0) = 0\). So, \[\mathbb P(S=1|C=0) = \mathbb P(S=1,C=0)/\mathbb P(C=0) = 0.\]

The second subcase, S = 0 given C = 0, means the points outside the cylinder and inside the cube (the length of edge is 2). Then, the volume of such space is \((4-\pi)*2\), which leads the probability \(\mathbb P(S=0,C=1)\) to \((4-\pi)*2/8 = 1-\pi/4\). Thus, \[\mathbb P(S=0|C=1) = \mathbb P(S=0,C=1)/\mathbb P(C=0) = (1-\pi/4)/(1-\pi/4) = 1.\]

(3) conditional expectation for each C

\[\begin{aligned} \mathbb E[S|C=1] &= \sum_{s\in\{0,1\}}s*\mathbb P(S=s|C=1) \\ &= 0*(\frac{1}{3}) + 1*(\frac{2}{3}) = \frac{2}{3}. \end{aligned}\]

Note that 2/3 can be expressed as (2/3)*1.

\[\begin{aligned} \mathbb E[S|C=0] &= \sum_{s\in\{0,1\}}s*\mathbb P(S=s|C=0) \\ &= 0*1 + 1*0 = 0. \end{aligned}\]

Note that 0 can be expressed as (2/3)*0.

(4) conclusion

To sum up, \(\mathbb E[S|C=c] = (2/3)*c\) for \(c = 0,1\). So, \(\mathbb E[S|C] = (2/3)*C\) Multiplying the integer 6 to each side, the following expression can be induced. \[4C = \mathbb E[6S|C].\]

The random variable C is identical to the indicator \(I(V_1^2 + V_2^2 \le 1)\), and the random variable S is same as \(I(V_1^2 + V_2^2 + V_3^2 \le 1)\). Therefore, the proper condition is the ‘C’ or the indicator\(I(V_1^2 + V_2^2 \le 1)\), and the expression above can be shows with respect to indicators as follows: \[4*I(V_1^2 + V_2^2 \le 1) = \mathbb E[6*I(V_1^2 + V_2^2 + V_3^2 \le 1)|I(V_1^2 + V_2^2 \le 1)].\]

i) Show that \(\mathbb E_Y[\hat \theta] = \mathbb P(X \in A)\), where \(X \sim F\).

\[\begin{aligned} \mathbb P(X \in A) &= \mathbb E_X [I(X\in A)]\\ &= \int_{-\infty}^{\infty} I(x \in A)f(x)dx &\qquad ①\\ &= \int_{-\infty}^{\infty} I(x \in A)\frac{f(x)}{g(x)}g(x)dx &\qquad ② \end{aligned}\]

The expression (1) is the expectation of a function of a variable that follows \(F\), but the expression (2) is the expectation of a function of a variable that follows \(G\). This problem denotes such variable as \(Y \sim G\), so changing ‘x’ into ‘y’ would be better notation. Then, \[\int_{-\infty}^{\infty} I(y \in A)\frac{f(y)}{g(y)}g(y)dy = \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)}\right].\]

Due to the law of total expectation, the below expression holds. \[\mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)}\right] \\ = \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)} \Big| Y \in A\right]\mathbb P(Y \in A) + \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)} \Big| Y \in A^c\right]\mathbb P(Y \in A^c).\]

When \(Y \in A^c\), the indicator \(I(Y \in A)\) becomes 0. Then, the latter conditional expectation \(\mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)} \Big| Y \in A^c\right]\) becomes the expectation of 0’s. Then, \[\begin{aligned} \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)}\right] &= \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)} \Big| Y \in A\right]\mathbb P(Y \in A) + 0 \\ &= \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)} \Big| Y \in A\right]*p &\qquad ③. \end{aligned}\]

The indicator function \(I(Y \in A)\) is always 1 when \(Y \in A\). According to the definition of \(\hat \theta\), each \(Y-i\) follows the distribution \(G_A\), the truncated distribution of \(G\). Then, every \(Y_i\) is included in \(A\). In other words, \(I(Y_i \in A)\) is always 1. After choosing any ‘n’ number of \(Y_i\), the above equation becomes \[p* \mathbb E_Y\left[I(Y \in A)*\frac{f(Y)}{g(Y)} \Big| Y \in A\right] = p*\mathbb E_Y\left[1*\frac{f(Y_i)}{g(Y_i)}\right].\]

Each \(Y_i\) is the identical random variable from the distribution \(G_A\), and each \(Y_i\) are independent. Then, each \(E\left[\frac{f(Y_i)}{g(Y_i)}\right]\) are identical and \[\begin{aligned} p*\mathbb E_Y\left[1*\frac{f(Y_i)}{g(Y_i)}\right] &= p*\frac{1}{n}\sum_{i = 1}^n \mathbb E_Y\left[\frac{f(Y_i)}{g(Y_i)}\right] \\ &= \frac{p}{n}\mathbb E_Y\left[\sum_{i = 1}^n\frac{f(Y_i)}{g(Y_i)} \right] \\ &= \mathbb E_Y\left[\frac{p}{n} \sum_{i = 1}^n\frac{f(Y_i)}{g(Y_i)} \right] \\ &= \mathbb E_Y[\hat \theta]. \end{aligned}\]

This proves \(\mathbb E_Y[\hat \theta] = \mathbb P(X \in A)\). Let this value be \(\theta\).

ii) Show that \(\hat \theta\) is constructed by combining the importance sampling, the conditioning, and the stratified sampling.

i) Plugging in the best \(c\) into \(g(U)\), find an expression for the variance of \(g(U)\) as a function of \(a\) and \(k\).

Let us abridge the indicator \(I(U \le a)\) as \(I\) simply from now on. The variance of \(g(U)\) is as follows: \[\begin{aligned} \text{Var}[g(U)] &= \text{Var}[I + c*U^k - c*\mathbb E[U^k]] \\ &= \text{Var}[I + c*U^k] \because \mathbb E[U^k] \text{ is constant} \\ &= \text{Var}[I] + c^2\text{Var}[U^k] + 2c*Cov[I, U^k] \\ &= \text{Var}[I] + \text{Var}[U^k]*\left(c^2 + 2*\frac{Cov[I,U^k]}{\text{Var}[U^k]}*c\right) \\ &= \text{Var}[I] + \text{Var}[U^k]*\left(c + \frac{Cov[I, U^k]}{\text{Var}[U^k]} \right)^2 - \frac{(Cov[I,U^k])^2}{\text{Var}[U^k]}. \end{aligned}\]

The above equation attains minimum when \(c = \frac{(Cov[I,U^k])^2}{\text{Var}[U^k]}\), which is the best c. After applying the best c to the above equation, \(\text{Var}[g(U)] = \text{Var}[I] - \frac{(Cov[I,U^k])^2}{\text{Var}[U^k]}\).

ii) With a given \(a\), we are certainly interested in \(k\) which minimizes the variance of \(g(U)\). Unfortunately, the expression for the variance involves both an exponential term and a polynomial term in \(k\), and hence the optimal \(k\) cannot be obtained analytically. We can still find the optimal \(k\) in a numerical way. For each of \(a = 0.2, 0.5, 0.8\), use the ‘optim’ function to find the optimal \(k\). Find the variance of \(g(U)\) for each \(‘a’\) and the corresponding \(k\), and calculate the percentage decrease compared to using \(I(U \le a)\).

Before using the ‘optim’, a function named ‘fn’ is defined to use it as a parameter of ‘optim’ function. This function works as the variance of \(g(U)\) for the given value ‘a’.

fn <- function(a, par){
  k <- par
  return(a*(1-a) - 1/k^2*(2*k+1)*a^2*(a^k - 1)^2)
}

The ‘optim’ function requires some parameters to operate. The ‘par’ is the initial parameter to be calculated, and the ‘fn’ is the function where optimization (in this case, the minimization) occurs.

For reference, the ‘optim’ suggests “Brent” method as a parameter for one dimensional problem, which requires minimum and maximum of the interval of k. However, setting any number that makes the expectation possible is allowed according to the notice. So, setting the initial ‘par’ to 1 was chosen, although it returns a warning with the optimal k and minimum variance. Note that returning a warning make no problem according to the notice.

The result of ‘optim’ when a = 0.2 is saved in ‘op.2’. The optimal k value when a = 0.2 (k.2) is 0.1466797 and the corresponding minimum variance of g(U) is 0.05367902. The value of k and the minimal variance is in the op.2$par and op.2$value respectively. Since the variance of \(I(U \le a)\) is \(a(1-a)\), \(\text{Var}[I(U \le 0.2)] = 0.2*0.8 = 0.16\). Thus, the percentage of variance reduction is (0.16 – 0.05367902)/0.16 = 66.45%.

op.2 <- optim(par = 1, fn = fn, a = 0.2)

## Warning in optim(par = 1, fn = fn, a = 0.2): Nelder-Mead를 이용한 1차원 최적화 문제는 신뢰할 수 없습니다:
## "Brent" 또는 optimize()를 이용해보세요

op.2$value

## [1] 0.05367902

k.2 <- op.2$par
k.2

## [1] 0.1466797

0.2*(1-0.2)

## [1] 0.16

(0.2*(1-0.2) - op.2$value)/(0.2*(1-0.2)) # 66.45%

## [1] 0.6645061

The result of ‘optim’ when a = 0.5 is saved in ‘op.5’. The optimal k value when a = 0.5 (k.5) is 1.164062 and the corresponding minimum variance of g(U) is 0.06171835. The value of k and the minimal variance is in the op.5$par and op.5$value respectively. The variance of \(I(U \le 0.5)\) is \(0.5*0.5 = 0.25\), so the percentage of variance reduction is (0.25 – 0.0617835)/0.25 = 75.31%

op.5 <- optim(par = 1, fn = fn, a = 0.5)

## Warning in optim(par = 1, fn = fn, a = 0.5): Nelder-Mead를 이용한 1차원 최적화 문제는 신뢰할 수 없습니다:
## "Brent" 또는 optimize()를 이용해보세요

op.5$value

## [1] 0.06171835

k.5 <- op.5$par
k.5

## [1] 1.164062

0.5*(1-0.5)

## [1] 0.25

(0.5*(1-0.5) - op.5$value)/(0.5*(1-0.5)) # 75.31%

## [1] 0.7531266

The result of ‘optim’ when a = 0.8 is saved in ‘op.8’. The optimal k value when a = 0.8 (k.8) is 4.97348 and the corresponding minimum variance of g(U) is 0.03271169. The value of k and the minimal variance is in the op.8$par and op.8$value respectively. The variance of \(I(U \le 0.8)\) is 0.8*0.2 = 0.16, so the percentage of variance reduction is (0.16 – 0.03271169)/0.16 = 79.555%

op.8 <- optim(par = 1, fn = fn, a = 0.8)

## Warning in optim(par = 1, fn = fn, a = 0.8): Nelder-Mead를 이용한 1차원 최적화 문제는 신뢰할 수 없습니다:
## "Brent" 또는 optimize()를 이용해보세요

op.8$value

## [1] 0.03271169

k.8 <- op.8$par
k.8

## [1] 4.973438

0.8*(1-0.8)

## [1] 0.16

(0.8*(1-0.8) - op.8$value)/(0.8*(1-0.8)) # 79.555%

## [1] 0.7955519

iii) For each of \(a = 0.2, 0.5, 0.8\) with the corresponding optimal \(k\), overlay graphs of \(I(U \le a)\) and \(u^k\) as functions of \(u, 0 < u < 1\), in order to see the relationship between the binary indicator \(I(U \le a)\) and the control variable \(U^k\). Also, as mentioned in class, another way of doing this is comparing the graph of \(I(V \le a^k)\) with a straight line by transforming 𝑣 = 𝑢𝑘. Draw the graph of 𝐼(𝑣 ≤𝑎𝑘), 0 < v < 1. Based on the two plots, explain why such a k must be obtained as the optimal value for each a.

To start with, a function that works as an indicator function has to be defined. It receives a vector and a number serving as a criterion whether the element of the vector is less than or equal to the number. If the element fits the criterion, it returns 1, and if not, it returns 0. The following code shows the function (indicator.a) and how it works.

indicator.a <- function(us, a){
  I <- rep(0,length(us))
  I[us <= a] <- 1
  return(I)
}
indicator.a(seq(0,1,0.1),0.2)

##  [1] 1 1 1 0 0 0 0 0 0 0 0

The following code shows how to make the plots of indicator function (unit step function) and the u^k. The plot shows when a = 0.2, a = 0.5, a = 0.8 from the left. The red, green, and blue graph in each plot means \(y = u^k\), where the value k is the optimal value of k for each a.

u <- seq(0,1,0.01)
par(mfrow = c(1,3))
plot(u, indicator.a(u, 0.2), type = 'l', ylab = 'y',
     main = 'plot of indicator and u^k when a = 0.2')
lines(u, u^k.2, col = 'red')

plot(u, indicator.a(u, 0.5), type = 'l', ylab = 'y',
     main = 'plot of indicator and u^k when a = 0.5')
lines(u, u^k.5, col = 'green')

plot(u, indicator.a(u, 0.8), type = 'l', ylab = 'y',
     main = 'plot of indicator and u^k when a = 0.8')
lines(u, u^k.8, col = 'blue')

The optimal k values for each a are as follows: k.2 = 0.147, k.5 = 1.16, and k.8 = 4.97, defined in previous question. The reason for the k being chosen is that such k can reduce the variance of g(U) as much as possible. The equation below can be obtained by manipulating the form of \(\text{Var}[g(U)]\): \[\text{Var}[g(U)] = \text{Var}[I] - \frac{(Cov[I,U^k])^2}{\text{Var} [U^k]} = \text{Var}[I] - \text{Var}[I]\left(Cov[I,U^k] \right)^2\]

The value of \(\text{Var}[I]\) is fixed for each a, so how much the variance reduces depends on the correlation between \(I\) and \(U^k\).

par(mfrow = c(1,3))
plot(u, 1-indicator.a(u, 0.2), type = 'l', ylab = 'y',
     main = '1-indicator and u^k when a = 0.2')
lines(u, u^k.2, col = 'red')

plot(u, 1-indicator.a(u, 0.5), type = 'l', ylab = 'y',
     main = '1-indicator and u^k when a = 0.5')
lines(u, u^k.5, col = 'green')

plot(u, 1-indicator.a(u, 0.8), type = 'l', ylab = 'y',
     main = '1-indicator and u^k when a = 0.8')
lines(u, u^k.8, col = 'blue')

v <- seq(0,1,0.01)
par(mfrow = c(1,3))
plot(v, indicator.a(v, 0.2^k.2), type = 'l', xlab = 'v', ylab = 'y',
     main = 'plot of indicator w.r.t v when a = 0.2')
lines(v, v, col = 'red')

plot(v, indicator.a(v, 0.5^k.5), type = 'l', xlab = 'v', ylab = 'y',
     main = 'plot of indicator w.r.t v when a = 0.5')
lines(v, v, col = 'green')


plot(v, indicator.a(v, 0.8^k.8), type = 'l', xlab = 'v', ylab = 'y',
     main = 'plot of indicator w.r.t v when a = 0.8')
lines(v, v, col = 'blue')

par(mfrow = c(1,3))
v2 <- seq(0.2^k.2,1,0.01)
ones.2 <- rep(1,length(v2))
plot(v, 1-indicator.a(v, 0.2^k.2), type = 'l', xlab = 'v', ylab = 'y',
     main = '1-indicator and v when a = 0.2')
lines(v, v, col = 'red')
lines(v2,v2, col = 'red', lwd = 4)
lines(v2,ones.2, col = 'black', lwd = 4)


plot(v, 1-indicator.a(v, 0.5^k.5), type = 'l', xlab = 'v', ylab = 'y',
     main = '1-indicator and v when a = 0.5')
lines(v, v, col = 'green', lwd = 2)

v8 <- seq(0,0.8^k.8, 0.01)
ones.8 <- rep(0, length(v8))
plot(v, 1-indicator.a(v, 0.8^k.8), type = 'l', xlab = 'v', ylab = 'y',
     main = '1-indicator and v when a = 0.8')
lines(v, v, col = 'blue')
lines(v8,v8, col = 'blue', lwd = 4)
lines(v8,ones.8, col = 'black', lwd = 4)

stu_id <- 2011142127
set.seed(stu_id)
us <- runif(10^6)
par(mfrow = c(1,3))
hist(us^k.2, freq = FALSE, main = 'histogram of u^k when a = 0.2')
hist(us^k.5, freq = FALSE, main = 'histogram of u^k when a = 0.5')
hist(us^k.8, freq = FALSE, main = 'histogram of u^k when a = 0.8')

iv) Let \(U_i, \ i=1,\dots,n\) be an i.i.d random sample of \(U\). Sampling \(n = 10^6\) uniform random numbers, estimate the variances of the estimators \(n^{-1}\sum_{i=1}^ng(U_i)\) for \(a = 0.2, 0.5, 0.8\). Calculating the estimated variance of \(n^{-1}\sum_{i=1}^nI(U_i \le a)\), estimate the percentage decrease for \(a = 0.2, 0.5, 0.8\). Compare the estimated percentage decreases with the theoretical ones obtained above.

Every built-in basic function in R is available in this homework. So functions such as ‘var’, ‘cov’ will be used.

The following code serves as a function of estimating the variance of \(n^{-1}\sum_{i=1}^ng(U_i)\), the sample mean of \(g(U)\), whose inputs should be a vector (consisting of uniform distribution), value ‘a’ and the corresponding ‘k’. The optimal value of ‘c’ can be obtained by estimating the covariance of \(I(U \le a)\) and \(U^k\) and estimating the variance of \(U^k\). Using the estimated ‘c’, the value of \(g(U_i)\) can be obtained for each \(U_i, \ i=1,\dots,10^6\). These values are saved in the ‘gs.ui’ vector, the length of which is same as the input vector (us). Since the estimated variance of a sample mean (corresponding to \(n^{-1}\sum_{i=1}^ng(U_i)\) is its sample variance divided by the number of elements (\(Var[\bar X] = S^2/n, \ S^2 = \frac{1}{n-1}\sum_{i=1}^n(X_i - \bar X)^2\)), the following code returns the (sample) variance of ‘gs.ui’ divided by n (= the length of ‘us’).

set.seed(stu_id)
n <- 10^6
us <- runif(n)

est_var_g <- function(us, a, k){
  n <- length(us)
  c <- -cov(indicator.a(us, a), us^k)/var(us^k) #estimated c
  gs.ui <- indicator.a(us, a) + c*(us^k - 1/(k+1)) # g(ui) from 1 to n
  S.2 <- var(gs.ui)
  return(S.2/n) #ch8 estimated var of X_bar
}

The variance of \(n^{-1}\sum_{i=1}^nI(U_i \le a)\) can be obtained by using ‘var’ function to an indicator function (indicator.a) and dividing by the number of elements.

When a = 0.2, the estimated variance of \(n^{-1}\sum_{i=1}^ng(U_i)\) is about \(5.35∗10^{−8}\) and the estimated variance of \(n^{-1}\sum_{i=1}^nI(U_i \le a)\) is about \(1.59∗10^{−7}\). Then, the percentage of variance reduction is \((1.59∗10^{−7}- 5.35∗10^{−8})/(1.59∗10^{−7}) \approx 66.41\%\) .

# for the case a = 0.2

var_g_2 <- est_var_g(us, 0.2, k.2) # 5.3539*10^(-8)
var_i_2 <- var(indicator.a(us, 0.2))/n # 1.594*10^(-7)

(var_i_2 - var_g_2)/var_i_2 # 66.41155 %

## [1] 0.6641155

The theoretical percentage of variance reduction in ii) was about 66.45%, so it is shown that the variance reduction using \(n^{-1}\sum_{i=1}^ng(U_i)\) is slightly lower than ideal one.

When a = 0.5, the estimated variance of \(n^{-1}\sum_{i=1}^ng(U_i)\) is about \(6.18∗10^{−8}\) and the estimated variance of \(n^{-1}\sum_{i=1}^nI(U_i \le a)\) is about \(2.50∗10^{−7}\). Then, the percentage of variance reduction is \((2.50∗10^{−7} - 6.18∗10^{−8})/(2.50∗10^{−7}) \approx 75.28\%\) .

# for the case a = 0.5

var_g_5 <- est_var_g(us, 0.5, k.5) # 6.18*10^(-8)
var_i_5 <- var(indicator.a(us, 0.5))/n # 2.500*10^(-7)

(var_i_5 - var_g_5)/var_i_5 # 75.27761 %

## [1] 0.7527761

The theoretical percentage of variance reduction in ii) was about 75.31%, so it is shown that the variance reduction using \(n^{-1}\sum_{i=1}^ng(U_i)\) is slightly lower than ideal one too.

When a = 0.8, the estimated variance of \(n^{-1}\sum_{i=1}^ng(U_i)\) is about \(3.27∗10^{−8}\) and the estimated variance of \(n^{-1}\sum_{i=1}^nI(U_i \le a)\) is about \(1.60∗10^{−7}\). Then, the percentage of variance reduction is \((1.60∗10^{−7} - 3.27∗10^{−8})/(1.60∗10^{−7}) \approx 79.546\%\).

# for the case a = 0.8

var_g_8 <- est_var_g(us, 0.8, k.8) # 3.27*10^(-8)
var_i_8 <- var(indicator.a(us, 0.8))/n # 1.59887*10^(-7)

(var_i_8 - var_g_8)/var_i_8 # 79.54615 %

## [1] 0.7954615

The theoretical percentage of variance reduction in ii) was about 79.555%, so it is shown that the variance reduction using \(n^{-1}\sum_{i=1}^ng(U_i)\) is slightly lower than ideal one too.

The reason for the estimated variance reduction being slightly lower than the ideal case is that the estimated one used the estimated optimal value of c, which would be different from ideal ‘c’. Note that the variance form of control variable is \[\text{Var}[I] + \text{Var}[U^k]*\left(c + \frac{Cov[I, U^k]}{\text{Var}[U^k]} \right)^2 - \frac{(Cov[I,U^k])^2}{\text{Var}[U^k]},\] the quadratic function of ‘c’. So, the value of c which is slightly different from the optimal \(c = -\frac{Cov[I, U^k]}{\text{Var}[U^k]}\) would make the variance slightly bigger than the minimum variance.

i) © We want to estimate 𝑃[ℎ(𝑈) ≤ −19.1]. A control variable estimator can be constructed by \[g(U) = I(h(U) \le -19.1) + c*(U^k - \mathbb E[U^k]).\] Now, an expression for the variance of \(g(U)\) is not available, and hence the optimal k is not easily obtained. However, we can still utilize the informal graphical approach in Q3 iii) (both types of graphs or only one). Using a graphical analysis through trial and error, find a reasonable \(k\).

To analyze it with graphical approach, the indicator function \(I(h(U) \le -19.1)\) should be defined. The input is a vector, and ‘h’ is defined to represent the h(u) function. For the indices where the values of ‘h’ is less than or equal to -19.1, the I vector filled with all 0’s will become ‘1’. The output will be the shifted ‘I’ vector. The plot below shows a function similar to the unit step function made by the indicator.h.

indicator.h <- function(us){
  I <- rep(0,length(us))
  
  h <- sin(20*pi*us)+20*sin(pi*(us-1)/2)
  I[h <= -19.1] <- 1
  
  return(I)
}
par(mfrow = c(1,1))

uh <- seq(0,1,1/100000)
plot(uh, indicator.h(uh), type = 'l', ylab = 'y', xlab = 'u',
     main = 'plot of indicator I(h(u) <= -19.1)',cex.main = 1)

Except for some points at \(0 < u < 2\), most of u in the interval make h(u) to be ‘1’. Analyzing the plot, it is assumed that the length of u which make h(u) ‘1’ is about 0.2. Then, this plot may be regarded as similar graph with the unit step function in Q3 iii) when a = 0.2.

The main point of this question is to find a proper ‘k’ that makes the variance of g(U) considerably small. As previously mentioned, the optimal k makes the absolute value of correlation between the indicator function \(I(h(U) \le -19.1)\) and \(U^k\) as large as possible. So, the following process attempts to find the best ‘k’ from the k.2 = 0.147, the value k when a = 0.2 in question i) through iv). As the k decreases (0.147 -> 0.14 -> 0.12), the correlation decreases as well. After letting the k increase, the correlation reaches maximum at k = 0.16. and decreases again. The maximum correlation is estimated as 0.6787737.

abs(cor(indicator.h(us), us^k.2)) # k.2 ~ 0.147

## [1] 0.678647

abs(cor(indicator.h(us), us^(0.14)))

## [1] 0.6785976

abs(cor(indicator.h(us), us^(0.12)))

## [1] 0.6782822

abs(cor(indicator.h(us), us^(0.15)))

## [1] 0.6786616

abs(cor(indicator.h(us), us^(0.16))) # chosen

## [1] 0.6786662

abs(cor(indicator.h(us), us^(0.17)))

## [1] 0.6786137

abs(cor(indicator.h(us), us^(0.18)))

## [1] 0.678506

Thus, k = 0.16 will be used as the reasonable value of k from now on (k.h).

k.h <- 0.16

ii) Sampling \(n = 10^6\) uniform random numbers, \(U_i \sim \text{Unif}(0,1)\), (i.i.d) \(i = 1,\dots, n\), estimate \(Cov[I(h(U) \le -19.1), U^k]\) with your \(k\) obtained in i) to find a suitable value of \(c\), and then estimate \(\mathbb P(I(h(U) \le -19.1)\) using both the control variable estimator and the naive average estimator. Give your estimates and the estimated variances of the two estimators. Find the estimated percentage decrease of variance reduction.

set.seed(stu_id)
us <- runif(n)

# 1: cov
cov(indicator.h(us), us^k.h) #  -0.03

## [1] -0.03002794

# 2 : c.h
c.h <- -cov(indicator.h(us), us^k.h)/var((us)^k.h) # 2.0865
c.h

## [1] 2.086522

ctrl_v_h <- function(us, k, c){
  I <- rep(0, length(us))
  I[sin(20*pi*us)+20*sin(pi*(us-1)/2) <= -19.1] <- 1
  
  Z <- I + c*(us^k - 1/(k+1))
  return(Z)
}

gs <- ctrl_v_h(us, k.h, c.h)
sum(gs)/n # 0.16276

## [1] 0.1627605

var(gs) # 0.073

## [1] 0.07337647

mean(indicator.h(us)) # 0.1624

## [1] 0.162406

var(indicator.h(us)) # 0.136

## [1] 0.1360304

(var(indicator.h(us)) - var(gs))/var(indicator.h(us)) # 46.06%

## [1] 0.4605878