Ctrl+K

Hints and Extra Materials for Lab 1: Sampling Random Variables

Hints and Extra Materials for Lab 1: Sampling Random Variables#

For the class on Wednesday, January 10th

See also

Go back to Lab 1

A. Inverse transform sampling#

Hints for Step 1

  • What is the inverse function of the CDF of an exponential distribution? This is known as the inverse cumulative distribution function (also known as the quantile function). You can find it on Wikipedia.

  • The inverse transform sampling method says that applying the inverse CDF of a random variable \(X\) on a stanfard uniform random variable \(U \sim \text{Unif}(0,1)\) gives you the random variable \(X\).

  • Here’s a snippet to get you started with the code structure:

    def exponential_rvs(l=1, size=None):
    # here l is the parameter lambda
    rng = np.random.default_rng()
    uniform_rvs = rng.random(size)
    # Put your implementation for Part A, Step 1 here!
Hints for Step 2

  • To make a two-sample quantile-quantile plot, you will first need to calculate a given set of quantiles on each sample. Here’s a function that calculate the quantiles for two samples.

    def qqplot_data(x, y):
    n = min(len(x), len(y))
    q = (np.arange(n) + 0.5) / n
    qx = np.quantile(x, q)
    qy = np.quantile(y, q)
    return qx, qy

  • Then, you will need to make a scatter plot between the calculated quantiles from the two samples. You will also need to add a \(y=x\) line to your plot to guide your eye. You can then try to interpret the q-q plot you make!

B. Rejection sampling#

Hints for Step 1

  • The circle is fully enclosed within the square whose vertices are (1,1), (-1,1), (-1,-1), and (1,-1). Can you generate points that are sampled uniformly at random within this square using the built-in standard uniform random number generator?

  • Then, you can check whether each point that you generate is within the circle. Reject the ones that are not.

  • Here’s a snippet to get you started with the code structure:

    def uniform_rvs_in_circle(size=None):
    rng = np.random.default_rng()

    # Here we are using np.ravel to make sure we get 1D arrays
    x = np.ravel(rng.random(size=size))
    y = np.ravel(rng.random(size=size))

    # Now x, y are both between 0 and 1.
    # How do you translate that to be within the square of (1,1), (-1,1), (-1,-1), and (1,-1)?
    # Add your implementation here!

    # Then, check if the points are within the circle:
    mask_within = # Add your implementation here!

    # Finally, let's return the points that are within the circle
    return x[mask_within], y[mask_within]
Hints for Step 2

Note that you need 1,000 points in the circle! If you start with generating 1,000 points in the square, how many will be left after the rejection sampling method?

Hints for Step 3

If you think in the polar coordinates, what would you expect the distributions of \(r\) and \(\theta\) to be if the points are sampled uniformly at random within the circle?

C. Uniformly sampling random points on a 2D sphere#

Hints

If you are not sure where to start, I would suggest starting with Google and see if you can find existing implementation that applies to this question! You can use the implementation you found. Just make sure that (1) you verify it does what you think it does and (2) cite your source.

Extra materials (Not hints! Open only after you finish your function)

Take a look at the following three methods and then answer the questions in the lab.

import numpy as np


def sphere_surface_rvs_method1(size=None):
    rng = np.random.default_rng()
    phi = rng.random(size=size) * 2 * np.pi
    z = rng.random(size=size) * 2 - 1
    s = np.sqrt(1 - z * z)
    x = s * np.cos(phi)
    y = s * np.sin(phi)
    return np.atleast_2d(np.stack([x, y, z]))


def sphere_surface_rvs_method2(size=None):
    rng = np.random.default_rng()
    new_size = (3, int(size or 1))
    xyz = rng.random(size=new_size) * 2 - 1
    within_mask = (xyz * xyz).sum(axis=0) < 1
    xyz = xyz[:,within_mask]
    xyz /= np.sqrt((xyz * xyz).sum(axis=0))
    return xyz


def sphere_surface_rvs_method3(size=None):
    rng = np.random.default_rng()
    try:
        size = tuple(size)
    except TypeError:
        size = (1,) if size is None else (size,)
    new_size = (3,) + size
    xyz = rng.standard_normal(size=new_size)
    xyz /= np.sqrt((xyz * xyz).sum(axis=0))
    return xyz
Contents