Simulation

Monte Carlo Integration

Suppose that we want to perform the integration of the function f(x) = x^2 exp(-2x) in the interval (0,6)

my_fun = function(x){ # declaration of the function

(x^2)*exp(-2*x) # body of the function

}

# It is always a good idea to plot the function

par(mfrow =c(1,1)) # figure plotting windows

curve(my_fun(x), 0, 6, col = "red", lwd=2)

# Other way to plot it (Mainly for Python user)

x = seq(0, 6, length.out = 100)

print(x)

val_fun = my_fun(x) # evaluate the function at x values (mesh)

plot(x, val_fun, col = "red", type = "l", lwd=2)

# Now let us consider doing the integration

# Use of inbuilt function.

val_integral = integrate(my_fun, lower = 0, upper = 6)

# Note: The function integrate is used to perform

# integration in R

print(val_integral) # the value can be considered exact

# let us do the same integral by using simulation

# We shall write down the integral as an expectation of

# some random variable

# So, we can consider X as a uniformly distributed

# random variable over the interval (0,6). Its density function is

# f(x) = 1/6, 0<x<6

# So we are able write down the integral as an expectation of th

# random variable X^2e^(-2X), where X follows U(0,1)

#

# So the problems boils down to approximating the expectation

# of a random variable g(X) = X^2e^(-2X), where X follows

# U(0,1)

# Weak Law of Large Numbers

# If X_1, X_2, ... X_n is a sample of size n from some population

# with mean value mu. Then sample mean (X_1+X_2+...+X_n)/n

# converges to population mean mu in probability.

# Note 1: The sample mean is a random variable.

# But as n tends to infinity, this random quantity acts like a

# constant. The constant is nothing but the population mean

# (mu).

# In our problem, we have random variable g(X) and we

# want to evaluate E(g(X)). So if, we simulate a large number

# of samples from the distribution of g(X) and average them

# then the average value should be close to the E(g(X)), by

# Weak Law of Large Numbers.

# However, the distribution of g(X) may not be known. So we

# will simulate random numbers from the distribution of X

# (uniform(0,6)). And will evaluate the function g(.) at

# those values

n = 1000 # number of observations to be simulated from U(0,6)

x_sim_val = runif(n, min = 0, max = 6)

# It is always a good idea to check the codes immediately

# after you write it.

print(x_sim_val)

hist(x_sim_val, probability = TRUE, col = "grey")

g = my_fun

g_sim_val = g(x_sim_val)

g_sim_val

hist(g_sim_val, probability = TRUE, col = "grey")

approx_integral = 6*sum(g_sim_val)/n

approx_integral

# What we have done: We have approximate the integral

# by using the Monte Carlo integration method. It is very

# important to mention that based on how many number of

# simuated observations, you have computed the value of the

# integral. In this case, we have n=1000

# As a programmer, you should always explore and try things

# Had we consider n= 50, would we get the same accuracy?

# It is important to see the accuracy of the approximation

# as a function of the sample size n

# We want to evaluate the integral as function of n as well

# We want to do the same exercise for different values of n

# We can use an array for values of n and run a loop

n_vals = 1:1000 # an array containing 1,2,3,...,1000

print(n_vals)

approx_integral = numeric(length = length(n_vals))

for (n in n_vals) {

x_sim_val = runif(n, min = 0, max = 6)

g_sim_val = g(x_sim_val)

approx_integral[n] = 6*sum(g_sim_val)/n

}

print(approx_integral) # important

# We will write the integral as a function of n

plot(n_vals, approx_integral, type = "l", col = "red", lwd=2)

abline(h = val_integral$value, col = "grey", lwd=2)

# How to write a legend nicely

legend("topright", legend = c("monte carlo integral", "using integrate()"),

col = c("red", "grey"), lwd=c(2,2))

# We can see that as n increases the accuracy of the monte

# carlo integration increases.

# What we have learnt is that any definite integral

# can be evaluated using Monte Carlo integration. The key

# step is that we need to write the integral as an expectation

# and use the property that the sample mean converges to the

# population mean (in probability ) as n increases.

# Exercise:

# (a) Integrate x^2e^(-2x) in the interval (0, infty)

# (b) Let us suppose that I_n is the value of the integral

# based on a sample of size n. Thus, as the sample values

# changes, the value of I_n also changes. Essentially,

# the quantity I_n is a random variable, not a fixed number.

# Show that the expectation of I_n is I (the exact Integral)

# and Compute the variance of I_n using simulation. Visualize

# the shape of the SAMPLING DISTRIBUTION of I_n using

# 1000 replicatoin. Do this for n=10, n= 30, n= 50, n=100

# in the same graph.

Statistical Simulation: Lecture - I (Date: 12th January 2018)

There are n persons and their hats are numbered as 1, 2, 3,...,n. Suppose they throw the hat in the center of the room and asked to randomly pick one hat. We are interested to compute the probability that at least one person will get his own hat. This problem is discussed last semester in the class and appeared in examination too. The probability of at least one correct assignment of a hat is 1 - (1-1/n)^n. The limit as n tends to infinity is 1-1/e. Let us verify this claim using Statistical Simulation. We are planning to do: we randomly select n numbers from the set {1, 2, ..., n} by simple random sampling without replacement, so that all numbers will be selected in every experiment. This experiment can be performed using R by using the sample(1:n, size = n, replace = FALSE). The replace = FALSE indicates that once a number is selected it will not appear again. Use set.seed(123) function so that the same answer can be obtained for all.

n = 10 # Number of hats

x = 1:n # Reference number to checked for each experiment

m = 1000 # Number of experiments

out = matrix(data = NA, ncol = n, nrow = m) # Experimental outcomes to be stored in out

for(i in 1:m){

out[i, ] = sample(x, size = length(x), replace = FALSE)

}

s = numeric(m) # Storing the number of success

for(i in 1:m){

if(sum(x==out[i,])>0)

s[i] = 1

}

prob = sum(s)/m

print(prob) # Probability of at least one success

Exercise: Write an R program to perform the above experiment with varying n, the number of hats. Plot the probability of at least one matching (obtained from simulation) with respect to n.

Solution:

n = 1:50 # number of persons

m = 10000 # number of times experiment to be performed to estimate the probability of at least on matching

prob_exact = 1-(1-1/n)^n

prob_sim = numeric(length = length(n))

for(i in 1:length(n)){

x = 1:n[i]

out = matrix(data = NA, ncol = n[i], nrow = m)

for(j in 1:m){

out[j, ] = sample(x, size = length(x), replace = FALSE)

}

s = numeric(m) # Storing the number of success

for(j in 1:m){

if(sum(x==out[j,])>0)

s[j] = 1

}

prob_sim[i] = sum(s)/m

}

plot(n, prob_sim, type = "l", col = "red", ylim = c(0.4,1), lwd=2)

lines(prob_exact, type= "l", col = "blue", lwd=2)

legend("topright", legend = c("empirical prob.", "exact prob."), lwd=c(2,2), col=c("red","blue"), bty="n")

While solving the above problem, I felt that a short introduction to R is required, hence the following codes have been developed during lecture and the basic things were explained. These will be rigorously covered in your Software Lab-II course. This page is only for statistical simulations. Find your valuable time to run these codes and understand their functionalities for a better understanding of R Programming. We shall explicitly write programs, so it is better you get familiar with vector, matrix, data.frame, list in R.

In the following we evaluate the probability that X lies between 0 to 1. The distribution of X is normal with mean 1 and sigma 1. The required area under the curve is done by means of simulation.

mu = 1

sigma = 1

fun = function(x){

(1/(sqrt(2*pi)*sigma))*exp(-(x-mu)^2/(2*sigma^2))

}

prob = integrate(fun, 0, 1)$value

curve(fun(x), -3, 5, col = "red", lwd=3)

segments(0,0,0,fun(0), lwd=2)

abline(h=0, lwd=2)

segments(1,0,1,fun(1), lwd=2)

rep = 100

count = numeric(rep)

for(i in 1:rep){

x = rnorm(i, mean = mu, sd = sqrt(sigma))

count[i] = sum(x>0 & x<1)/i

}

plot(1:rep, count, col = "red", pch = ".")

abline(h = prob, col = "green", lwd=2)

Lecture - III (Date: 17th January 2018)

We start our discussion with some plot options. Few students had doubts about some plot options. First, we clarify those, in particular, type = option and how to put legends in a plot. The code which is developed during the lecture follows: Concepts are self-explanatory. Use help(plot), help(legend), help(lines), help(points), help(seq) to know more about those functions.

x = seq(1, 10 ,length.out = 100)

print(x)

y = sin(x)

z = cos(x)

plot(x, y, col = "red" )

points(x, z, col = "blue")

plot(x, y , col = "red", type="l")

lines(x, z, col = "blue", type = "l")

plot(x, y , col = "red", type="b", lwd =2)

lines(x, z, col = "blue", type = "b" ,lwd=2)

# Adding legend

plot(x, y , col = "red", type="l", lwd =2, xlim = c(1, 10),

ylim=c(-1,1.5))

lines(x, z, col = "blue", type = "l", lwd=2)

legend("topright", legend = c("sin(x)","cos(x)"),

col = c("red", "blue"), lwd =c(2,2), bty = "n")

Then we have recalled what we have done in the previous lecture. In the previous lecture, we have written a function generate_exp() to generate n random numbers from the exponential distribution with a fixed rate parameter lambda. We do the process for the logistic(mu, beta) distribution. Note that for the logistic distribution mu is the location parameter and beta is the scale parameter. It is important to introduce the concept of location and scale family.

Location Family: Let f(x) be any probability density function. Then the family of probability density functions f(x-mu), indexed by the parameter mu, -inf < mu < inf, is called location family with standard pdf f(x) and mu is called the location parameter for the family.

Scale Family: Let f(x) be any probability density function. The for any sigma>0, the family of probability density function (1/sigma)f(x/sigma), indexed by the parameter sigma, is called the scale family with standard pdf f(x) and sigma is called the scale parameter of the family.

Location-Scale Family: Let f(x) be any probability density function. The for any mu, -inf < mu < inf, and any sigma>0, the family of probability density functions (1/sigma)f((x-mu)/sigma), indexed by the parameter (mu, sigma), is called the location-scale family with standard pdf f(x); mu is called the location parameter and sigma is called the scale parameter.

Recall that in the first semester we have studied the logistic distribution with parameter mu and beta. The logistic density function is given by f(x) = (1/beta)*(exp(-(x-mu)/beta))/(1 + exp(-(x-mu)/beta))^2, -inf < mu < inf, beta > 0. Note that the logistic family is a location-scale family with the location parameter mu and scale parameter beta. Now we shall write our own function to generate random numbers from the logistic distribution with parameter mu and sigma. Since the distribution function can be explicitly computed from the logistic density function, it is easy to verify the probability integral transform method. Explicit computation (done during the lecture) provides the logistic distribution function as F(x) = 1/(1 + exp(-(x-mu)/beta)). So if u is a uniform random number from [0,1] then the inverse of F at u will give a logistic random variable. Inv(F)(u) = mu + beta * log(u/(1-u)). The following function generate_log(n, mu, beta) utilizes this fact to generate random numbers from logistic(mu, beta) distribution. We also compare our results with the generator rlogis() in R by means of histograms. The corresponding probability density function is overlayed on the histogram for verification.

generate_log = function(n=10, mu=1, beta=2){

n = n # number of random values

mu = mu # location parameter

beta = beta # scale parameter

if(beta <= 0){

return("Invalid scale parameter")

}

u = runif(n, min = 0, max = 1)

x = mu + beta * log(u/(1-u))

return(x)

}

par(mfrow = c(1,2))

mu = 5

beta =3

x = generate_log(n=10000, mu=5, beta = 3)

hist(x, probability = TRUE, ylim = c(0,0.08),

main = "Generate using Inv tranform")

logistic_pdf = function(x){

(1/beta)*(exp(-(x-mu)/beta))/(1+exp(-(x-mu)/beta))^2

}

curve(logistic_pdf(x), -20, 30, lwd=2,

col = "red", add = TRUE)

y = rlogis(n = 10000, location = 5, scale = 3)

hist(y, probability = TRUE, main = "Generate using R")

# Checking for probability density function

exp_pdf = function(x){

exp(-x)

}

integrate(exp_pdf, 0, Inf)

f= function(x){

x^2

}

integrate(f, -1, 1)

Convergence in Probability: This is a very important concept in statistics and plays a very crucial role in statistical simulation. In the lectures, we have already checked how to show the convergence in probability (recall the use of Chebychev's inequality). We shall develop R codes to visualize this concept. There is one magical theorem, we studied in the class, is called the Weak Law of Large Numbers (WLLN). WLLN states that, suppose we generate a random sample of size n from a population distribution which has a finite mean (mu, say) and finite variance, then as n increases the sample mean converges to the population mean (mu) in probability. In layman language, it means that for large n, the sample mean values will be highly concentrated about mu. In the following code, we visualize this fact by simulating random numbers from the uniform distribution in the interval [a,b]. The mean (mu) = (a+b)/2 and var = (b-a)^2/12. So the population distribution has finite mean and finite variance, hence WLLN can be applied. For each n we generate n uniform random numbers from [a,b] and compute the sample mean bar(x_n). Plot n versus bar(x_n) for a range of n values.

par(mfrow= c(1,1))

a = 1; b = 10

n = 10000 # sample size

x_mean = numeric(n)

for(i in 1:n){

x = runif(i, min = a, max = b)

x_mean[i] = mean(x)

}

plot(1:n, x_mean, type= "p", col = "red")

abline(h = (a+b)/2, col = "blue") # mean values of are concentrated about this line

We can visualize the same thing in terms of distribution as well. Basically, the random variables bar(X_n) will be asymptotically degenerate and mu. So here, we plot the distribution of the sample mean for different sample sizes. Please understand the logic of each loop. First, visualize them and draw them on a piece of paper and then start writing the code. Understand the graphical output clearly and match with the lecture notes (theory class).