-Since \( \widehat{\beta} = \widehat{\beta}(\boldsymbol{X}) \) is a function of random variables, -\( \widehat{\beta} \) itself must be a random variable. Thus it has -a pdf, call this function \( p(\boldsymbol{t}) \). The aim of the bootstrap is to -estimate \( p(\boldsymbol{t}) \) by the relative frequency of -\( \widehat{\beta} \). You can think of this as using a histogram -in the place of \( p(\boldsymbol{t}) \). If the relative frequency closely -resembles \( p(\vec{t}) \), then using numerics, it is straight forward to -estimate all the interesting parameters of \( p(\boldsymbol{t}) \) using point -estimators. +Suppose we have a PDF \( p(x) \) from which we generate a series \( N \) +of averages \( \langle x_i \rangle \). Each mean value \( \langle x_i \rangle \) +is viewed as the average of a specific measurement, e.g., throwing +dice 100 times and then taking the average value, or producing a certain +amount of random numbers. +For notational ease, we set \( \langle x_i \rangle=x_i \) in the discussion +which follows. + +
+If we compute the mean \( z \) of \( m \) such mean values \( x_i \) +$$ + z=\frac{x_1+x_2+\dots+x_m}{m}, +$$ + +the question we pose is which is the PDF of the new variable \( z \).
@@ -292,7 +318,7 @@ estimators.
-In the case that \( \widehat{\beta} \) has -more than one component, and the components are independent, we use the -same estimator on each component separately. If the probability -density function of \( X_i \), \( p(x) \), had been known, then it would have -been straightforward to do this by: +The probability of obtaining an average value \( z \) is the product of the +probabilities of obtaining arbitrary individual mean values \( x_i \), +but with the constraint that the average is \( z \). We can express this through +the following expression +$$ + \tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m) + \delta(z-\frac{x_1+x_2+\dots+x_m}{m}), +$$ -
@@ -298,7 +316,7 @@ idea is to use the relative frequency of \( \widehat{\beta}^* \)
-But -unless there is enough information available about the process that -generated \( X_1,X_2,\cdots,X_n \), \( p(x) \) is in general -unknown. Therefore, Efron in 1979 asked the -question: What if we replace \( p(x) \) by the relative frequency -of the observation \( X_i \)? +If we use the integral expression for the \( \delta \)-function -
-If we draw observations in accordance with -the relative frequency of the observations, will we obtain the same -result in some asymptotic sense? The answer is yes. +$$ + \delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty} + dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)}, +$$ + +and inserting \( e^{i\mu q-i\mu q} \) where \( \mu \) is the mean value +we arrive at +$$ + \tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty} + dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty} + dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m, +$$ + +with the integral over \( x \) resulting in + +$$ + \int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}= + \int_{-\infty}^{\infty}dxp(x) + \left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right]. +$$
@@ -294,7 +325,7 @@ result in some asymptotic sense? The answer is yes.
-The independent bootstrap works like this: +The second term on the rhs disappears since this is just the mean and +employing the definition of \( \sigma^2 \) we have +$$ + \int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}= + 1-\frac{q^2\sigma^2}{2m^2}+\dots, +$$ -
@@ -301,7 +327,7 @@ example, if you are interested in estimating the variance of \( \widehat
-The following code starts with a Gaussian distribution with mean value -\( \mu =100 \) and variance \( \sigma=15 \). We use this to generate the data -used in the bootstrap analysis. The bootstrap analysis returns a data -set after a given number of bootstrap operations (as many as we have -data points). This data set consists of estimated mean values for each -bootstrap operation. The histogram generated by the bootstrap method -shows that the distribution for these mean values is also a Gaussian, -centered around the mean value \( \mu=100 \) but with standard deviation -\( \sigma/\sqrt{n} \), where \( n \) is the number of bootstrap samples (in -this case the same as the number of original data points). The value -of the standard deviation is what we expect from the central limit -theorem. +Thus, the central limit theorem states that the PDF \( \tilde{p}(z) \) of +the average of \( m \) random values corresponding to a PDF \( p(x) \) +is a normal distribution whose mean is the +mean value of the PDF \( p(x) \) and whose variance is the variance +of the PDF \( p(x) \) divided by \( m \), the number of values used to compute \( z \).
+The theorem is satisfied by a large class of PDFs. Note however that for a +finite \( m \), it is not always possible to find a closed expression for +\( \tilde{p}(x) \). +The central limit theorem leads then to the well-known expression for the +standard deviation, given by - -
from numpy import *
-from numpy.random import randint, randn
-from time import time
-import matplotlib.mlab as mlab
-import matplotlib.pyplot as plt
+$$
+ \sigma_m=
+\frac{\sigma}{\sqrt{m}}.
+$$
-# Returns mean of bootstrap samples
-def stat(data):
- return mean(data)
+
+The latter is true only if the average value is known exactly. This is obtained in the limit
+\( m\rightarrow \infty \) only. Because the mean and the variance are measured quantities we obtain
+the familiar expression in statistics
+$$
+ \sigma_m\approx
+\frac{\sigma}{\sqrt{m-1}}.
+$$
-# Bootstrap algorithm
-def bootstrap(data, statistic, R):
- t = zeros(R); n = len(data); inds = arange(n); t0 = time()
- # non-parametric bootstrap
- for i in range(R):
- t[i] = statistic(data[randint(0,n,n)])
+
+In many cases however the above estimate for the standard deviation,
+in particular if correlations are strong, may be too simplistic.
- # analysis
- print("Runtime: %g sec" % (time()-t0)); print("Bootstrap Statistics :")
- print("original bias std. error")
- print("%8g %8g %14g %15g" % (statistic(data), std(data),mean(t),std(t)))
- return t
-
-
-mu, sigma = 100, 15
-datapoints = 10000
-x = mu + sigma*random.randn(datapoints)
-# bootstrap returns the data sample
-t = bootstrap(x, stat, datapoints)
-# the histogram of the bootstrapped data
-n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
-
-# add a 'best fit' line
-y = mlab.normpdf( binsboot, mean(t), std(t))
-lt = plt.plot(binsboot, y, 'r--', linewidth=1)
-plt.xlabel('Smarts')
-plt.ylabel('Probability')
-plt.axis([99.5, 100.6, 0, 3.0])
-plt.grid(True)
-
-plt.show()
-
@@ -340,7 +333,7 @@ plt.show()
-We will discuss the bias-variance tradeoff in the context of -continuous predictions such as regression. However, many of the -intuitions and ideas discussed here also carry over to classification -tasks. Consider a dataset \( \mathcal{L} \) consisting of the data -\( \mathbf{X}_\mathcal{L}=\{(y_j, \boldsymbol{x}_j), j=0\ldots n-1\} \). - -
-Let us assume that the true data is generated from a noisy model - -$$ -\boldsymbol{y}=f(\boldsymbol{x}) + \boldsymbol{\epsilon} -$$ - -
-where \( \epsilon \) is normally distributed with mean zero and standard deviation \( \sigma^2 \). - -
-In our derivation of the ordinary least squares method we defined then -an approximation to the function \( f \) in terms of the parameters -\( \boldsymbol{\beta} \) and the design matrix \( \boldsymbol{X} \) which embody our model, -that is \( \boldsymbol{\tilde{y}}=\boldsymbol{X}\boldsymbol{\beta} \). - -
-Thereafter we found the parameters \( \boldsymbol{\beta} \) by optimizing the means squared error via the so-called cost function -$$ -C(\boldsymbol{X},\boldsymbol{\beta}) =\frac{1}{n}\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2=\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]. -$$ - -
-We can rewrite this as -$$ -\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\frac{1}{n}\sum_i(f_i-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2+\frac{1}{n}\sum_i(\tilde{y}_i-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2+\sigma^2. -$$ - -
-The three terms represent the square of the bias of the learning -method, which can be thought of as the error caused by the simplifying -assumptions built into the method. The second term represents the -variance of the chosen model and finally the last terms is variance of -the error \( \boldsymbol{\epsilon} \). - -
-To derive this equation, we need to recall that the variance of \( \boldsymbol{y} \) and \( \boldsymbol{\epsilon} \) are both equal to \( \sigma^2 \). The mean value of \( \boldsymbol{\epsilon} \) is by definition equal to zero. Furthermore, the function \( f \) is not a stochastics variable, idem for \( \boldsymbol{\tilde{y}} \). -We use a more compact notation in terms of the expectation value -$$ -\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathbb{E}\left[(\boldsymbol{f}+\boldsymbol{\epsilon}-\boldsymbol{\tilde{y}})^2\right], -$$ - -and adding and subtracting \( \mathbb{E}\left[\boldsymbol{\tilde{y}}\right] \) we get -$$ -\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathbb{E}\left[(\boldsymbol{f}+\boldsymbol{\epsilon}-\boldsymbol{\tilde{y}}+\mathbb{E}\left[\boldsymbol{\tilde{y}}\right]-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2\right], -$$ - -which, using the abovementioned expectation values can be rewritten as -$$ -\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathbb{E}\left[(\boldsymbol{y}-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2\right]+\mathrm{Var}\left[\boldsymbol{\tilde{y}}\right]+\sigma^2, -$$ - -that is the rewriting in terms of the so-called bias, the variance of the model \( \boldsymbol{\tilde{y}} \) and the variance of \( \boldsymbol{\epsilon} \). +
@@ -342,7 +301,7 @@ that is the rewriting in terms of the so-called bias, the variance of the model
+Since \( \widehat{\beta} = \widehat{\beta}(\boldsymbol{X}) \) is a function of random variables, +\( \widehat{\beta} \) itself must be a random variable. Thus it has +a pdf, call this function \( p(\boldsymbol{t}) \). The aim of the bootstrap is to +estimate \( p(\boldsymbol{t}) \) by the relative frequency of +\( \widehat{\beta} \). You can think of this as using a histogram +in the place of \( p(\boldsymbol{t}) \). If the relative frequency closely +resembles \( p(\vec{t}) \), then using numerics, it is straight forward to +estimate all the interesting parameters of \( p(\boldsymbol{t}) \) using point +estimators. - -
import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.linear_model import LinearRegression, Ridge, Lasso
-from sklearn.preprocessing import PolynomialFeatures
-from sklearn.model_selection import train_test_split
-from sklearn.pipeline import make_pipeline
-from sklearn.utils import resample
-
-np.random.seed(2018)
-
-n = 500
-n_boostraps = 100
-degree = 18 # A quite high value, just to show.
-noise = 0.1
-
-# Make data set.
-x = np.linspace(-1, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + np.random.normal(0, 0.1, x.shape)
-
-# Hold out some test data that is never used in training.
-x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
-
-# Combine x transformation and model into one operation.
-# Not neccesary, but convenient.
-model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))
-
-# The following (m x n_bootstraps) matrix holds the column vectors y_pred
-# for each bootstrap iteration.
-y_pred = np.empty((y_test.shape[0], n_boostraps))
-for i in range(n_boostraps):
- x_, y_ = resample(x_train, y_train)
-
- # Evaluate the new model on the same test data each time.
- y_pred[:, i] = model.fit(x_, y_).predict(x_test).ravel()
-
-# Note: Expectations and variances taken w.r.t. different training
-# data sets, hence the axis=1. Subsequent means are taken across the test data
-# set in order to obtain a total value, but before this we have error/bias/variance
-# calculated per data point in the test set.
-# Note 2: The use of keepdims=True is important in the calculation of bias as this
-# maintains the column vector form. Dropping this yields very unexpected results.
-error = np.mean( np.mean((y_test - y_pred)**2, axis=1, keepdims=True) )
-bias = np.mean( (y_test - np.mean(y_pred, axis=1, keepdims=True))**2 )
-variance = np.mean( np.var(y_pred, axis=1, keepdims=True) )
-print('Error:', error)
-print('Bias^2:', bias)
-print('Var:', variance)
-print('{} >= {} + {} = {}'.format(error, bias, variance, bias+variance))
-
-plt.plot(x[::5, :], y[::5, :], label='f(x)')
-plt.scatter(x_test, y_test, label='Data points')
-plt.scatter(x_test, np.mean(y_pred, axis=1), label='Pred')
-plt.legend()
-plt.show()
-
@@ -338,7 +312,7 @@ plt.show()
+In the case that \( \widehat{\beta} \) has +more than one component, and the components are independent, we use the +same estimator on each component separately. If the probability +density function of \( X_i \), \( p(x) \), had been known, then it would have +been straightforward to do this by: - -
import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.linear_model import LinearRegression, Ridge, Lasso
-from sklearn.preprocessing import PolynomialFeatures
-from sklearn.model_selection import train_test_split
-from sklearn.pipeline import make_pipeline
-from sklearn.utils import resample
+
+- Drawing lots of numbers from \( p(x) \), suppose we call one such set of numbers \( (X_1^*, X_2^*, \cdots, X_n^*) \).
+- Then using these numbers, we could compute a replica of \( \widehat{\beta} \) called \( \widehat{\beta}^* \).
+
-np.random.seed(2018)
+By repeated use of the above two points, many
+estimates of \( \widehat{\beta} \) can be obtained. The
+idea is to use the relative frequency of \( \widehat{\beta}^* \)
+(think of a histogram) as an estimate of \( p(\boldsymbol{t}) \).
-n = 40
-n_boostraps = 100
-maxdegree = 14
-
-
-# Make data set.
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-error = np.zeros(maxdegree)
-bias = np.zeros(maxdegree)
-variance = np.zeros(maxdegree)
-polydegree = np.zeros(maxdegree)
-x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
-
-for degree in range(maxdegree):
- model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))
- y_pred = np.empty((y_test.shape[0], n_boostraps))
- for i in range(n_boostraps):
- x_, y_ = resample(x_train, y_train)
- y_pred[:, i] = model.fit(x_, y_).predict(x_test).ravel()
-
- polydegree[degree] = degree
- error[degree] = np.mean( np.mean((y_test - y_pred)**2, axis=1, keepdims=True) )
- bias[degree] = np.mean( (y_test - np.mean(y_pred, axis=1, keepdims=True))**2 )
- variance[degree] = np.mean( np.var(y_pred, axis=1, keepdims=True) )
- print('Polynomial degree:', degree)
- print('Error:', error[degree])
- print('Bias^2:', bias[degree])
- print('Var:', variance[degree])
- print('{} >= {} + {} = {}'.format(error[degree], bias[degree], variance[degree], bias[degree]+variance[degree]))
-
-plt.plot(polydegree, error, label='Error')
-plt.plot(polydegree, bias, label='bias')
-plt.plot(polydegree, variance, label='Variance')
-plt.legend()
-plt.show()
-
@@ -329,6 +317,8 @@ plt.show()
- + -
-The bias-variance tradeoff summarizes the fundamental tension in -machine learning, particularly supervised learning, between the -complexity of a model and the amount of training data needed to train -it. Since data is often limited, in practice it is often useful to -use a less-complex model with higher bias, that is a model whose asymptotic -performance is worse than another model because it is easier to -train and less sensitive to sampling noise arising from having a -finite-sized training dataset (smaller variance). +But +unless there is enough information available about the process that +generated \( X_1,X_2,\cdots,X_n \), \( p(x) \) is in general +unknown. Therefore, Efron in 1979 asked the +question: What if we replace \( p(x) \) by the relative frequency +of the observation \( X_i \)?
-The above equations tell us that in -order to minimize the expected test error, we need to select a -statistical learning method that simultaneously achieves low variance -and low bias. Note that variance is inherently a nonnegative quantity, -and squared bias is also nonnegative. Hence, we see that the expected -test MSE can never lie below \( Var(\epsilon) \), the irreducible error. - -
-What do we mean by the variance and bias of a statistical learning -method? The variance refers to the amount by which our model would change if we -estimated it using a different training data set. Since the training -data are used to fit the statistical learning method, different -training data sets will result in a different estimate. But ideally the -estimate for our model should not vary too much between training -sets. However, if a method has high variance then small changes in -the training data can result in large changes in the model. In general, more -flexible statistical methods have higher variance. - -
-You may also find this recent article of interest. +If we draw observations in accordance with +the relative frequency of the observations, will we obtain the same +result in some asymptotic sense? The answer is yes.
@@ -311,6 +312,9 @@ You may also find this recent 39
+The independent bootstrap works like this: - -
"""
-============================
-Underfitting vs. Overfitting
-============================
+
+- Draw with replacement \( n \) numbers for the observed variables \( \boldsymbol{x} = (x_1,x_2,\cdots,x_n) \).
+- Define a vector \( \boldsymbol{x}^* \) containing the values which were drawn from \( \boldsymbol{x} \).
+- Using the vector \( \boldsymbol{x}^* \) compute \( \widehat{\beta}^* \) by evaluating \( \widehat \beta \) under the observations \( \boldsymbol{x}^* \).
+- Repeat this process \( k \) times.
+
-This example demonstrates the problems of underfitting and overfitting and
-how we can use linear regression with polynomial features to approximate
-nonlinear functions. The plot shows the function that we want to approximate,
-which is a part of the cosine function. In addition, the samples from the
-real function and the approximations of different models are displayed. The
-models have polynomial features of different degrees. We can see that a
-linear function (polynomial with degree 1) is not sufficient to fit the
-training samples. This is called **underfitting**. A polynomial of degree 4
-approximates the true function almost perfectly. However, for higher degrees
-the model will **overfit** the training data, i.e. it learns the noise of the
-training data.
-We evaluate quantitatively **overfitting** / **underfitting** by using
-cross-validation. We calculate the mean squared error (MSE) on the validation
-set, the higher, the less likely the model generalizes correctly from the
-training data.
-"""
+When you are done, you can draw a histogram of the relative frequency
+of \( \widehat \beta^* \). This is your estimate of the probability
+distribution \( p(t) \). Using this probability distribution you can
+estimate any statistics thereof. In principle you never draw the
+histogram of the relative frequency of \( \widehat{\beta}^* \). Instead
+you use the estimators corresponding to the statistic of interest. For
+example, if you are interested in estimating the variance of \( \widehat
+\beta \), apply the etsimator \( \widehat \sigma^2 \) to the values
+\( \widehat \beta^* \).
-print(__doc__)
-
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.pipeline import Pipeline
-from sklearn.preprocessing import PolynomialFeatures
-from sklearn.linear_model import LinearRegression
-from sklearn.model_selection import cross_val_score
-
-
-def true_fun(X):
- return np.cos(1.5 * np.pi * X)
-
-np.random.seed(0)
-
-n_samples = 30
-degrees = [1, 4, 15]
-
-X = np.sort(np.random.rand(n_samples))
-y = true_fun(X) + np.random.randn(n_samples) * 0.1
-
-plt.figure(figsize=(14, 5))
-for i in range(len(degrees)):
- ax = plt.subplot(1, len(degrees), i + 1)
- plt.setp(ax, xticks=(), yticks=())
-
- polynomial_features = PolynomialFeatures(degree=degrees[i],
- include_bias=False)
- linear_regression = LinearRegression()
- pipeline = Pipeline([("polynomial_features", polynomial_features),
- ("linear_regression", linear_regression)])
- pipeline.fit(X[:, np.newaxis], y)
-
- # Evaluate the models using crossvalidation
- scores = cross_val_score(pipeline, X[:, np.newaxis], y,
- scoring="neg_mean_squared_error", cv=10)
-
- X_test = np.linspace(0, 1, 100)
- plt.plot(X_test, pipeline.predict(X_test[:, np.newaxis]), label="Model")
- plt.plot(X_test, true_fun(X_test), label="True function")
- plt.scatter(X, y, edgecolor='b', s=20, label="Samples")
- plt.xlabel("x")
- plt.ylabel("y")
- plt.xlim((0, 1))
- plt.ylim((-2, 2))
- plt.legend(loc="best")
- plt.title("Degree {}\nMSE = {:.2e}(+/- {:.2e})".format(
- degrees[i], -scores.mean(), scores.std()))
-plt.show()
-
@@ -352,6 +318,10 @@ plt.show()
- + -
-When the repetitive splitting of the data set is done randomly, -samples may accidently end up in a fast majority of the splits in -either training or test set. Such samples may have an unbalanced -influence on either model building or prediction evaluation. To avoid -this \( k \)-fold cross-validation structures the data splitting. The -samples are divided into \( k \) more or less equally sized exhaustive and -mutually exclusive subsets. In turn (at each split) one of these -subsets plays the role of the test set while the union of the -remaining subsets constitutes the training set. Such a splitting -warrants a balanced representation of each sample in both training and -test set over the splits. Still the division into the \( k \) subsets -involves a degree of randomness. This may be fully excluded when -choosing \( k=n \). This particular case is referred to as leave-one-out -cross-validation (LOOCV). +The following code starts with a Gaussian distribution with mean value +\( \mu =100 \) and variance \( \sigma=15 \). We use this to generate the data +used in the bootstrap analysis. The bootstrap analysis returns a data +set after a given number of bootstrap operations (as many as we have +data points). This data set consists of estimated mean values for each +bootstrap operation. The histogram generated by the bootstrap method +shows that the distribution for these mean values is also a Gaussian, +centered around the mean value \( \mu=100 \) but with standard deviation +\( \sigma/\sqrt{n} \), where \( n \) is the number of bootstrap samples (in +this case the same as the number of original data points). The value +of the standard deviation is what we expect from the central limit +theorem. + +
+ + +
from numpy import *
+from numpy.random import randint, randn
+from time import time
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+
+# Returns mean of bootstrap samples # Alternatively, we can run it using Scikit-Learn's function resample # See the examples below
+
+def statistics(data):
+ return mean(data)
+
+
+# Bootstrap algorithm
+def bootstrap(data, statistic, R):
+ t = zeros(R); n = len(data); inds = arange(n); t0 = time()
+ # non-parametric bootstrap
+ for i in range(R):
+ t[i] = statistic(data[randint(0,n,n)])
+
+ # analysis
+ print("Runtime: %g sec" % (time()-t0)); print("Bootstrap Statistics :")
+ print("original bias std. error")
+ print("%8g %8g %14g %15g" % (statistic(data), std(data),mean(t),std(t)))
+ return t
+
+
+mu, sigma = 100, 15
+datapoints = 10000
+x = mu + sigma*random.randn(datapoints)
+# bootstrap returns the data sample
+t = bootstrap(x, statistics, datapoints)
++We see that our new variance and from that the standard deviation, agrees with the central limit theorem.
@@ -293,6 +349,11 @@ cross-validation (LOOCV).
- + -
-
# the histogram of the bootstrapped data
+n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
-$$
-\begin{align*}
-\boldsymbol{\beta}_{-i}(\lambda) & = ( \boldsymbol{X}_{-i, \ast}^{T}
-\boldsymbol{X}_{-i, \ast} + \lambda \boldsymbol{I}_{pp})^{-1}
-\boldsymbol{X}_{-i, \ast}^{T} \boldsymbol{y}_{-i}
-\end{align*}
-$$
-
-
-
@@ -303,6 +311,12 @@ $$
-For the various values of \( k \) +We will discuss the bias-variance tradeoff in the context of +continuous predictions such as regression. However, many of the +intuitions and ideas discussed here also carry over to classification +tasks. Consider a dataset \( \mathcal{L} \) consisting of the data +\( \mathbf{X}_\mathcal{L}=\{(y_j, \boldsymbol{x}_j), j=0\ldots n-1\} \). -
+Let us assume that the true data is generated from a noisy model -
+where \( \epsilon \) is normally distributed with mean zero and standard deviation \( \sigma^2 \). +
+In our derivation of the ordinary least squares method we defined then +an approximation to the function \( f \) in terms of the parameters +\( \boldsymbol{\beta} \) and the design matrix \( \boldsymbol{X} \) which embody our model, +that is \( \boldsymbol{\tilde{y}}=\boldsymbol{X}\boldsymbol{\beta} \). + +
+Thereafter we found the parameters \( \boldsymbol{\beta} \) by optimizing the means squared error via the so-called cost function +$$ +C(\boldsymbol{X},\boldsymbol{\beta}) =\frac{1}{n}\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2=\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]. +$$ + +
+We can rewrite this as +$$ +\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\frac{1}{n}\sum_i(f_i-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2+\frac{1}{n}\sum_i(\tilde{y}_i-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2+\sigma^2. +$$ + +
+The three terms represent the square of the bias of the learning +method, which can be thought of as the error caused by the simplifying +assumptions built into the method. The second term represents the +variance of the chosen model and finally the last terms is variance of +the error \( \boldsymbol{\epsilon} \). + +
+To derive this equation, we need to recall that the variance of \( \boldsymbol{y} \) and \( \boldsymbol{\epsilon} \) are both equal to \( \sigma^2 \). The mean value of \( \boldsymbol{\epsilon} \) is by definition equal to zero. Furthermore, the function \( f \) is not a stochastics variable, idem for \( \boldsymbol{\tilde{y}} \). +We use a more compact notation in terms of the expectation value +$$ +\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathbb{E}\left[(\boldsymbol{f}+\boldsymbol{\epsilon}-\boldsymbol{\tilde{y}})^2\right], +$$ + +and adding and subtracting \( \mathbb{E}\left[\boldsymbol{\tilde{y}}\right] \) we get +$$ +\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathbb{E}\left[(\boldsymbol{f}+\boldsymbol{\epsilon}-\boldsymbol{\tilde{y}}+\mathbb{E}\left[\boldsymbol{\tilde{y}}\right]-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2\right], +$$ + +which, using the abovementioned expectation values can be rewritten as +$$ +\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathbb{E}\left[(\boldsymbol{y}-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2\right]+\mathrm{Var}\left[\boldsymbol{\tilde{y}}\right]+\sigma^2, +$$ + +that is the rewriting in terms of the so-called bias, the variance of the model \( \boldsymbol{\tilde{y}} \) and the variance of \( \boldsymbol{\epsilon} \). + +
-The code here uses Ridge regression with cross-validation (CV) resampling and \( k \)-fold CV in order to fit a specific polynomial. +
-
import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
+import matplotlib.pyplot as plt
+import numpy as np
+from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
+from sklearn.utils import resample
-# A seed just to ensure that the random numbers are the same for every run.
-# Useful for eventual debugging.
-np.random.seed(3155)
+np.random.seed(2018)
-# Generate the data.
-nsamples = 100
-x = np.random.randn(nsamples)
-y = 3*x**2 + np.random.randn(nsamples)
+n = 500
+n_boostraps = 100
+degree = 18 # A quite high value, just to show.
+noise = 0.1
-## Cross-validation on Ridge regression using KFold only
+# Make data set.
+x = np.linspace(-1, 3, n).reshape(-1, 1)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 6)
+# Hold out some test data that is never used in training.
+x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
+# Combine x transformation and model into one operation.
+# Not neccesary, but convenient.
+model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
+# The following (m x n_bootstraps) matrix holds the column vectors y_pred
+# for each bootstrap iteration.
+y_pred = np.empty((y_test.shape[0], n_boostraps))
+for i in range(n_boostraps):
+ x_, y_ = resample(x_train, y_train)
-# Perform the cross-validation to estimate MSE
-scores_KFold = np.zeros((nlambdas, k))
+ # Evaluate the new model on the same test data each time.
+ y_pred[:, i] = model.fit(x_, y_).predict(x_test).ravel()
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- j = 0
- for train_inds, test_inds in kfold.split(x):
- xtrain = x[train_inds]
- ytrain = y[train_inds]
-
- xtest = x[test_inds]
- ytest = y[test_inds]
-
- Xtrain = poly.fit_transform(xtrain[:, np.newaxis])
- ridge.fit(Xtrain, ytrain[:, np.newaxis])
-
- Xtest = poly.fit_transform(xtest[:, np.newaxis])
- ypred = ridge.predict(Xtest)
-
- scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)
-
- j += 1
- i += 1
-
-
-estimated_mse_KFold = np.mean(scores_KFold, axis = 1)
-
-## Cross-validation using cross_val_score from sklearn along with KFold
-
-# kfold is an instance initialized above as:
-# kfold = KFold(n_splits = k)
-
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
-
- X = poly.fit_transform(x[:, np.newaxis])
- estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)
-
- # cross_val_score return an array containing the estimated negative mse for every fold.
- # we have to the the mean of every array in order to get an estimate of the mse of the model
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-
- i += 1
-
-## Plot and compare the slightly different ways to perform cross-validation
-
-plt.figure()
-
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')
-
-plt.xlabel('log10(lambda)')
-plt.ylabel('mse')
+# Note: Expectations and variances taken w.r.t. different training
+# data sets, hence the axis=1. Subsequent means are taken across the test data
+# set in order to obtain a total value, but before this we have error/bias/variance
+# calculated per data point in the test set.
+# Note 2: The use of keepdims=True is important in the calculation of bias as this
+# maintains the column vector form. Dropping this yields very unexpected results.
+error = np.mean( np.mean((y_test - y_pred)**2, axis=1, keepdims=True) )
+bias = np.mean( (y_test - np.mean(y_pred, axis=1, keepdims=True))**2 )
+variance = np.mean( np.var(y_pred, axis=1, keepdims=True) )
+print('Error:', error)
+print('Bias^2:', bias)
+print('Var:', variance)
+print('{} >= {} + {} = {}'.format(error, bias, variance, bias+variance))
+plt.plot(x[::5, :], y[::5, :], label='f(x)')
+plt.scatter(x_test, y_test, label='Data points')
+plt.scatter(x_test, np.mean(y_pred, axis=1), label='Pred')
plt.legend()
-
plt.show()
@@ -369,6 +351,14 @@ plt.show()
-
# Common imports
-import os
+import matplotlib.pyplot as plt
import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge, Lasso
+from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
+from sklearn.pipeline import make_pipeline
from sklearn.utils import resample
-from sklearn.metrics import mean_squared_error
-# Where to save the figures and data files
-PROJECT_ROOT_DIR = "Results"
-FIGURE_ID = "Results/FigureFiles"
-DATA_ID = "DataFiles/"
-if not os.path.exists(PROJECT_ROOT_DIR):
- os.mkdir(PROJECT_ROOT_DIR)
+np.random.seed(2018)
-if not os.path.exists(FIGURE_ID):
- os.makedirs(FIGURE_ID)
+n = 40
+n_boostraps = 100
+maxdegree = 14
-if not os.path.exists(DATA_ID):
- os.makedirs(DATA_ID)
-def image_path(fig_id):
- return os.path.join(FIGURE_ID, fig_id)
+# Make data set.
+x = np.linspace(-3, 3, n).reshape(-1, 1)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
+error = np.zeros(maxdegree)
+bias = np.zeros(maxdegree)
+variance = np.zeros(maxdegree)
+polydegree = np.zeros(maxdegree)
+x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
-def data_path(dat_id):
- return os.path.join(DATA_ID, dat_id)
+for degree in range(maxdegree):
+ model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))
+ y_pred = np.empty((y_test.shape[0], n_boostraps))
+ for i in range(n_boostraps):
+ x_, y_ = resample(x_train, y_train)
+ y_pred[:, i] = model.fit(x_, y_).predict(x_test).ravel()
-def save_fig(fig_id):
- plt.savefig(image_path(fig_id) + ".png", format='png')
+ polydegree[degree] = degree
+ error[degree] = np.mean( np.mean((y_test - y_pred)**2, axis=1, keepdims=True) )
+ bias[degree] = np.mean( (y_test - np.mean(y_pred, axis=1, keepdims=True))**2 )
+ variance[degree] = np.mean( np.var(y_pred, axis=1, keepdims=True) )
+ print('Polynomial degree:', degree)
+ print('Error:', error[degree])
+ print('Bias^2:', bias[degree])
+ print('Var:', variance[degree])
+ print('{} >= {} + {} = {}'.format(error[degree], bias[degree], variance[degree], bias[degree]+variance[degree]))
-infile = open(data_path("EoS.csv"),'r')
-
-# Read the EoS data as csv file and organize the data into two arrays with density and energies
-EoS = pd.read_csv(infile, names=('Density', 'Energy'))
-EoS['Energy'] = pd.to_numeric(EoS['Energy'], errors='coerce')
-EoS = EoS.dropna()
-Energies = EoS['Energy']
-Density = EoS['Density']
-# The design matrix now as function of various polytrops
-
-Maxpolydegree = 30
-X = np.zeros((len(Density),Maxpolydegree))
-X[:,0] = 1.0
-testerror = np.zeros(Maxpolydegree)
-trainingerror = np.zeros(Maxpolydegree)
-polynomial = np.zeros(Maxpolydegree)
-
-trials = 100
-for polydegree in range(1, Maxpolydegree):
- polynomial[polydegree] = polydegree
- for degree in range(polydegree):
- X[:,degree] = Density**(degree/3.0)
-
-# loop over trials in order to estimate the expectation value of the MSE
- testerror[polydegree] = 0.0
- trainingerror[polydegree] = 0.0
- for samples in range(trials):
- x_train, x_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2)
- model = LinearRegression(fit_intercept=True).fit(x_train, y_train)
- ypred = model.predict(x_train)
- ytilde = model.predict(x_test)
- testerror[polydegree] += mean_squared_error(y_test, ytilde)
- trainingerror[polydegree] += mean_squared_error(y_train, ypred)
-
- testerror[polydegree] /= trials
- trainingerror[polydegree] /= trials
- print("Degree of polynomial: %3d"% polynomial[polydegree])
- print("Mean squared error on training data: %.8f" % trainingerror[polydegree])
- print("Mean squared error on test data: %.8f" % testerror[polydegree])
-
-plt.plot(polynomial, np.log10(trainingerror), label='Training Error')
-plt.plot(polynomial, np.log10(testerror), label='Test Error')
-plt.xlabel('Polynomial degree')
-plt.ylabel('log10[MSE]')
+plt.plot(polydegree, error, label='Error')
+plt.plot(polydegree, bias, label='bias')
+plt.plot(polydegree, variance, label='Variance')
plt.legend()
plt.show()
@@ -355,6 +342,13 @@ plt.show()
+The bias-variance tradeoff summarizes the fundamental tension in +machine learning, particularly supervised learning, between the +complexity of a model and the amount of training data needed to train +it. Since data is often limited, in practice it is often useful to +use a less-complex model with higher bias, that is a model whose asymptotic +performance is worse than another model because it is easier to +train and less sensitive to sampling noise arising from having a +finite-sized training dataset (smaller variance). - -
# Common imports
-import os
-import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-from sklearn.linear_model import LinearRegression, Ridge, Lasso
-from sklearn.metrics import mean_squared_error
-from sklearn.model_selection import KFold
-from sklearn.model_selection import cross_val_score
+
+The above equations tell us that in
+order to minimize the expected test error, we need to select a
+statistical learning method that simultaneously achieves low variance
+and low bias. Note that variance is inherently a nonnegative quantity,
+and squared bias is also nonnegative. Hence, we see that the expected
+test MSE can never lie below \( Var(\epsilon) \), the irreducible error.
+
+What do we mean by the variance and bias of a statistical learning
+method? The variance refers to the amount by which our model would change if we
+estimated it using a different training data set. Since the training
+data are used to fit the statistical learning method, different
+training data sets will result in a different estimate. But ideally the
+estimate for our model should not vary too much between training
+sets. However, if a method has high variance then small changes in
+the training data can result in large changes in the model. In general, more
+flexible statistical methods have higher variance.
-# Where to save the figures and data files
-PROJECT_ROOT_DIR = "Results"
-FIGURE_ID = "Results/FigureFiles"
-DATA_ID = "DataFiles/"
+
+You may also find this recent article of interest.
-if not os.path.exists(PROJECT_ROOT_DIR):
- os.mkdir(PROJECT_ROOT_DIR)
-
-if not os.path.exists(FIGURE_ID):
- os.makedirs(FIGURE_ID)
-
-if not os.path.exists(DATA_ID):
- os.makedirs(DATA_ID)
-
-def image_path(fig_id):
- return os.path.join(FIGURE_ID, fig_id)
-
-def data_path(dat_id):
- return os.path.join(DATA_ID, dat_id)
-
-def save_fig(fig_id):
- plt.savefig(image_path(fig_id) + ".png", format='png')
-
-infile = open(data_path("EoS.csv"),'r')
-
-# Read the EoS data as csv file and organize the data into two arrays with density and energies
-EoS = pd.read_csv(infile, names=('Density', 'Energy'))
-EoS['Energy'] = pd.to_numeric(EoS['Energy'], errors='coerce')
-EoS = EoS.dropna()
-Energies = EoS['Energy']
-Density = EoS['Density']
-# The design matrix now as function of various polytrops
-
-Maxpolydegree = 30
-X = np.zeros((len(Density),Maxpolydegree))
-X[:,0] = 1.0
-estimated_mse_sklearn = np.zeros(Maxpolydegree)
-polynomial = np.zeros(Maxpolydegree)
-k =5
-kfold = KFold(n_splits = k)
-
-for polydegree in range(1, Maxpolydegree):
- polynomial[polydegree] = polydegree
- for degree in range(polydegree):
- X[:,degree] = Density**(degree/3.0)
- OLS = LinearRegression()
-# loop over trials in order to estimate the expectation value of the MSE
- estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold)
-#[:, np.newaxis]
- estimated_mse_sklearn[polydegree] = np.mean(-estimated_mse_folds)
-
-plt.plot(polynomial, np.log10(estimated_mse_sklearn), label='Test Error')
-plt.xlabel('Polynomial degree')
-plt.ylabel('log10[MSE]')
-plt.legend()
-plt.show()
-
@@ -343,6 +324,13 @@ plt.show()
+Suppose we have a PDF \( p(x) \) from which we generate a series \( N \) +of averages \( \langle x_i \rangle \). Each mean value \( \langle x_i \rangle \) +is viewed as the average of a specific measurement, e.g., throwing +dice 100 times and then taking the average value, or producing a certain +amount of random numbers. +For notational ease, we set \( \langle x_i \rangle=x_i \) in the discussion +which follows. + +
+If we compute the mean \( z \) of \( m \) such mean values \( x_i \) + +
+$$
+ z=\frac{x_1+x_2+\dots+x_m}{m},
+$$
+
+
+the question we pose is which is the PDF of the new variable \( z \).
+
+The probability of obtaining an average value \( z \) is the product of the +probabilities of obtaining arbitrary individual mean values \( x_i \), +but with the constraint that the average is \( z \). We can express this through +the following expression +
+$$
+ \tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
+$$
+
+
+where the \( \delta \)-function enbodies the constraint that the mean is \( z \).
+All measurements that lead to each individual \( x_i \) are expected to
+be independent, which in turn means that we can express \( \tilde{p} \) as the
+product of individual \( p(x_i) \). The independence assumption is important in the derivation of the central limit theorem.
+
+If we use the integral expression for the \( \delta \)-function + +
+$$
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
+$$
+
+
+and inserting \( e^{i\mu q-i\mu q} \) where \( \mu \) is the mean value
+we arrive at
+
+$$
+ \tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
+ dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
+$$
+
+
+with the integral over \( x \) resulting in
+
+
+$$
+ \int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
+ \int_{-\infty}^{\infty}dxp(x)
+ \left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
+$$
+
+
+The second term on the rhs disappears since this is just the mean and +employing the definition of \( \sigma^2 \) we have +
+$$
+ \int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
+ 1-\frac{q^2\sigma^2}{2m^2}+\dots,
+$$
+
+
+resulting in
+
+
+$$
+ \left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
+ \left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
+$$
+
+
+and in the limit \( m\rightarrow \infty \) we obtain
+
+
+$$
+ \tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
+ \exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
+$$
+
+
+which is the normal distribution with variance
+\( \sigma^2_m=\sigma^2/m \), where \( \sigma \) is the variance of the PDF \( p(x) \)
+and \( \mu \) is also the mean of the PDF \( p(x) \).
+
+Thus, the central limit theorem states that the PDF \( \tilde{p}(z) \) of +the average of \( m \) random values corresponding to a PDF \( p(x) \) +is a normal distribution whose mean is the +mean value of the PDF \( p(x) \) and whose variance is the variance +of the PDF \( p(x) \) divided by \( m \), the number of values used to compute \( z \). + +
+The theorem is satisfied by a large class of PDFs. Note however that for a +finite \( m \), it is not always possible to find a closed expression for +\( \tilde{p}(x) \). +The central limit theorem leads then to the well-known expression for the +standard deviation, given by + +
+$$
+ \sigma_m=
+\frac{\sigma}{\sqrt{m}}.
+$$
+
+
+
+The latter is true only if the average value is known exactly. This is obtained in the limit +\( m\rightarrow \infty \) only. Because the mean and the variance are measured quantities we obtain +the familiar expression in statistics +
+$$
+ \sigma_m\approx
+\frac{\sigma}{\sqrt{m-1}}.
+$$
+
+
+
+In many cases however the above estimate for the standard deviation, +in particular if correlations are strong, may be too simplistic. +
+We see that our new variance and from that the standard deviation, agrees with the central limit theorem.
+
+
+
+
+
+
+
+Suppose we have a PDF \( p(x) \) from which we generate a series \( N \)
+of averages \( \langle x_i \rangle \). Each mean value \( \langle x_i \rangle \)
+is viewed as the average of a specific measurement, e.g., throwing
+dice 100 times and then taking the average value, or producing a certain
+amount of random numbers.
+For notational ease, we set \( \langle x_i \rangle=x_i \) in the discussion
+which follows.
+
+
+If we compute the mean \( z \) of \( m \) such mean values \( x_i \)
+$$
+ z=\frac{x_1+x_2+\dots+x_m}{m},
+$$
+
+the question we pose is which is the PDF of the new variable \( z \).
+
+
+
+The probability of obtaining an average value \( z \) is the product of the
+probabilities of obtaining arbitrary individual mean values \( x_i \),
+but with the constraint that the average is \( z \). We can express this through
+the following expression
+$$
+ \tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
+$$
+
+where the \( \delta \)-function enbodies the constraint that the mean is \( z \).
+All measurements that lead to each individual \( x_i \) are expected to
+be independent, which in turn means that we can express \( \tilde{p} \) as the
+product of individual \( p(x_i) \). The independence assumption is important in the derivation of the central limit theorem.
+
+
+
+If we use the integral expression for the \( \delta \)-function
+
+$$
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
+$$
+
+and inserting \( e^{i\mu q-i\mu q} \) where \( \mu \) is the mean value
+we arrive at
+$$
+ \tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
+ dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
+$$
+
+with the integral over \( x \) resulting in
+
+$$
+ \int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
+ \int_{-\infty}^{\infty}dxp(x)
+ \left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
+$$
+
+
+
+The second term on the rhs disappears since this is just the mean and
+employing the definition of \( \sigma^2 \) we have
+$$
+ \int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
+ 1-\frac{q^2\sigma^2}{2m^2}+\dots,
+$$
+
+resulting in
+
+$$
+ \left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
+ \left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
+$$
+
+and in the limit \( m\rightarrow \infty \) we obtain
+
+$$
+ \tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
+ \exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
+$$
+
+which is the normal distribution with variance
+\( \sigma^2_m=\sigma^2/m \), where \( \sigma \) is the variance of the PDF \( p(x) \)
+and \( \mu \) is also the mean of the PDF \( p(x) \).
+
+
+
+Thus, the central limit theorem states that the PDF \( \tilde{p}(z) \) of
+the average of \( m \) random values corresponding to a PDF \( p(x) \)
+is a normal distribution whose mean is the
+mean value of the PDF \( p(x) \) and whose variance is the variance
+of the PDF \( p(x) \) divided by \( m \), the number of values used to compute \( z \).
+
+
+The theorem is satisfied by a large class of PDFs. Note however that for a
+finite \( m \), it is not always possible to find a closed expression for
+\( \tilde{p}(x) \).
+The central limit theorem leads then to the well-known expression for the
+standard deviation, given by
+
+$$
+ \sigma_m=
+\frac{\sigma}{\sqrt{m}}.
+$$
+
+
+The latter is true only if the average value is known exactly. This is obtained in the limit
+\( m\rightarrow \infty \) only. Because the mean and the variance are measured quantities we obtain
+the familiar expression in statistics
+$$
+ \sigma_m\approx
+\frac{\sigma}{\sqrt{m-1}}.
+$$
+
+
+In many cases however the above estimate for the standard deviation,
+in particular if correlations are strong, may be too simplistic.
+
+
+
+
@@ -1024,10 +1181,12 @@ theorem.
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
-# Returns mean of bootstrap samples
-def stat(data):
+# Returns mean of bootstrap samples # Alternatively, we can run it using Scikit-Learn's function resample # See the examples below
+
+def statistics(data):
return mean(data)
+
# Bootstrap algorithm
def bootstrap(data, statistic, R):
t = zeros(R); n = len(data); inds = arange(n); t0 = time()
@@ -1046,8 +1205,19 @@ mu, sigma = 100, 10000
x = mu + sigma*random.randn(datapoints)
# bootstrap returns the data sample
-t = bootstrap(x, stat, datapoints)
-# the histogram of the bootstrapped data
+t = bootstrap(x, statistics, datapoints)
+
+We see that our new variance and from that the standard deviation, agrees with the central limit theorem.
+
+
+
+
+
+
+Suppose we have a PDF \( p(x) \) from which we generate a series \( N \)
+of averages \( \langle x_i \rangle \). Each mean value \( \langle x_i \rangle \)
+is viewed as the average of a specific measurement, e.g., throwing
+dice 100 times and then taking the average value, or producing a certain
+amount of random numbers.
+For notational ease, we set \( \langle x_i \rangle=x_i \) in the discussion
+which follows.
+
+
+If we compute the mean \( z \) of \( m \) such mean values \( x_i \)
+$$
+ z=\frac{x_1+x_2+\dots+x_m}{m},
+$$
+
+the question we pose is which is the PDF of the new variable \( z \).
+
+
+
+The probability of obtaining an average value \( z \) is the product of the
+probabilities of obtaining arbitrary individual mean values \( x_i \),
+but with the constraint that the average is \( z \). We can express this through
+the following expression
+$$
+ \tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
+$$
+
+where the \( \delta \)-function enbodies the constraint that the mean is \( z \).
+All measurements that lead to each individual \( x_i \) are expected to
+be independent, which in turn means that we can express \( \tilde{p} \) as the
+product of individual \( p(x_i) \). The independence assumption is important in the derivation of the central limit theorem.
+
+
+
+If we use the integral expression for the \( \delta \)-function
+
+$$
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
+$$
+
+and inserting \( e^{i\mu q-i\mu q} \) where \( \mu \) is the mean value
+we arrive at
+$$
+ \tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
+ dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
+$$
+
+with the integral over \( x \) resulting in
+
+$$
+ \int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
+ \int_{-\infty}^{\infty}dxp(x)
+ \left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
+$$
+
+
+
+The second term on the rhs disappears since this is just the mean and
+employing the definition of \( \sigma^2 \) we have
+$$
+ \int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
+ 1-\frac{q^2\sigma^2}{2m^2}+\dots,
+$$
+
+resulting in
+
+$$
+ \left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
+ \left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
+$$
+
+and in the limit \( m\rightarrow \infty \) we obtain
+
+$$
+ \tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
+ \exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
+$$
+
+which is the normal distribution with variance
+\( \sigma^2_m=\sigma^2/m \), where \( \sigma \) is the variance of the PDF \( p(x) \)
+and \( \mu \) is also the mean of the PDF \( p(x) \).
+
+
+
+Thus, the central limit theorem states that the PDF \( \tilde{p}(z) \) of
+the average of \( m \) random values corresponding to a PDF \( p(x) \)
+is a normal distribution whose mean is the
+mean value of the PDF \( p(x) \) and whose variance is the variance
+of the PDF \( p(x) \) divided by \( m \), the number of values used to compute \( z \).
+
+
+The theorem is satisfied by a large class of PDFs. Note however that for a
+finite \( m \), it is not always possible to find a closed expression for
+\( \tilde{p}(x) \).
+The central limit theorem leads then to the well-known expression for the
+standard deviation, given by
+
+$$
+ \sigma_m=
+\frac{\sigma}{\sqrt{m}}.
+$$
+
+
+The latter is true only if the average value is known exactly. This is obtained in the limit
+\( m\rightarrow \infty \) only. Because the mean and the variance are measured quantities we obtain
+the familiar expression in statistics
+$$
+ \sigma_m\approx
+\frac{\sigma}{\sqrt{m-1}}.
+$$
+
+
+In many cases however the above estimate for the standard deviation,
+in particular if correlations are strong, may be too simplistic.
+
+
+
+
@@ -1029,10 +1186,12 @@ theorem.
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
-# Returns mean of bootstrap samples
-def stat(data):
+# Returns mean of bootstrap samples # Alternatively, we can run it using Scikit-Learn's function resample # See the examples below
+
+def statistics(data):
return mean(data)
+
# Bootstrap algorithm
def bootstrap(data, statistic, R):
t = zeros(R); n = len(data); inds = arange(n); t0 = time()
@@ -1051,8 +1210,19 @@ mu, sigma = 100
datapoints = 10000
x = mu + sigma*random.randn(datapoints)
# bootstrap returns the data sample
-t = bootstrap(x, stat, datapoints)
-# the histogram of the bootstrapped data
+t = bootstrap(x, statistics, datapoints)
+
+We see that our new variance and from that the standard deviation, agrees with the central limit theorem.
+
+
+
+
+
+Plotting the Histogram
+# the histogram of the bootstrapped data
n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
# add a 'best fit' line
diff --git a/doc/pub/week37/html/week37-solarized.html b/doc/pub/week37/html/week37-solarized.html
index 1a7823d9f..baf9dd49b 100644
--- a/doc/pub/week37/html/week37-solarized.html
+++ b/doc/pub/week37/html/week37-solarized.html
@@ -115,6 +115,18 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'resampling-methods-bootstrap'),
+ ('The Central Limit Theorem',
+ 2,
+ None,
+ 'the-central-limit-theorem'),
+ ('Finding the Limit', 2, None, 'finding-the-limit'),
+ ('Rewriting the $\\delta$-function',
+ 2,
+ None,
+ 'rewriting-the-delta-function'),
+ ('Identifying Terms', 2, None, 'identifying-terms'),
+ ('Wrapping it up', 2, None, 'wrapping-it-up'),
+ ('Confidence Intervals', 2, None, 'confidence-intervals'),
('Resampling methods: Bootstrap background',
2,
None,
@@ -135,6 +147,7 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'code-example-for-the-bootstrap-method'),
+ ('Plotting the Histogram', 2, None, 'plotting-the-histogram'),
('The bias-variance tradeoff',
2,
None,
@@ -918,6 +931,150 @@ Before we proceed however, we need to remind ourselves about a central theorem i
+The Central Limit Theorem
+
+
+
+Finding the Limit
+
+
+
+Rewriting the \( \delta \)-function
+
+
+
+Identifying Terms
+
+
+
+Wrapping it up
+
+
+
+Confidence Intervals
+
+
+
Resampling methods: Bootstrap background
+
+Plotting the Histogram
+# the histogram of the bootstrapped data
n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
# add a 'best fit' line
diff --git a/doc/pub/week37/html/week37.html b/doc/pub/week37/html/week37.html
index 55b22f101..4155098ea 100644
--- a/doc/pub/week37/html/week37.html
+++ b/doc/pub/week37/html/week37.html
@@ -120,6 +120,18 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'resampling-methods-bootstrap'),
+ ('The Central Limit Theorem',
+ 2,
+ None,
+ 'the-central-limit-theorem'),
+ ('Finding the Limit', 2, None, 'finding-the-limit'),
+ ('Rewriting the $\\delta$-function',
+ 2,
+ None,
+ 'rewriting-the-delta-function'),
+ ('Identifying Terms', 2, None, 'identifying-terms'),
+ ('Wrapping it up', 2, None, 'wrapping-it-up'),
+ ('Confidence Intervals', 2, None, 'confidence-intervals'),
('Resampling methods: Bootstrap background',
2,
None,
@@ -140,6 +152,7 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'code-example-for-the-bootstrap-method'),
+ ('Plotting the Histogram', 2, None, 'plotting-the-histogram'),
('The bias-variance tradeoff',
2,
None,
@@ -923,6 +936,150 @@ Before we proceed however, we need to remind ourselves about a central theorem i
+The Central Limit Theorem
+
+
+
+Finding the Limit
+
+
+
+Rewriting the \( \delta \)-function
+
+
+
+Identifying Terms
+
+
+
+Wrapping it up
+
+
+
+Confidence Intervals
+
+
+
Resampling methods: Bootstrap background
+
+Plotting the Histogram
+# the histogram of the bootstrapped data
n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
# add a 'best fit' line
diff --git a/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz b/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz
index 788a81b6a..bc36c717d 100644
Binary files a/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz and b/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz differ
diff --git a/doc/pub/week37/ipynb/week37.ipynb b/doc/pub/week37/ipynb/week37.ipynb
index ba7ad38df..e18b4bcde 100644
--- a/doc/pub/week37/ipynb/week37.ipynb
+++ b/doc/pub/week37/ipynb/week37.ipynb
@@ -881,6 +881,231 @@
"\n",
"Before we proceed however, we need to remind ourselves about a central theorem in statistics, namely the so-called **central limit theorem**.\n",
"\n",
+ "## The Central Limit Theorem\n",
+ "\n",
+ "\n",
+ "Suppose we have a PDF $p(x)$ from which we generate a series $N$\n",
+ "of averages $\\langle x_i \\rangle$. Each mean value $\\langle x_i \\rangle$\n",
+ "is viewed as the average of a specific measurement, e.g., throwing \n",
+ "dice 100 times and then taking the average value, or producing a certain\n",
+ "amount of random numbers. \n",
+ "For notational ease, we set $\\langle x_i \\rangle=x_i$ in the discussion\n",
+ "which follows. \n",
+ "\n",
+ "If we compute the mean $z$ of $m$ such mean values $x_i$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "z=\\frac{x_1+x_2+\\dots+x_m}{m},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "the question we pose is which is the PDF of the new variable $z$.\n",
+ "\n",
+ "## Finding the Limit\n",
+ "\n",
+ "The probability of obtaining an average value $z$ is the product of the \n",
+ "probabilities of obtaining arbitrary individual mean values $x_i$,\n",
+ "but with the constraint that the average is $z$. We can express this through\n",
+ "the following expression"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\tilde{p}(z)=\\int dx_1p(x_1)\\int dx_2p(x_2)\\dots\\int dx_mp(x_m)\n",
+ " \\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m}),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where the $\\delta$-function enbodies the constraint that the mean is $z$.\n",
+ "All measurements that lead to each individual $x_i$ are expected to\n",
+ "be independent, which in turn means that we can express $\\tilde{p}$ as the \n",
+ "product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.\n",
+ "\n",
+ "\n",
+ "## Rewriting the $\\delta$-function\n",
+ "\n",
+ "If we use the integral expression for the $\\delta$-function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m})=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
+ " dq\\exp{\\left(iq(z-\\frac{x_1+x_2+\\dots+x_m}{m})\\right)},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and inserting $e^{i\\mu q-i\\mu q}$ where $\\mu$ is the mean value\n",
+ "we arrive at"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\tilde{p}(z)=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
+ " dq\\exp{\\left(iq(z-\\mu)\\right)}\\left[\\int_{-\\infty}^{\\infty}\n",
+ " dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "with the integral over $x$ resulting in"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}=\n",
+ " \\int_{-\\infty}^{\\infty}dxp(x)\n",
+ " \\left[1+\\frac{iq(\\mu-x)}{m}-\\frac{q^2(\\mu-x)^2}{2m^2}+\\dots\\right].\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Identifying Terms\n",
+ "\n",
+ "The second term on the rhs disappears since this is just the mean and \n",
+ "employing the definition of $\\sigma^2$ we have"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\int_{-\\infty}^{\\infty}dxp(x)e^{\\left(iq(\\mu-x)/m\\right)}=\n",
+ " 1-\\frac{q^2\\sigma^2}{2m^2}+\\dots,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "resulting in"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\left[\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m\\approx\n",
+ " \\left[1-\\frac{q^2\\sigma^2}{2m^2}+\\dots \\right]^m,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and in the limit $m\\rightarrow \\infty$ we obtain"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\tilde{p}(z)=\\frac{1}{\\sqrt{2\\pi}(\\sigma/\\sqrt{m})}\n",
+ " \\exp{\\left(-\\frac{(z-\\mu)^2}{2(\\sigma/\\sqrt{m})^2}\\right)},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which is the normal distribution with variance\n",
+ "$\\sigma^2_m=\\sigma^2/m$, where $\\sigma$ is the variance of the PDF $p(x)$\n",
+ "and $\\mu$ is also the mean of the PDF $p(x)$. \n",
+ "\n",
+ "## Wrapping it up\n",
+ "\n",
+ "Thus, the central limit theorem states that the PDF $\\tilde{p}(z)$ of\n",
+ "the average of $m$ random values corresponding to a PDF $p(x)$ \n",
+ "is a normal distribution whose mean is the \n",
+ "mean value of the PDF $p(x)$ and whose variance is the variance\n",
+ "of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.\n",
+ "\n",
+ "The theorem is satisfied by a large class of PDFs. Note however that for a\n",
+ "finite $m$, it is not always possible to find a closed expression for\n",
+ "$\\tilde{p}(x)$.\n",
+ "The central limit theorem leads then to the well-known expression for the\n",
+ "standard deviation, given by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\sigma_m=\n",
+ "\\frac{\\sigma}{\\sqrt{m}}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The latter is true only if the average value is known exactly. This is obtained in the limit\n",
+ "$m\\rightarrow \\infty$ only. Because the mean and the variance are measured quantities we obtain \n",
+ "the familiar expression in statistics"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\sigma_m\\approx \n",
+ "\\frac{\\sigma}{\\sqrt{m-1}}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In many cases however the above estimate for the standard deviation,\n",
+ "in particular if correlations are strong, may be too simplistic.\n",
+ "\n",
+ "## Confidence Intervals\n",
+ "\n",
+ "\n",
"## Resampling methods: Bootstrap background\n",
"\n",
"Since $\\widehat{\\beta} = \\widehat{\\beta}(\\boldsymbol{X})$ is a function of random variables,\n",
@@ -979,10 +1204,12 @@
"import matplotlib.mlab as mlab\n",
"import matplotlib.pyplot as plt\n",
"\n",
- "# Returns mean of bootstrap samples \n",
- "def stat(data):\n",
+ "# Returns mean of bootstrap samples # Alternatively, we can run it using Scikit-Learn's function resample # See the examples below\n",
+ "\n",
+ "def statistics(data):\n",
" return mean(data)\n",
"\n",
+ "\n",
"# Bootstrap algorithm\n",
"def bootstrap(data, statistic, R):\n",
" t = zeros(R); n = len(data); inds = arange(n); t0 = time()\n",
@@ -1001,7 +1228,27 @@
"datapoints = 10000\n",
"x = mu + sigma*random.randn(datapoints)\n",
"# bootstrap returns the data sample \n",
- "t = bootstrap(x, stat, datapoints)\n",
+ "t = bootstrap(x, statistics, datapoints)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We see that our new variance and from that the standard deviation, agrees with the central limit theorem.\n",
+ "\n",
+ "## Plotting the Histogram"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "editable": true
+ },
+ "outputs": [],
+ "source": [
"# the histogram of the bootstrapped data \n",
"n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)\n",
"\n",
diff --git a/doc/src/week37/test.py b/doc/src/week37/test.py
deleted file mode 100644
index 958bceb5f..000000000
--- a/doc/src/week37/test.py
+++ /dev/null
@@ -1,38 +0,0 @@
-import numpy as np
-import pandas as pd
-from IPython.display import display
-import matplotlib.pyplot as plt
-from sklearn.model_selection import train_test_split
-from sklearn import linear_model
-
-# Make data set.
-n = 1000
-x = np.random.rand(n)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)#+ np.random.randn(n)
-
-Maxpolydegree = 5
-X = np.zeros((len(x),Maxpolydegree))
-X[:,0] = 1.0
-
-for polydegree in range(1, Maxpolydegree):
- for degree in range(polydegree):
- X[:,degree] = x**(degree)
-
-
-# We split the data in test and training data
-X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
-
-# Decide which values of lambda to use
-nlambdas = 2
-lambdas = np.logspace(-3, -1, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- # Make the fit using Ridge only
- RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
- RegRidge.fit(X_train,y_train)
- # and then make the prediction
- ypredictRidge = RegRidge.predict(X_test)
- Coeffs = np.array(RegRidge.coef_)
- BetaValues = pd.DataFrame(Coeffs)
- BetaValues.columns = ['beta']
- display(BetaValues)
diff --git a/doc/src/week37/week37.do.txt b/doc/src/week37/week37.do.txt
index 74f5eae2f..2b8043666 100644
--- a/doc/src/week37/week37.do.txt
+++ b/doc/src/week37/week37.do.txt
@@ -636,6 +636,145 @@ The textbook by "Davison on the Bootstrap Methods and their Applications":"https
Before we proceed however, we need to remind ourselves about a central theorem in statistics, namely the so-called _central limit theorem_.
+!split
+===== The Central Limit Theorem =====
+
+
+Suppose we have a PDF $p(x)$ from which we generate a series $N$
+of averages $\langle x_i \rangle$. Each mean value $\langle x_i \rangle$
+is viewed as the average of a specific measurement, e.g., throwing
+dice 100 times and then taking the average value, or producing a certain
+amount of random numbers.
+For notational ease, we set $\langle x_i \rangle=x_i$ in the discussion
+which follows.
+
+If we compute the mean $z$ of $m$ such mean values $x_i$
+!bt
+\[
+ z=\frac{x_1+x_2+\dots+x_m}{m},
+\]
+!et
+the question we pose is which is the PDF of the new variable $z$.
+
+!split
+===== Finding the Limit =====
+
+The probability of obtaining an average value $z$ is the product of the
+probabilities of obtaining arbitrary individual mean values $x_i$,
+but with the constraint that the average is $z$. We can express this through
+the following expression
+!bt
+\[
+ \tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
+\]
+!et
+where the $\delta$-function enbodies the constraint that the mean is $z$.
+All measurements that lead to each individual $x_i$ are expected to
+be independent, which in turn means that we can express $\tilde{p}$ as the
+product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.
+
+
+!split
+===== Rewriting the $\delta$-function =====
+
+If we use the integral expression for the $\delta$-function
+
+!bt
+\[
+ \delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
+\]
+!et
+and inserting $e^{i\mu q-i\mu q}$ where $\mu$ is the mean value
+we arrive at
+!bt
+\[
+ \tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
+ dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
+ dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
+\]
+!et
+with the integral over $x$ resulting in
+
+!bt
+\[
+ \int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
+ \int_{-\infty}^{\infty}dxp(x)
+ \left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
+\]
+!et
+
+!split
+===== Identifying Terms =====
+
+The second term on the rhs disappears since this is just the mean and
+employing the definition of $\sigma^2$ we have
+!bt
+\[
+ \int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
+ 1-\frac{q^2\sigma^2}{2m^2}+\dots,
+\]
+!et
+resulting in
+
+!bt
+\[
+ \left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
+ \left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
+\]
+!et
+and in the limit $m\rightarrow \infty$ we obtain
+
+!bt
+\[
+ \tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
+ \exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
+\]
+!et
+which is the normal distribution with variance
+$\sigma^2_m=\sigma^2/m$, where $\sigma$ is the variance of the PDF $p(x)$
+and $\mu$ is also the mean of the PDF $p(x)$.
+
+!split
+===== Wrapping it up =====
+
+Thus, the central limit theorem states that the PDF $\tilde{p}(z)$ of
+the average of $m$ random values corresponding to a PDF $p(x)$
+is a normal distribution whose mean is the
+mean value of the PDF $p(x)$ and whose variance is the variance
+of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.
+
+The theorem is satisfied by a large class of PDFs. Note however that for a
+finite $m$, it is not always possible to find a closed expression for
+$\tilde{p}(x)$.
+The central limit theorem leads then to the well-known expression for the
+standard deviation, given by
+
+!bt
+\[
+ \sigma_m=
+\frac{\sigma}{\sqrt{m}}.
+\]
+!et
+
+The latter is true only if the average value is known exactly. This is obtained in the limit
+$m\rightarrow \infty$ only. Because the mean and the variance are measured quantities we obtain
+the familiar expression in statistics
+!bt
+\[
+ \sigma_m\approx
+\frac{\sigma}{\sqrt{m-1}}.
+\]
+!et
+
+In many cases however the above estimate for the standard deviation,
+in particular if correlations are strong, may be too simplistic.
+
+!split
+===== Confidence Intervals =====
+
+
!split
===== Resampling methods: Bootstrap background =====
@@ -727,10 +866,12 @@ from time import time
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
-# Returns mean of bootstrap samples
-def stat(data):
+# Returns mean of bootstrap samples # Alternatively, we can run it using Scikit-Learn's function resample # See the examples below
+
+def statistics(data):
return mean(data)
+
# Bootstrap algorithm
def bootstrap(data, statistic, R):
t = zeros(R); n = len(data); inds = arange(n); t0 = time()
@@ -749,7 +890,13 @@ mu, sigma = 100, 15
datapoints = 10000
x = mu + sigma*random.randn(datapoints)
# bootstrap returns the data sample
-t = bootstrap(x, stat, datapoints)
+t = bootstrap(x, statistics, datapoints)
+!ec
+We see that our new variance and from that the standard deviation, agrees with the central limit theorem.
+
+!split
+===== Plotting the Histogram =====
+!bc pycod
# the histogram of the bootstrapped data
n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)