Extension version homework4.2

ext_homework.jl

# Author: Sen Mu
# Extension of Homework 4.2 27-Sept-2018

# to install Distributions.jl, please use "julia>Pkg.add("Distributions")"
using Compat, Random, Distributions

# Generate seed of random numbers
Random.seed!(241)

#(1) sample from uniform random variable (a, b)
function sample_uniform(sample_size, a, b)
    println("parameter a=$(a), b=$(b)")
    deviates = []
    #create uniform random variable with a,b
    d_generate = Uniform(a, b)
    for i=1:11
        # generate deviates
        x = rand(d_generate, sample_size)
        x_mean = mean(x)
        push!(deviates, x_mean)
    end
    E_x = mean(deviates)
    Var_x = var(deviates)
    sigma = sqrt(Var_x)
    Ur = 2.23*(sigma/sqrt(11))
    #compute confidence interval
    lower_bound = E_x-Ur
    upper_bound = E_x+Ur
    println("Sample size = $(sample_size)")
    println("Mean:$(E_x), Deviation:$(sigma)")
    println("Confidence Interval 95%:($(lower_bound), $(upper_bound))\n")
end

#(2) using inverse transform sampling to generate Exponential random variable
function inverse2exp(lambda, sample_size)
    println("parameter lambda=$(lambda)")
    # first we should use a uniform distribution
    deviates = []
    d_generate = Uniform()
    for i=1:11
        # sampling from (0,1) as the value of Pareto distribution's cdf
        x = rand(d_generate, sample_size)
        # use inverse function to calculate points in Pareto distribution
        x1 = -(1/lambda)*log.(x)
        x_mean = mean(x1)
        # compute mean
        push!(deviates, x_mean)
    end
    E_x = mean(deviates)
    Var_x = var(deviates)
    sigma = sqrt(Var_x)
    Ur = 2.23*(sigma/sqrt(11))
    #compute confidence interval
    lower_bound = E_x-Ur
    upper_bound = E_x+Ur
    println("Sample size = $(sample_size)")
    println("Mean:$(E_x), Deviation:$(sigma)")
    println("Confidence Interval 95%:($(lower_bound), $(upper_bound))\n")
end

#(3) using inverse transform sampling to generate Pareto random variable
function inverse2pareto(alpha, xm, sample_size) 
    println("parameter alpha=$(alpha), xm=$(xm)")
    # first we should use a uniform distribution
    d_generate = Uniform()
    deviates = []
    for i=1:11
        # sampling from (0,1) as the value of Pareto distribution's cdf
        x = rand(d_generate, sample_size)
        # use inverse function to calculate points in Pareto distribution
        x1 = xm*(x.^(-1)).^(1/alpha)
        x_mean = mean(x1)
        # compute mean
        push!(deviates, x_mean)
    end
    E_x = mean(deviates)
    Var_x = var(deviates)
    sigma = sqrt(Var_x)
    Ur = 2.23*(sigma/sqrt(11))
    #compute confidence interval
    lower_bound = E_x-Ur
    upper_bound = E_x+Ur
    println("Sample size = $(sample_size)")
    println("Mean:$(E_x), Deviation:$(sigma)")
    println("Confidence Interval 95%:($(lower_bound), $(upper_bound))\n")
end

# set different sample sizes
sample_size = [10, 100, 1000, 10000, 100000, 1000000]

#main function of simulation
function main()
    println("Uniform distribution\n")
    for n in sample_size
        sample_uniform(n, 1,10)
    end
    for n in sample_size
        sample_uniform(n, 10, 20)
    end

    println("Exponential distribution\n")
    #try with different parameter lambda 
    for n in sample_size
        inverse2exp(3, n)
    end
    for n in sample_size
        inverse2exp(5, n)
    end
    for n in sample_size
        inverse2exp(12, n)
    end
    for n in sample_size
        inverse2exp(0.9, n)
    end

    println("Pareto distribution\n")
    for n in sample_size
        inverse2pareto(3,1, n)
    end
    #try when 1<alpha<2
    for n in sample_size
        inverse2pareto(1.5, 3, n)
    end
    # try when alpha is smaller than 1
    for n in sample_size
        inverse2pareto(0.9,3, n)
    end
end

main()