Syllabus

Syllabus

Resources

Week 1

Sampling Functions

Sampling Functions

Julia provides several functions for creating samples, depending on the type of sampling you need to perform. Here’s a breakdown of common functions and their uses:

1. rand():

• Purpose: Generates random numbers. Can be used for simple random sampling with replacement.
• Usage:
  • rand(): Generates a random floating-point number between 0 and 1.
  • rand(n): Generates an array of n random floating-point numbers between 0 and 1.
  • rand(dims...): Generates an array with the specified dimensions filled with random floating-point numbers between 0 and 1.
  • rand(rng, n): Generates n random numbers using the random number generator rng.
  • rand(1:10, n): Generates n random integers between 1 and 10 (inclusive). This is useful for simulating dice rolls or drawing numbers from a hat with replacement.

2. sample():

• Purpose: Provides more general sampling functionality, including with and without replacement, and weighted sampling.
• Usage:
  • sample(population, n): Samples n elements from population without replacement.
  • sample(population, n, replace=true): Samples n elements from population with replacement.
  • sample(population, Weights(weights), n): Samples n elements from population with probabilities proportional to the given weights. weights should be a vector of the same length as population.
  • sample(rng, population, n): Samples using the random number generator rng.

3. shuffle():

• Purpose: Randomly shuffles the elements of an array. This can be useful for creating a random permutation of your data, which can then be used for sampling or splitting into training/testing sets.
• Usage:
  • shuffle(arr): Returns a new array with the elements of arr randomly shuffled.
  • shuffle!(arr): Shuffles the elements of arr in place.

4. randperm():

• Purpose: Generates a random permutation of integers from 1 to n. Useful for creating random indices for sampling.
• Usage:
  • randperm(n): Returns a random permutation of the integers 1 to n.

5. StatsBase functions:

The StatsBase package provides more advanced sampling methods, including:

• bootstrap: For bootstrap resampling.
• stratifiedsample: For stratified sampling.

You’ll need to install the package using ] add StatsBase in the Julia REPL.

Choosing the right function:

• Simple random sampling with replacement: Use rand() with a range for integers (e.g., rand(1:10, n)) or sample(population, n, replace=true).
• Simple random sampling without replacement: Use sample(population, n).
• Weighted sampling: Use sample(population, Weights(weights), n).
• Creating random permutations: Use randperm() or shuffle().
• Advanced sampling techniques (bootstrap, stratified, etc.): Use functions from the StatsBase package.

Example: Sampling with replacement

population = ['A', 'B', 'C', 'D', 'E']
n_samples = 10

samples_with_replacement = sample(population, n_samples, replace=true)
println(samples_with_replacement)

Example: Sampling without replacement

population = ['A', 'B', 'C', 'D', 'E']
n_samples = 3

samples_without_replacement = sample(population, n_samples)
println(samples_without_replacement)

Example: Weighted sampling

population = ['A', 'B', 'C']
weights = [0.2, 0.5, 0.3]  # Probabilities for each element
n_samples = 10

weighted_samples = sample(population, Weights(weights), n_samples)
println(weighted_samples)

Remember to install StatsBase if you need its functions: ] add StatsBase in the Julia REPL. Then, using StatsBase in your code.

Simulating Exponential Random Variates using Julia

Simulating Exponential Random Variates

This demonstration provides a more comprehensive and practical example of how to simulate exponential random numbers in Julia.

It covers both common methods, includes verification steps, and offers optional visualization to help you understand and trust the results. Remember to install the Plots package if you intend to use the histogram functionality.

# Method 1: Using the built-in `rand` function with the Exponential distribution

# Specify the rate parameter (lambda).  The mean of the Exponential distribution is 1/lambda.
lambda = 2.0  # Example rate

# Generate a single exponential random number
exponential_rn = rand(Exponential(1/lambda))
println("Single Exponential random number (Method 1): ", exponential_rn)

# Generate an array of exponential random numbers
num_samples = 1000
exponential_rns = rand(Exponential(1/lambda), num_samples)

# Verify the mean and variance (should be approximately 1/lambda and 1/lambda^2 respectively)
mean_simulated = mean(exponential_rns)
var_simulated = var(exponential_rns)
println("Mean of simulated (Method 1): ", mean_simulated)
println("Variance of simulated (Method 1): ", var_simulated)

# Method 2: Using the inverse transform sampling method

# Generate a uniform random number between 0 and 1
uniform_rn = rand()

# Apply the inverse of the CDF of the exponential distribution: -log(1-U)/lambda
exponential_rn_inverse = -log(1 - uniform_rn) / lambda
println("Single Exponential random number (Method 2): ", exponential_rn_inverse)


# Generate an array of exponential random numbers using inverse transform sampling
uniform_rns = rand(num_samples)
exponential_rns_inverse = -log.(1 .- uniform_rns) ./ lambda  # Use . for element-wise operations

# Verify the mean and variance
mean_simulated_inverse = mean(exponential_rns_inverse)
var_simulated_inverse = var(exponential_rns_inverse)
println("Mean of simulated (Method 2): ", mean_simulated_inverse)
println("Variance of simulated (Method 2): ", var_simulated_inverse)

Visualization the Distribution

# Plotting a histogram to visualize the distribution (optional, requires Plots.jl)
using Plots

histogram(exponential_rns, bins = 50,  xlabel = "Random Number", ylabel = "Frequency",
          title = "Histogram of Exponential Random Numbers (Method 1)",
          normalize = true, label = "Simulated") # Normalize to probability density

# Overlay the theoretical PDF for comparison
x = range(0, 5/lambda, length=100) # Adjust range as needed
plot!(x, lambda * exp.(-lambda*x), label = "Theoretical PDF", linewidth=2)

savefig("exponential_histogram_method1.png")


histogram(exponential_rns_inverse, bins = 50,  xlabel = "Random Number", ylabel = "Frequency",
          title = "Histogram of Exponential Random Numbers (Method 2)",
          normalize = true, label = "Simulated") # Normalize to probability density

# Overlay the theoretical PDF for comparison
x = range(0, 5/lambda, length=100) # Adjust range as needed
plot!(x, lambda * exp.(-lambda*x), label = "Theoretical PDF", linewidth=2)

savefig("exponential_histogram_method2.png")

Explanation and Key Improvements:

1. Two Methods: The code demonstrates two common ways to generate exponential random numbers:

   • Method 1: Using the built-in rand(Exponential(1/lambda)) function. This is the most straightforward and efficient way in Julia. Note that Exponential(θ) in Julia uses the mean (θ) parameterization, which is the inverse of the rate (λ) often used in other contexts.
   • Method 2: Using the inverse transform sampling method. This method is more general and can be applied to other distributions as well. It involves generating a uniform random number and then transforming it using the inverse of the cumulative distribution function (CDF) of the exponential distribution.

2. Clearer Parameterization: The code emphasizes the use of the rate parameter (lambda) and clearly explains the relationship between the rate and the mean of the exponential distribution (mean = 1/lambda).

3. Verification: The code calculates and prints the mean and variance of the generated random numbers to verify that they are close to the expected theoretical values (1/lambda and 1/lambda^2, respectively). This is a good practice to ensure the simulation is working correctly.

4. Histograms (Optional): The code includes an optional section that uses the Plots.jl package to create histograms of the generated random numbers and overlay the theoretical probability density function (PDF). This allows you to visually compare the simulated distribution to the expected distribution. You’ll need to install the Plots package if you want to use this part: ] add Plots in the Julia REPL. The savefig function is used to save the plot as a PNG image.

5. Element-wise operations: The code uses the dot operator (.) for element-wise operations when generating arrays of random numbers using the inverse transform method. This is crucial for efficiency in Julia.

6. Comments and Clarity: The code is well-commented to explain each step, making it easier to understand and modify.

Week 2

Bootstrap Sampling

Bootstrap Sampling

Bootstrap sampling is a powerful statistical technique for estimating the distribution of a sample statistic by resampling with replacement. It’s especially useful when the theoretical distribution of the statistic is complex or unknown. Julia, a high-performance programming language, is well-suited for implementing this method. Let’s dive into a step-by-step tutorial on how to perform bootstrap sampling in Julia.

Step-by-Step Tutorial: Bootstrap Sampling with Julia

using Pkg
Pkg.add("Distributions")
Pkg.add("Statistics")

using Distributions
using Statistics

# Generate a sample dataset of 100 observations from a normal distribution
n = 100
sample_data = rand(Normal(0, 1), n)

function bootstrap_samples(data, n_samples)
    n = length(data)
    bootstrap_samples = zeros(n_samples, n)
    for i in 1:n_samples
        bootstrap_samples[i, :] = sample(data, n, replace=true)
    end
    return bootstrap_samples
end

n_bootstrap_samples = 1000
bootstrap_data = bootstrap_samples(sample_data, n_bootstrap_samples)

bootstrap_means = mean(bootstrap_data, dims=2)

# Bootstrap estimate of the mean
bootstrap_mean_estimate = mean(bootstrap_means)

# Standard error of the bootstrap estimate
bootstrap_se = std(bootstrap_means)

# Construct a 95% confidence interval
confidence_interval = quantile(bootstrap_means, [0.025, 0.975])

println("Bootstrap Mean Estimate: $bootstrap_mean_estimate")
println("Bootstrap Standard Error: $bootstrap_se")
println("95% Confidence Interval: $confidence_interval")

Summary

In this tutorial, we have covered the basics of bootstrap sampling in Julia. We generated a sample dataset, defined a bootstrap sampling function, created bootstrap samples, calculated the statistic of interest, and analyzed the results. Bootstrap sampling is a versatile technique that can be applied to various statistical problems, and Julia makes it easy to implement and analyze.

Feel free to experiment with different datasets and statistics. Happy coding!

Is there anything specific you’d like to explore further in this tutorial?

Week 5

• Monte Carlo Simulation in Julia (Week 5 Lesson Plan)

https://rpubs.com/JuliaWorkshop/Markov_Chains_Julia_Tutorial https://rpubs.com/JuliaWorkshop/Markov_Chain_Example_Julia

Markov Chains with Julia

Several Julia packages are well-suited for modeling Markov chains, each with its strengths and focuses. Here’s a breakdown of the most commonly used options:

1. MarkovChains.jl: This package is specifically designed for working with discrete-time Markov chains (DTMCs). It provides a comprehensive set of tools for:

• Defining Markov chains: You can create Markov chains by specifying the transition probability matrix.
• Analyzing Markov chains: MarkovChains.jl offers functions for calculating:
  • Stationary distributions (the long-run probabilities of being in each state).
  • Transient probabilities (probabilities of being in a state at a specific time).
  • Mean first passage times (expected time to reach a state).
  • Recurrence and transience properties of states.
• Simulating Markov chains: You can generate random sequences of states according to the transition probabilities.
• Higher-order Markov chains: The package also supports higher-order Markov chains where the next state depends on a fixed number of previous states.

using MarkovChains

# Define the transition probability matrix
P = [0.7 0.3; 0.2 0.8]  # Example 2-state Markov chain

# Create a Markov chain object
mc = MarkovChain(P)

# Calculate the stationary distribution
π = stationary_distribution(mc)
println("Stationary Distribution: ", π)

# Simulate the Markov chain for 10 steps
states = simulate(mc, 10)
println("Simulated States: ", states)

# Calculate the probability of being in state 1 after 5 steps, starting in state 1
p = transient_probabilities(mc, 5, 1)
println("Transient probability: ", p[1])

# Calculate the mean first passage time from state 1 to state 2
m = mean_first_passage_time(mc, 1, 2)
println("Mean First Passage Time: ", m)

2. StatsBase.jl: While not exclusively for Markov chains, StatsBase.jl provides some useful functions that can be applied to Markov chains, especially when dealing with discrete data and probabilities. You can use it in conjunction with other packages.

3. Distributions.jl: This package is essential for working with probability distributions, which are fundamental to Markov chains. You’ll often use Distributions.jl to define the underlying probability distributions that govern transitions between states (especially for continuous-time Markov chains or more complex models).

4. QuantEcon.jl: As mentioned before, QuantEcon.jl is a powerful package for quantitative economics, and it includes tools for working with Markov chains, particularly in the context of dynamic programming and other economic models.

5. SparseArrays.jl: If you’re dealing with very large Markov chains with sparse transition matrices (where most entries are zero), using SparseArrays.jl can be crucial for memory efficiency and performance. You can represent the transition matrix as a sparse array, which can significantly reduce storage requirements and speed up computations.

Choosing the Right Package:

• MarkovChains.jl: This is generally the best choice for most standard discrete-time Markov chain analysis. It’s specifically designed for this purpose and offers a complete set of tools.
• QuantEcon.jl: If you are working on more complex economic models that involve Markov chains, then QuantEcon.jl may be a more suitable choice, as it integrates well with other tools for dynamic programming and economic modeling.
• SparseArrays.jl: Use this when dealing with large, sparse Markov chains.
• StatsBase.jl and Distributions.jl: These are useful supporting packages, especially when you need to work with specific probability distributions or perform general statistical calculations related to your Markov chain analysis.

For most common Markov chain analysis tasks, MarkovChains.jl will be the most direct and convenient option. If you are working on more specialized applications (e.g., large-scale systems, connections to economic models), then you might want to explore the other packages mentioned above. Remember to install the chosen package(s) using the Julia package manager (e.g., ] add MarkovChains).

Week 10

Numerical Integration

Numerical Integration

Numerical integration is a technique for approximating the definite integral of a function when it is difficult or impossible to find an analytical solution. It involves using numerical methods to estimate the area under the curve of the function within the given limits of integration.

Here’s a breakdown of the key aspects of numerical integration:

What is a definite integral?

• A definite integral represents the net area between the graph of a function and the x-axis, within specified limits on the x-axis.
• It has numerous applications in various fields, including physics, engineering, economics, and statistics.

Why use numerical integration?

• Analytical solutions may not exist: For some functions, finding an antiderivative (the function whose derivative is the original function) in a closed form may be impossible.
• Complex functions: Even if an antiderivative exists, it might be too complex to evaluate easily.
• Data-driven functions: When the function is defined by a set of data points rather than an explicit formula, numerical integration is necessary.

Common numerical integration methods:

• Rectangle rule: Approximates the area under the curve using rectangles.
• Trapezoidal rule: Approximates the area using trapezoids, which generally provide a better approximation than rectangles.
• Simpson’s rule: Uses quadratic functions to approximate the curve, resulting in even higher accuracy.
• Gaussian quadrature: Selects specific points and weights to evaluate the function, optimizing accuracy for a given number of function evaluations.

Factors affecting accuracy:

• Choice of method: Different methods have varying levels of accuracy.
• Number of intervals/points: Increasing the number of intervals or points used in the approximation generally improves accuracy.
• Smoothness of the function: Smoother functions are typically easier to approximate accurately.

Applications of numerical integration:

• Calculating areas and volumes: Finding the area of irregular shapes or the volume of complex objects.
• Solving differential equations: Many differential equations cannot be solved analytically, and numerical integration is used to approximate their solutions.
• Probability and statistics: Calculating probabilities associated with continuous random variables.
• Financial modeling: Pricing financial instruments and evaluating risk.
• Computer graphics: Rendering images and simulating physical phenomena.

Advantages of numerical integration:

• Applicable to a wide range of functions: Can handle functions that are difficult or impossible to integrate analytically.
• Relatively easy to implement: Many numerical integration methods are straightforward to program.
• Controllable accuracy: The accuracy of the approximation can be adjusted by changing the number of intervals or points.

Limitations of numerical integration:

• Approximation, not exact: Numerical integration provides an approximation of the definite integral, not the exact value.
• Computational cost: Increasing the accuracy often requires more function evaluations, which can be computationally expensive.
• Potential for error: Numerical methods can introduce errors, such as round-off errors and truncation errors.

In summary, numerical integration is a powerful tool for approximating definite integrals when analytical solutions are not feasible. It has wide-ranging applications in various fields and provides a practical way to estimate areas, volumes, and solutions to complex problems.

Numerical Integration

Numerical integration by simulation refers to a specific type of numerical integration technique that leverages the power of simulation to estimate the value of a definite integral. It’s also known as Monte Carlo integration. This method is particularly useful when dealing with integrals that are difficult or impossible to solve analytically, or when the function being integrated is complex or defined by a set of data points.

Here’s how it works:

1. Random Sampling: Instead of dividing the integration interval into smaller subintervals like in traditional numerical integration methods, Monte Carlo integration relies on generating random samples within the integration domain.
2. Function Evaluation: The function to be integrated is evaluated at each of these randomly sampled points.
3. Averaging: The average of the function values obtained in step 2 is calculated.
4. Scaling: This average is then multiplied by the size of the integration domain to obtain an estimate of the definite integral.

The underlying principle behind Monte Carlo integration is the Law of Large Numbers, which states that the average of a large number of random samples will converge towards the true expected value. In this context, the expected value is the definite integral we are trying to estimate.

Advantages of Monte Carlo Integration:

• Handles complex integrals: It can be applied to integrals with complex integrands or integration domains, where traditional methods might struggle.
• Deals with high-dimensional integrals: It is particularly effective for high-dimensional integrals, where other numerical integration techniques become computationally expensive.
• Simple to implement: The basic algorithm is relatively straightforward to implement.

Limitations of Monte Carlo Integration:

• Accuracy: The accuracy of the estimate depends on the number of random samples used. Increasing the number of samples improves accuracy but also increases computational cost.
• Convergence: Convergence can be slow compared to other numerical integration methods, especially for low-dimensional integrals.
• Randomness: The results can vary slightly due to the random nature of the sampling process.

Applications of Monte Carlo Integration:

• Finance: Pricing complex financial instruments and options.
• Physics: Simulating particle transport and interactions.
• Engineering: Reliability analysis and risk assessment.
• Computer graphics: Rendering images and simulating light transport.

In summary, Monte Carlo integration is a powerful technique for estimating definite integrals, especially in situations where traditional methods are not feasible. It leverages the power of random sampling and the Law of Large Numbers to provide approximate solutions to complex integration problems.

Week 11

Box-Muller Transform

The algorithm used in the previous example is the Box-Muller transform, specifically the polar form of the Box-Muller transform. It’s a method for generating pseudo-random numbers from the standard normal distribution (mean 0 and variance 1) given a source of uniform random numbers. Here’s a detailed breakdown:

1. The Goal:

The objective is to transform two independent random variables, U1 and U2, which are uniformly distributed between 0 and 1, into two independent random variables, X and Y, which follow a standard normal distribution.

2. The Box-Muller Transform (Polar Form):

The polar form avoids the calculation of trigonometric functions, which can be computationally expensive. It proceeds as follows:

• Generate Two Uniform Random Numbers: Two independent random numbers, U1 and U2, are generated from a uniform distribution on the interval (0, 1). In the code, this is done using rand().

• Transform to [-1, 1]: The uniform random numbers are transformed to the interval [-1, 1]:

  • V1 = 2 * U1 - 1
  • V2 = 2 * U2 - 1

• Calculate s: Calculate the sum of squares:

  • s = V1² + V2²

• Rejection Sampling (The “Polar” Part): This is the key part that makes this the polar form. We only accept the pair (V1, V2) if it falls within the unit circle (i.e., s < 1). If s ≥ 1, we reject the pair and generate new U1 and U2 until we find a pair that satisfies s < 1. This is a form of rejection sampling. It ensures that the distribution of the accepted points is uniform within the unit circle. The condition s > 0 is also important to prevent issues with the log function if s is very close to zero.

• The Transformation (If Accepted): If the pair is accepted (s < 1), then the following transformation is applied:

  • X = sqrt(-2 * log(s) / s) * V1
  • Y = sqrt(-2 * log(s) / s) * V2

  These X and Y will be independent standard normal random variables.

3. Why it Works:

The Box-Muller transform is based on the fact that if X and Y are independent standard normal variables, then X² + Y² follows a chi-squared distribution with two degrees of freedom, which is equivalent to an exponential distribution with a mean of 2. The polar form is an optimized version of the Box-Muller method that avoids the use of trigonometric functions by using a rejection sampling step. The rejection sampling ensures that the distribution of the points is uniform within the unit circle. This uniformity is crucial for the transformation to produce standard normal variates.

4. In the Code:

The while loop and the if s < 1 condition implement the rejection sampling. The lines inside the if block perform the Box-Muller transformation to produce the normal variates. The code generates two normal random variables per loop iteration, which is why the loop counter i is incremented twice and the if i <= n check is in place.

5. Advantages of the Polar Method:

• No Trigonometric Functions: It avoids the use of sine and cosine, which can be computationally expensive.
• Simpler: Conceptually, it’s a bit simpler to understand than the original Box-Muller transform.

6. Disadvantages:

• Rejection Sampling: The rejection sampling step means that some pairs of uniform random numbers are discarded. This can make the polar method slightly less efficient than other normal random number generators (although it’s still quite fast). The probability of rejection is approximately (1 - π/4) ≈ 0.215. This means that, on average, you’ll need to generate about 4/π ≈ 1.27 times more uniform random number pairs than normal random numbers you want.

In summary, the polar form of the Box-Muller transform is a clever and efficient way to generate standard normal random variates using uniform random numbers. It’s widely used in Monte Carlo simulations and other applications that require normally distributed random numbers.

The Box-Muller transform is a method for generating pairs of independent, standard, normally distributed random numbers given a source of uniformly distributed random numbers. It’s a fundamental algorithm in computational statistics and Monte Carlo simulations.

Here’s a breakdown of the Box-Muller transform:

1. The Goal:

The objective is to transform two independent random variables, U1 and U2, which are uniformly distributed between 0 and 1, into two independent random variables, X and Y, which follow a standard normal distribution (mean 0 and variance 1).

2. The Transformation:

The Box-Muller transform involves the following steps:

• Generate Two Uniform Random Numbers: Two independent random numbers, U1 and U2, are generated from a uniform distribution on the interval (0, 1).

• Calculate Intermediate Values:

  • R = sqrt(-2 * ln(U1))
  • Θ = 2 * π * U2

• Calculate Normal Variates:

  • X = R * cos(Θ)
  • Y = R * sin(Θ)

These X and Y will be independent standard normal random variables.

3. Why it Works:

The Box-Muller transform is based on the fact that if X and Y are independent standard normal variables, then X² + Y² follows a chi-squared distribution with two degrees of freedom, which is equivalent to an exponential distribution. The transformation essentially maps points from a uniform distribution in polar coordinates (R and Θ) to Cartesian coordinates (X and Y), where the resulting X and Y are normally distributed.

4. In Simpler Terms:

Imagine you have a dartboard where you throw darts randomly, and they land uniformly across the board. The Box-Muller transform is like a recipe that takes where the dart landed (represented by two random numbers) and converts it into two new numbers that would have landed on a bullseye if the throws were normally distributed.

5. Key Points:

• Independence: The Box-Muller transform generates two independent standard normal random variables.
• Uniform to Normal: It transforms uniform random variables into normal random variables.
• No Approximations: The Box-Muller transform is an exact method (not an approximation).

6. Applications:

The Box-Muller transform is widely used in:

• Monte Carlo simulations: To generate random numbers from a normal distribution, which is often needed in these simulations.
• Statistics: For various statistical computations and simulations.
• Computer graphics: To create realistic-looking random effects.

7. Variations:

There is a variation called the polar form of the Box-Muller transform, which is computationally more efficient as it avoids the use of trigonometric functions (sine and cosine).

In Summary:

The Box-Muller transform is a fundamental and widely used algorithm for generating normally distributed random numbers. It provides an exact method for converting uniform random numbers into standard normal random variables, making it a valuable tool in various fields.

Week 12