<jupyter_start><jupyter_text>Homework 2**Due: 01/30/2018** Instructions+ In any case, develop the code and generate the figures you need to solve the problems using this notebook.+ For the answers that require a mathematical proof or derivation you can either:        - Type the answer using the built-in latex capabilities. In this case, simply export the notebook as a pdf and upload it on gradescope; or    - you can print the notebook (after you are done with all the code), write your answers by hand, scan, turn your response to a single pdf, and upload on gradescope.+ The total homework points are 100. Please note that the problems are not weighed equally.  Student details+ **First Name:**+ **Last Name:**+ **Email:**<jupyter_code># Here are some modules that you may need - please run this block of code:
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set_context('paper')
import numpy as np
import scipy.stats as st<jupyter_output><empty_output><jupyter_text>Problem 1 Part ALet $X$ be a continuous random variable with CDF:$$F(x) = p(X \le x).$$Show that the random variable$$Z = F(X)$$    is distributed uniformly in $[0,1]$.    Hint: Show that:$$    F_Z(z) := p(Z \le z) = z.$$**Proof:** *Type your proof here. Delete that ```` line (it just makes some white space).*  Part B The theorem you have just proved is very useful when you want to test if a set of observations $x_1,x_2,\dots,x_N$ has been indeed independent realizations of a random variable $X$ with CDF $F(x)$.Part A proves that, if this hypothesis is valid, then the transformed dataset $z_1,z_2,\dots,z_N$, where$$z_i = F(x_i),$$should be distributed uniformely in $[0,1]$.In other words, the empirical histrogram of $z_1,z_2,\dots,z_N$ should be a flat line.Use this observation to find the distribution from which this dataset was sampled:<jupyter_code>data = np.loadtxt('hw2_p1b_data.txt')
fig, ax = plt.subplots()
ax.hist(data, normed=True, alpha=0.5)
ax.set_xlabel('$x$')
ax.set_ylabel('Empirical PDF')
ax.set_title('Histogram of data')<jupyter_output><empty_output><jupyter_text>The correct distribution is one of the following:1. Standard normal, $\mathcal{N}(0,1)$;2. Normal with mean 2 and variance 2, $\mathcal{N}(2,2); or3. Exponential with rate parameter $1$, $\mathcal{E}(1)$; or4. Exponential with rate parameter $2$, $\mathcal{E}(2)$; or5. Exponential with rate parameter $10$, $\mathcal{E}(10)$; or6. Gamma distribution with parameters $\alpha=2.$ and $\beta=3.$.Systematically go over these distributions and try to determine which one generated the data.All the required CDF's and inverse CDF's are implemented in [scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html).Check also [scipy.stats.rv_continuous](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.htmlscipy.stats.rv_continuous).)Please pay special attention to the defintiion of the probability distributions of the various random variables and how you can control their parameters.As a hint, here is how you can test for $\mathcal{N}(0,1)$:<jupyter_code># Testing for N(0,1):
transformed_data = st.norm(loc=0, scale=1).cdf(data) # cdf() gives the CDF of a random variable
# If the data came from N(0,1), the histogram of the transformed_data should match that of a uniform:
fig, ax = plt.subplots()
ax.hist(transformed_data, normed=True, alpha=0.5)
ax.plot(np.linspace(0, 1,50), np.ones(50), lw=2, color='r') # This is the line that you should try to much.
ax.set_xlabel('$z = F(x)$')
ax.set_ylabel('Empirical PDF')
ax.set_title('Testing the $\mathcal{N}(0,1)$')<jupyter_output><empty_output><jupyter_text>Problem 2 Part ALet $U\sim \mathcal{U}(0,1)$ and let $F(x)$ be a function with the following properties:+ $F(-\infty) = 0$,+ $F(\infty)=1$, + and $x_1 \le x_2\Rightarrow F(x_1)\le F(x_2)$.Consider the random variable:$$X = F^{-1}(U),$$where $F^{-1}$ is the inverse of $F$.Show that the CDF of $X$ is $F(x)$, i.e.,$$F(x) = P(X \le x).$$This is called *inverse transform sampling*.In this part (Part A), you have to prove it.In a while, you will use it to learn how to sample from a continuous univariate random variable when all you have is samples from it.**Proof:** *Type your proof here. Delete that ```` line (it just makes some white space).*  Part BConsider the following data set $x_1,\dots,x_N$:<jupyter_code>data = np.loadtxt('hw2_p2_data.txt')
data = np.array(data)
fig, ax = plt.subplots()
ax.hist(data, alpha=0.5)
ax.set_xlabel('$x$')
ax.set_ylabel('Empirical PDF')<jupyter_output><empty_output><jupyter_text>Your goal is to generate a procedure that samples from the same distribution as the observed data. This is a variation of the standard *density estimation* problem.In general, this is a very difficult problem and we will see various ways to solve it later on.In this problem, you will develop a simple method that relies on the empirical CDF of the observed data.Needless to say, this method works only for one dimensional cases in which you have a lot of observations.The [empirical CDF](https://en.wikipedia.org/wiki/Empirical_distribution_function) of our data set, $x_1,\dots,x_N$ is defined to be:$$\hat{F}_N(x) = \frac{\text{Number of observations}\;\le x}{N} = \frac{1}{N}\sum_{i=1}^N 1_{[x_i,+\infty]}(x),$$where $1_A(x)$ is the [indicator function](https://en.wikipedia.org/wiki/Indicator_function) of the set $A$.Using the, so called, [strong law of large numbers](https://en.wikipedia.org/wiki/Law_of_large_numbersStrong_law), we can show that $\hat{F}_N(x)$ converges to the true CDF of the data.Complete the code that calculates the empirical CDF:<jupyter_code>def myECDF_base(x):
    """
    Make this code work if ``x`` is a simple scalar.
    
    :param x:    The point at which you want to observe the PDF.
    :param data: All the data you have observed (1D numpy array).
    :returns:    The value of the empirical CDF at ``x``.
    """
    N = data.shape[0]
    # Write your code here (delete the next line and return the right value)
    raise NotImplementedError('Implement me!')

# Vectorize your function (i.e., make it work with 1D numpy arrays).
# See this: https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.vectorize.html
myECDF = np.vectorize(myECDF_base)
# You can test your results by comparing the empirical CDF you can get with of scipy.stats
# The two should match almost exactly
hist_rv = st.rv_histogram(np.histogram(data, bins=1000))
fig, ax = plt.subplots()
# The range in which the x's takes values:
x_min = data.min()
x_max = data.max()
xx = np.linspace(x_min, x_max, 100)
ax.plot(xx, myECDF(xx), label='My empirical CDF')
ax.plot(xx, hist_rv.cdf(xx), '--', label='Empirical CDF from scipy.stats')
ax.set_xlabel('$x$')
ax.set_ylabel('Empirical CDF')
plt.legend(loc='best')<jupyter_output><empty_output><jupyter_text>Now complete the code that computes the inverse of the empirical CDF $\hat{F}^{-1}$.There are may ways of doing this.Let's do it in a way that will teach us something about the root finding toolbox of numpy (see [this]()).Mathematically, we wish to find a function $F^{-1}$ such that$$F(F^{-1}(u))) = u,$$for any $u\in[0,1]$ (the domain in which $F(x)$ takes values).It is obvious that $F^{-1}(u)$ is the solution to the *root finding* problem:$$F(x^*) = u.$$Since we know that $F$ is increasing, this problem must have a unique solution for any $u\in[0,1]$.To find this solution, we can use [Brent's method](https://en.wikipedia.org/wiki/Brent%27s_method).Please note, that the problem that this code solves is of the form:$$g(x^*) = 0.$$So, you will have to reformulate the original problem as:$$F(x^*) - u = 0.$$Study the [numpy implementation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.htmlscipy.optimize.brentq) of Brent's method and complete the following code:<jupyter_code>from scipy import optimize    # Gives you access to optimize.brentq



def myiECDF_base(u):
    """
    Evaluates the inverse of the empirical CDF.
    
    :param u:   A scalar at which to evaluate the function.
    :param data: All the data you have observed (1D numpy array).
    :returns:    The value of the inverse of the empirical CDF at ``u``.
    """
    # Write your code here
    # You will have to define a functioin to pass to optimize.brentq
    # You can define this function in here
    # (delete the next line and return the right value)
    raise NotImplementedError('Implement me!')

# Vectorize your function (i.e., make it work with 1D numpy arrays).
# See this: https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.vectorize.html
myiECDF = np.vectorize(myiECDF_base)
# You can test your results by comparing the inverse of the empirical CDF you can get with scipy.stats
# The two should match almost exactly
hist_rv = st.rv_histogram(np.histogram(data, bins=1000))
fig, ax = plt.subplots()
uu = np.linspace(1e-3, 1, 100)    # For convergence issues we cannot start at 0
ax.plot(uu, myiECDF(uu), label='Inverse of my empirical CDF')
ax.plot(uu, hist_rv.ppf(uu), '--', label='Inverse of empirical CDF from scipy.stats')
ax.set_xlabel('$x$')
ax.set_ylabel('Inverse of empirical CDF')
plt.legend(loc='best')<jupyter_output><empty_output><jupyter_text>Part CNow use the *inverse transform sampling* method to generate samples from same distribution as the original data.That is, you can now generate uniform samples:$$u_i \sim U([0,1]),$$and transform them as:$$\hat{x}_i = {\hat{F}}^{-1}(u_i).$$The $\hat{x}_i$'s generated in this way should have the same distribution of the data you started with.Verify this by comparing the histrogram of $1,000$ $\hat{x}_i$ samples with the original data of this problem.<jupyter_code># Write your code here.
# Feel free to copy paste code from above.<jupyter_output><empty_output><jupyter_text>Problem 3Demonstate that the central limit theorem holds for the data-generated distribution of Problem 2.Feel free to copy paste code from [handout 5](http://localhost:8888/notebooks/uq-course/handouts/handout_05.ipynb).At the very least you should1. Consider the sum:$$S_M = \frac{X_1+\dots+X_M}{M} = \frac{\hat{F}^{-1}(U_1)+\dots+\hat{F}^{-1}{U_M}}{M},$$and plot its histogram for various M;2. You should calculate $\mu$ and $\frac{\sigma^2}{M}$ from the right-hand-side of the central limit theorem and draw the corresponding Gaussian along with the histogram for 1;3. You should repeat 1 and 2 for various values of $M$ until you start getting convergent results.<jupyter_code># Put your code here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12333965844402277, 'fraction_numerical': 0.0150853889943074, 'mean_word_length': 4.865887590428492, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 9, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '431874', 'license': 'MIT', 'n_tokens_mistral': 3480, 'n_tokens_neox': 3328, 'n_words': 1280, 'path': 'homeworks/hw_02.ipynb', 'repo_name': 'bdlafleur/uq-course', 'size': '10540'}