import numpy as np

# set a seed for reproducibility
np.random.seed(1234)

# generate 1000 random numbers (between 0 and 1) for each model, representing 1000 observations
mod1 = np.random.rand(1000)
mod2 = np.random.rand(1000)
mod3 = np.random.rand(1000)
mod4 = np.random.rand(1000)
mod5 = np.random.rand(1000)

# each model independently predicts 1 (the "correct response") if random number was at least 0.3
preds1 = np.where(mod1 > 0.3, 1, 0)
preds2 = np.where(mod2 > 0.3, 1, 0)
preds3 = np.where(mod3 > 0.3, 1, 0)
preds4 = np.where(mod4 > 0.3, 1, 0)
preds5 = np.where(mod5 > 0.3, 1, 0)

# print the first 20 predictions from each model
print(preds1[:20])
print(preds2[:20])
print(preds3[:20])
print(preds4[:20])
print(preds5[:20])

# average the predictions and then round to 0 or 1
ensemble_preds = np.round((preds1 + preds2 + preds3 + preds4 + preds5)/5.0).astype(int)

# print the ensemble's first 20 predictions
print(ensemble_preds[:20])

# how accurate was each individual model?
print(preds1.mean())
print(preds2.mean())
print(preds3.mean())
print(preds4.mean())
print(preds5.mean())

# how accurate was the ensemble?
print(ensemble_preds.mean())

"""
Explanation: 08 - Ensemble Methods - Bagging
by Alejandro Correa Bahnsen & Iván Torroledo
version 1.2, Feb 2018
Part of the class Machine Learning for Risk Management
This notebook is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License. Special thanks goes to Kevin Markham
Why are we learning about ensembling?

Very popular method for improving the predictive performance of machine learning models
Provides a foundation for understanding more sophisticated models

Lesson objectives
Students will be able to:

Define ensembling and its requirements
Identify the two basic methods of ensembling
Decide whether manual ensembling is a useful approach for a given problem
Explain bagging and how it can be applied to decision trees
Explain how out-of-bag error and feature importances are calculated from bagged trees
Explain the difference between bagged trees and Random Forests
Build and tune a Random Forest model in scikit-learn
Decide whether a decision tree or a Random Forest is a better model for a given problem

Part 1: Introduction
Ensemble learning is a widely studied topic in the machine learning community. The main idea behind 
the ensemble methodology is to combine several individual base classifiers in   order to have a 
classifier that outperforms each of them.
Nowadays,   ensemble methods are  one 
of the most popular and well studied machine learning techniques, and it can be 
noted that since 2009 all the first-place and   second-place winners of the KDD-Cup https://www.sigkdd.org/kddcup/   used  ensemble methods. The core 
principle in ensemble learning, is to induce random perturbations into  the learning procedure in 
order to produce several different base classifiers from a single  training set, then combining the 
base classifiers in order to make the final prediction.  In order to induce the random permutations 
and therefore create the different base classifiers,   several methods have been proposed, in 
particular: 
* bagging
* pasting
* random forests 
* random patches  
Finally, after  the base   classifiers 
are trained, they are typically   combined using either:
* majority voting
* weighted  voting
* stacking
There are three main reasons regarding why ensemble 
methods perform better than single models: statistical, computational and representational . First, from a statistical point of view, when the learning set is too 
small, an algorithm can find several good models within the search space, that arise to the same 
performance on the training set $\mathcal{S}$. Nevertheless, without a validation set, there is 
a risk of choosing the wrong model. The second reason is computational; in general, algorithms 
rely on some local search optimization and may get stuck in a local optima. Then, an ensemble may 
solve this by focusing different algorithms to different spaces across the training set. The last 
reason is representational. In most cases, for a learning set of finite size, the  true function 
$f$ cannot be represented by any of the candidate models. By combining several  models in an 
ensemble, it may be possible to obtain a model with a larger coverage across the  space of 
representable functions.

Example
Let's pretend that instead of building a single model to solve a binary classification problem, you created five independent models, and each model was correct about 70% of the time. If you combined these models into an "ensemble" and used their majority vote as a prediction, how often would the ensemble be correct?
End of explanation
"""


# read in and prepare the vehicle training data
import zipfile
import pandas as pd
with zipfile.ZipFile('../datasets/vehicles_train.csv.zip', 'r') as z:
    f = z.open('vehicles_train.csv')
    train = pd.io.parsers.read_table(f, index_col=False, sep=',')
with zipfile.ZipFile('../datasets/vehicles_test.csv.zip', 'r') as z:
    f = z.open('vehicles_test.csv')
    test = pd.io.parsers.read_table(f, index_col=False, sep=',')

train['vtype'] = train.vtype.map({'car':0, 'truck':1})
# read in and prepare the vehicle testing data
test['vtype'] = test.vtype.map({'car':0, 'truck':1})

train.head()

"""
Explanation: Note: As you add more models to the voting process, the probability of error decreases, which is known as Condorcet's Jury Theorem.
What is ensembling?
Ensemble learning (or "ensembling") is the process of combining several predictive models in order to produce a combined model that is more accurate than any individual model.

Regression: take the average of the predictions
Classification: take a vote and use the most common prediction, or take the average of the predicted probabilities

For ensembling to work well, the models must have the following characteristics:

Accurate: they outperform the null model
Independent: their predictions are generated using different processes

The big idea: If you have a collection of individually imperfect (and independent) models, the "one-off" mistakes made by each model are probably not going to be made by the rest of the models, and thus the mistakes will be discarded when averaging the models.
There are two basic methods for ensembling:

Manually ensemble your individual models
Use a model that ensembles for you

Theoretical performance of an ensemble
If we assume that each one of the $T$ base classifiers has a probability $\rho$ of 
  being correct, the probability of an ensemble making the correct decision, assuming independence, 
  denoted by $P_c$, can be calculated using the binomial distribution
$$P_c = \sum_{j>T/2}^{T} {{T}\choose{j}} \rho^j(1-\rho)^{T-j}.$$
Furthermore, as shown, if $T\ge3$ then:
$$
  \lim_{T \to  \infty} P_c= \begin{cases} 
            1  &\mbox{if } \rho>0.5 \ 
            0  &\mbox{if } \rho<0.5 \ 
            0.5  &\mbox{if } \rho=0.5 ,
            \end{cases}
$$
    leading to the conclusion that 
$$
  \rho \ge 0.5 \quad \text{and} \quad T\ge3 \quad \Rightarrow \quad P_c\ge \rho.
$$
Part 2: Manual ensembling
What makes a good manual ensemble?

Different types of models
Different combinations of features
Different tuning parameters


Machine learning flowchart created by the winner of Kaggle's CrowdFlower competition
End of explanation
"""


from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsRegressor

models = {'lr': LinearRegression(),
          'dt': DecisionTreeRegressor(),
          'nb': GaussianNB(),
          'nn': KNeighborsRegressor()}

# Train all the models
X_train = train.iloc[:, 1:]
X_test = test.iloc[:, 1:]
y_train = train.price
y_test = test.price

for model in models.keys():
    models[model].fit(X_train, y_train)

# predict test for each model
y_pred = pd.DataFrame(index=test.index, columns=models.keys())
for model in models.keys():
    y_pred[model] = models[model].predict(X_test)
    

# Evaluate each model
from sklearn.metrics import mean_squared_error

for model in models.keys():
    print(model,np.sqrt(mean_squared_error(y_pred[model], y_test)))

"""
Explanation: Train different models
End of explanation
"""


np.sqrt(mean_squared_error(y_pred.mean(axis=1), y_test))

"""
Explanation: Evaluate the error of the mean of the predictions
End of explanation
"""


# set a seed for reproducibility
np.random.seed(1)

# create an array of 1 through 20
nums = np.arange(1, 21)
print(nums)

# sample that array 20 times with replacement
print(np.random.choice(a=nums, size=20, replace=True))

"""
Explanation: Comparing manual ensembling with a single model approach
Advantages of manual ensembling:

Increases predictive accuracy
Easy to get started

Disadvantages of manual ensembling:

Decreases interpretability
Takes longer to train
Takes longer to predict
More complex to automate and maintain
Small gains in accuracy may not be worth the added complexity

Part 3: Bagging
The primary weakness of decision trees is that they don't tend to have the best predictive accuracy. This is partially due to high variance, meaning that different splits in the training data can lead to very different trees.
Bagging is a general purpose procedure for reducing the variance of a machine learning method, but is particularly useful for decision trees. Bagging is short for bootstrap aggregation, meaning the aggregation of bootstrap samples.
What is a bootstrap sample? A random sample with replacement:
End of explanation
"""


# set a seed for reproducibility
np.random.seed(123)

n_samples = train.shape[0]
n_B = 10

# create ten bootstrap samples (will be used to select rows from the DataFrame)
samples = [np.random.choice(a=n_samples, size=n_samples, replace=True) for _ in range(1, n_B +1 )]
samples

# show the rows for the first decision tree
train.iloc[samples[0], :]

"""
Explanation: How does bagging work (for decision trees)?

Grow B trees using B bootstrap samples from the training data.
Train each tree on its bootstrap sample and make predictions.
Combine the predictions:
Average the predictions for regression trees
Take a vote for classification trees



Notes:

Each bootstrap sample should be the same size as the original training set.
B should be a large enough value that the error seems to have "stabilized".
The trees are grown deep so that they have low bias/high variance.

Bagging increases predictive accuracy by reducing the variance, similar to how cross-validation reduces the variance associated with train/test split (for estimating out-of-sample error) by splitting many times an averaging the results.
End of explanation
"""


from sklearn.tree import DecisionTreeRegressor

# grow each tree deep
treereg = DecisionTreeRegressor(max_depth=None, random_state=123)

# DataFrame for storing predicted price from each tree
y_pred = pd.DataFrame(index=test.index, columns=[list(range(n_B))])

# grow one tree for each bootstrap sample and make predictions on testing data
for i, sample in enumerate(samples):
    X_train = train.iloc[sample, 1:]
    y_train = train.iloc[sample, 0]
    treereg.fit(X_train, y_train)
    y_pred[i] = treereg.predict(X_test)

y_pred

"""
Explanation: Build one tree for each sample
End of explanation
"""


for i in range(n_B):
    print(i, np.sqrt(mean_squared_error(y_pred[i], y_test)))

"""
Explanation: Results of each tree
End of explanation
"""


y_pred.mean(axis=1)

np.sqrt(mean_squared_error(y_test, y_pred.mean(axis=1)))

"""
Explanation: Results of the ensemble
End of explanation
"""


# define the training and testing sets
X_train = train.iloc[:, 1:]
y_train = train.iloc[:, 0]
X_test = test.iloc[:, 1:]
y_test = test.iloc[:, 0]

# instruct BaggingRegressor to use DecisionTreeRegressor as the "base estimator"
from sklearn.ensemble import BaggingRegressor
bagreg = BaggingRegressor(DecisionTreeRegressor(), n_estimators=500, 
                          bootstrap=True, oob_score=True, random_state=1)

# fit and predict
bagreg.fit(X_train, y_train)
y_pred = bagreg.predict(X_test)
y_pred

# calculate RMSE
np.sqrt(mean_squared_error(y_test, y_pred))

"""
Explanation: Bagged decision trees in scikit-learn (with B=500)
End of explanation
"""


# show the first bootstrap sample
samples[0]

# show the "in-bag" observations for each sample
for sample in samples:
    print(set(sample))

# show the "out-of-bag" observations for each sample
for sample in samples:
    print(sorted(set(range(n_samples)) - set(sample)))

"""
Explanation: Estimating out-of-sample error
For bagged models, out-of-sample error can be estimated without using train/test split or cross-validation!
On average, each bagged tree uses about two-thirds of the observations. For each tree, the remaining observations are called "out-of-bag" observations.
End of explanation
"""


# compute the out-of-bag R-squared score (not MSE, unfortunately!) for B=500
bagreg.oob_score_

"""
Explanation: How to calculate "out-of-bag error":

For every observation in the training data, predict its response value using only the trees in which that observation was out-of-bag. Average those predictions (for regression) or take a vote (for classification).
Compare all predictions to the actual response values in order to compute the out-of-bag error.

When B is sufficiently large, the out-of-bag error is an accurate estimate of out-of-sample error.
End of explanation
"""


# read in the data
with zipfile.ZipFile('../datasets/hitters.csv.zip', 'r') as z:
    f = z.open('hitters.csv')
    hitters = pd.read_csv(f, sep=',', index_col=False)

# remove rows with missing values
hitters.dropna(inplace=True)
hitters.head()

# encode categorical variables as integers
hitters['League'] = pd.factorize(hitters.League)[0]
hitters['Division'] = pd.factorize(hitters.Division)[0]
hitters['NewLeague'] = pd.factorize(hitters.NewLeague)[0]
hitters.head()

# allow plots to appear in the notebook
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

# scatter plot of Years versus Hits colored by Salary
hitters.plot(kind='scatter', x='Years', y='Hits', c='Salary', colormap='jet', xlim=(0, 25), ylim=(0, 250))

# define features: exclude career statistics (which start with "C") and the response (Salary)
feature_cols = hitters.columns[hitters.columns.str.startswith('C') == False].drop('Salary')
feature_cols

# define X and y
X = hitters[feature_cols]
y = hitters.Salary

"""
Explanation: Estimating feature importance
Bagging increases predictive accuracy, but decreases model interpretability because it's no longer possible to visualize the tree to understand the importance of each feature.
However, we can still obtain an overall summary of feature importance from bagged models:

Bagged regression trees: calculate the total amount that MSE is decreased due to splits over a given feature, averaged over all trees
Bagged classification trees: calculate the total amount that Gini index is decreased due to splits over a given feature, averaged over all trees

Part 4: Random Forests
Random Forests is a slight variation of bagged trees that has even better performance:

Exactly like bagging, we create an ensemble of decision trees using bootstrapped samples of the training set.
However, when building each tree, each time a split is considered, a random sample of m features is chosen as split candidates from the full set of p features. The split is only allowed to use one of those m features.
A new random sample of features is chosen for every single tree at every single split.
For classification, m is typically chosen to be the square root of p.
For regression, m is typically chosen to be somewhere between p/3 and p.



What's the point?

Suppose there is one very strong feature in the data set. When using bagged trees, most of the trees will use that feature as the top split, resulting in an ensemble of similar trees that are highly correlated.
Averaging highly correlated quantities does not significantly reduce variance (which is the entire goal of bagging).
By randomly leaving out candidate features from each split, Random Forests "decorrelates" the trees, such that the averaging process can reduce the variance of the resulting model.

Part 5: Building and tuning decision trees and Random Forests

Major League Baseball player data from 1986-87: data, data dictionary (page 7)
Each observation represents a player
Goal: Predict player salary
End of explanation
"""


# list of values to try for max_depth
max_depth_range = range(1, 21)

# list to store the average RMSE for each value of max_depth
RMSE_scores = []

# use 10-fold cross-validation with each value of max_depth
from sklearn.model_selection import cross_val_score
for depth in max_depth_range:
    treereg = DecisionTreeRegressor(max_depth=depth, random_state=1)
    MSE_scores = cross_val_score(treereg, X, y, cv=10, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# plot max_depth (x-axis) versus RMSE (y-axis)
plt.plot(max_depth_range, RMSE_scores)
plt.xlabel('max_depth')
plt.ylabel('RMSE (lower is better)')

# show the best RMSE and the corresponding max_depth
sorted(zip(RMSE_scores, max_depth_range))[0]

# max_depth=2 was best, so fit a tree using that parameter
treereg = DecisionTreeRegressor(max_depth=2, random_state=1)
treereg.fit(X, y)

# compute feature importances
pd.DataFrame({'feature':feature_cols, 'importance':treereg.feature_importances_}).sort_values('importance')

"""
Explanation: Predicting salary with a decision tree
Find the best max_depth for a decision tree using cross-validation:
End of explanation
"""


from sklearn.ensemble import RandomForestRegressor
rfreg = RandomForestRegressor()
rfreg

"""
Explanation: Predicting salary with a Random Forest
End of explanation
"""


# list of values to try for n_estimators
estimator_range = range(10, 310, 10)

# list to store the average RMSE for each value of n_estimators
RMSE_scores = []

# use 5-fold cross-validation with each value of n_estimators (WARNING: SLOW!)
for estimator in estimator_range:
    rfreg = RandomForestRegressor(n_estimators=estimator, random_state=1, n_jobs=-1)
    MSE_scores = cross_val_score(rfreg, X, y, cv=5, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# plot n_estimators (x-axis) versus RMSE (y-axis)
plt.plot(estimator_range, RMSE_scores)
plt.xlabel('n_estimators')
plt.ylabel('RMSE (lower is better)')

"""
Explanation: Tuning n_estimators
One important tuning parameter is n_estimators, which is the number of trees that should be grown. It should be a large enough value that the error seems to have "stabilized".
End of explanation
"""


# list of values to try for max_features
feature_range = range(1, len(feature_cols)+1)

# list to store the average RMSE for each value of max_features
RMSE_scores = []

# use 10-fold cross-validation with each value of max_features (WARNING: SLOW!)
for feature in feature_range:
    rfreg = RandomForestRegressor(n_estimators=150, max_features=feature, random_state=1, n_jobs=-1)
    MSE_scores = cross_val_score(rfreg, X, y, cv=10, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# plot max_features (x-axis) versus RMSE (y-axis)
plt.plot(feature_range, RMSE_scores)
plt.xlabel('max_features')
plt.ylabel('RMSE (lower is better)')

# show the best RMSE and the corresponding max_features
sorted(zip(RMSE_scores, feature_range))[0]

"""
Explanation: Tuning max_features
The other important tuning parameter is max_features, which is the number of features that should be considered at each split.
End of explanation
"""


# max_features=8 is best and n_estimators=150 is sufficiently large
rfreg = RandomForestRegressor(n_estimators=150, max_features=8, oob_score=True, random_state=1)
rfreg.fit(X, y)

# compute feature importances
pd.DataFrame({'feature':feature_cols, 'importance':rfreg.feature_importances_}).sort_values('importance')

# compute the out-of-bag R-squared score
rfreg.oob_score_

"""
Explanation: Fitting a Random Forest with the best parameters
End of explanation
"""


# check the shape of X
X.shape

rfreg

# set a threshold for which features to include
from sklearn.feature_selection import SelectFromModel
print(SelectFromModel(rfreg, threshold=0.1, prefit=True).transform(X).shape)
print(SelectFromModel(rfreg, threshold='mean', prefit=True).transform(X).shape)
print(SelectFromModel(rfreg, threshold='median', prefit=True).transform(X).shape)

# create a new feature matrix that only includes important features
X_important = SelectFromModel(rfreg, threshold='mean', prefit=True).transform(X)

# check the RMSE for a Random Forest that only includes important features
rfreg = RandomForestRegressor(n_estimators=150, max_features=3, random_state=1)
scores = cross_val_score(rfreg, X_important, y, cv=10, scoring='neg_mean_squared_error')
np.mean(np.sqrt(-scores))

"""
Explanation: Reducing X to its most important features
End of explanation
"""


