<jupyter_start><jupyter_text>Energy Usage in Saudi Arabia<jupyter_code>import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import os
DATASET_PATH = "dataset"<jupyter_output><empty_output><jupyter_text>Contents  [comment]:  (This is how to do a comment. the  is the syntax for a hyperlink tag in html. The ... is a link to the corresponding tag in markdown)1. [Load and Prepare the Data](load)2. [Explore the Data](explore_the_data)3. [Data Preparation](data_preparation)4. [Train a Model](train_a_model)5. [Evaluate](evaluate)  1. Load and Prepare the Data<jupyter_code>def load_architecture_data(dataset_path=DATASET_PATH):
    csv_path = os.path.join(dataset_path, "data.csv")
    return pd.read_csv(csv_path)
architecture = load_architecture_data()
arch = architecture.drop(columns=['Unnamed: 12', 'Unnamed: 14', 'Unnamed: 15', 'Unnamed: 16', 'Unnamed: 17'])
arch = arch.drop(index=[868, 869, 870, 871, 872, 873, 874])
print(arch.columns)  # columns come out weirdly (especially 'total family members '). Rename them
arch.columns = ['timestamp', 'location', 'family_member', 'age', 'male_members',
       'female_members', 'building_age', 'land_size', 'male_guest_room',
       'female_guest_room', 'living_room', 'outside_room', 'total_family_members']

# for some reason 'outside_room' was of type object and I couldn't change it by editing the csv?
arch['outside_room'] = arch['outside_room'].astype(str).astype(float)
arch.tail(2)
arch.info()
arch.describe()
arch.hist(bins=50, figsize=(20, 15)) # plots value vs frequency density (~ frequency for most of these)
plt.show()<jupyter_output><empty_output><jupyter_text>Things to note about the data:* Guest rooms, land size, living room and building age have small but significant outliers, far larger/more frequent than would be expected from any standard distribution. Most of the distributions have **heavy tails** - defined as a tail probability decaying less quickly than an exponential distribution - i.e. they have high outlier probabilities. This may cause problems for machine learning algorithms, however, it is often easy to transform variables like this to get more bell-shaped distributions.* Location is discrete (and doesnt have any resembelence to where it actually is). Therefore a one-hot encoding should be used.* Living room is the "dodgy data" mentioned earlier. The survey meant the hall, but most people thought it was the lounge* There are no nan/empty data columns, unless you count living room. The solution to this could just be to remove the values at 356. Alternatively as > 90% of people thought this was the lounge, the data could be relabelled (however, this is pointless as there is only one value!)* Building age and land size both have very large values. It would be sensible to normalise these to a more managable scale for a machine learning algorithm.The first thing to do for any machine learning task is to split the data into a training and a test set. Since this dataset is quite small (868 entries) there is a significant amount of sampling bias (e.g. only 1.7% of values in catagory 8, yet there are 5% in reality), therefore we could use stratified sampling. Generally the most *important* attribute would be used for stratifying, however, since I'm not sure what that is, I will use age as age is both the closest to a uniform distribution and a continuous distribution (i.e. it has a good spread over the data).<jupyter_code>arch['age_category'] = pd.cut(arch['age'], 
                              bins = [0.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, np.inf], # min=20.0, max = 89.0
                              labels = [2, 3, 4, 5, 6, 7, 8])
arch['age_category'].hist(bins=6, density=True)
from sklearn.model_selection import train_test_split
train_set, test_set = train_test_split(arch, test_size=0.2, random_state=42)  # seed=42 so test set doesnt change each run
print(len(train_set), len(test_set))
test_set['age_category'].value_counts(sort=False) / len(test_set)
# stratified sampling
from sklearn.model_selection import StratifiedShuffleSplit
split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)

for train_index, test_index in split.split(arch, arch['age_category']):
    strat_train_set = arch.iloc[train_index]
    strat_test_set = arch.iloc[test_index]
strat_test_set['age_category'].value_counts(sort=False) / len(strat_test_set)  # much better!<jupyter_output><empty_output><jupyter_text>Finally, drop the age_category column<jupyter_code>strat_train_set = strat_train_set.drop('age_category', axis=1)
strat_test_set = strat_test_set.drop('age_category', axis=1)<jupyter_output><empty_output><jupyter_text>2. Explore the Data [back to top](top)<jupyter_code>arch = strat_train_set.copy()  # create a copy to work on so we dont mess things up!<jupyter_output><empty_output><jupyter_text>2.1 Correlations<jupyter_code>corr_matrix = arch.corr()
corr_matrix<jupyter_output><empty_output><jupyter_text>Note, the correlation coefficient (pearson's r) only measures linear combinations (and therefore misses non-linear ones). There are a few things to note about this.* The most correlated attributes are `male_members` and `female_members` with `total_family_members`. However, `male_members` are not strongly correlated with `female_members` (odd?)* `female_guest_room` and `male_guest_room` are reasonable correlated.* There is a slight correlation between `age` and `total_family_members`.* Perhaps surprisingly, there is a slight correlation between the number of `female members` and the `location`.* The `living_room`, `outside_room`, `building_age`, `land_size`, `family_member` and `location` are very uncorrelated with everything.We can plot these attributes against eachother to see what is happening more clearly. This would be an 11x11 matrix, generally we would select 3-4 "promising" attributes. However, 11 isn't a massive number so why not just plot them against eachother first to see what the promising attributes are?<jupyter_code>pd.plotting.scatter_matrix(arch, figsize=(17, 12))
plt.show()<jupyter_output><empty_output><jupyter_text>From the above it is obvious that `family_member`, `loaction`, and `building_age` do not reveal anything in their current state. For `location` this is due to the catagorical nature of the data, for `building_age` this is because the data is far too succeptable to outliers and needs to be normalised, and for `family_member` it may be becasue there are only 4 different values. **Not actually sure what `family_member` means. Is it the number that live in the house or the person asked? I.e. Father, Wife, Son, Daughter.**There are a couple more interesting things to glean from these plots. `male_guest_room`, `female_guest_room`,`outside_room` and `living_room` have been filled in by survey participants as guesses. Participants were more likely to round their answer to the nearest 10 than to give an exact number, and therefore it may be a good idea to bin (via pandas `.cut()` method) these attributes to remove this bias towards factors of 10. The attributes with the most promise are `total_family_members`, `male_members` and `female_members` as these have more interesting scatter plots than other attributes in general. We may take a closer look.<jupyter_code>attributes = ['male_members', 'female_members', 'total_family_members']
pd.plotting.scatter_matrix(arch[attributes], figsize=(12, 8))
plt.show()<jupyter_output><empty_output><jupyter_text>`male_members` has almost no correlation with `female_members`. We can repeat this using the target data to see how they are correlated.<jupyter_code>attributes = ['male_guest_room', 'female_guest_room', 'living_room']
pd.plotting.scatter_matrix(arch[attributes], figsize=(12, 8))
plt.show()<jupyter_output><empty_output><jupyter_text>Not much correlation at all. Therefore, it is a good idea to train a different model for each traget variable.  2.2 Combining Attributes It is always a good idea to have a look at sensible combinations of attributes. Note that we can only use combinations that don't include the three y-attributes, which limits our possibilities. Furthermore, location is catagorical, so that can't be used.<jupyter_code>arch['people_per_area'] = arch['total_family_members']/arch['land_size']
arch['males_in_family_ratio'] = arch['male_members']/arch['total_family_members']
arch['mtf_ratio'] = arch['male_members']/arch['female_members']
arch['land_per_building_age'] = arch['land_size']/arch['building_age']
arch['living_to_outside'] = arch['living_room']/arch['outside_room']<jupyter_output><empty_output><jupyter_text>we may have ended up with infinite values if any of the denominators are zero, so replace these with np.nan, and then drop. Now we can deal with the np.nan values. Having a look at the data, the only attributes that don't have a minimum value of 0 are `building_age` and `land_size`. Therefore we need a way to deal with these infinite values. The easiest method (which we will use in the meantime) is just to set them to zero.<jupyter_code>arch.replace([np.inf, -np.inf], 0.0, inplace=True)
# Note we are only interest in finding new correlations in the top right (or bottom left) quadrant.
attributes = ['people_per_area', 'males_in_family_ratio', 'mtf_ratio', 'land_per_building_age', 'living_to_outside',
              'male_guest_room', 'female_guest_room', 'living_room']
arch[attributes].corr()<jupyter_output><empty_output><jupyter_text>Unfortunately, none of these combinations helped too much  3. Data Preparation   [Back to top](top) Data preparation usually involves 3 steps: 1. Clean the data from null/inf values, or any biases2. Handle text and categorical data3. Rescale any unweildy featuresFirstly, we need to separate the training set into the input variables/features and the target variables/features (or labels).<jupyter_code>train_inputs = strat_train_set.drop(['male_guest_room', 'female_guest_room', 'living_room'], axis=1)
train_targets = strat_train_set[['male_guest_room', 'female_guest_room', 'living_room']].copy()<jupyter_output><empty_output><jupyter_text>3.1 Cleaning the data The data supplied has a number of quirks, as previously shown. These include the rounding of the number of days that each room has been used, and the ambiguity of the `living_room` responses (the responses were split between people who assumed that the questions meant the lounge, and the entrance hall - which was the intended meaning - hence there is a massive skew to 365 days if the lounge interpretation was made). There are a number of possible ways to deal with this:For the `living_room` we could:* Remove the entire attribute.* Split the data in two, and have a number of unfilled values for each of those catagories. However, having a `lounge` attribute will give no interesting data as all answers will be 365. This leaves a partially filled `living_room` attribute.* With a partially filled living room attribute we could remove the empty rows, set the empty values to zero or to the median value (use `from sklearn.impute import SimpleImputer`).For the rounded attributes we could:* Add some noise to each datapoint, e.g. $x_{new} = x_{old} + \mathbb{U}(-5, 5)$, a uniform random number in the range $-5 < X < 5$.* Round the data using the .cut() method in pandas, and therefore treat it as a catagorical feature.* Leave the data as-is.However, as I am not sure what to do with these now, they have just been left as is.<jupyter_code># every custom Transformer class must have a .fit(), a .transform() and a .fit_transform() method.
# Inheriting from BaseEstimator gives you the extra methods .set_params() and .get_params()
# Inheriting from TransformerMixin gives you the .fit_transform() method for free (as long as you avoid *args and **kwargs)
from sklearn.base import BaseEstimator, TransformerMixin
class DataCleaner(BaseEstimator, TransformerMixin):
    
    def __init__(self, add_extra_attributes=False):
        self.add_extra_attributes = add_extra_attributes
        
    def fit(self, X, y=None):
        return self
    
    def transform(self, X, y=None):
        '''
        here we could add code to add additional columns, or remove attributes, 
        or add noise to the data, as described above.        
        '''
        return self<jupyter_output><empty_output><jupyter_text>3.2 Categorical Locations If the location attribute was given as text we could us sklearn's `OrdinalEncoder` class, however, this assigns each catagory a number and this does not make sense for the location category as this is an **unordered category** (an ordered category would be "bad", "average", "good", "excellent", where each category is "close" to its neighbours). The standard way to deal with unordered categories is to use a one-hot encoding. We can use sklearn's `OneHotEncoder` class to easily do this.Note `OneHotEncoder`'s `fit_transform` method creates a SciPy *sparse martix*, which stores the value and location of each entry rather than the full matrix (which can be very memory intensive).<jupyter_code>from sklearn.preprocessing import OneHotEncoder
cat_encoder = OneHotEncoder(categories='auto')  # create an instance of the OneHotEncoder class
locations = [[i] for i in train_inputs['location'].values]  # need values to be a list of lists
location_one_hot = cat_encoder.fit_transform(locations)
print(location_one_hot.toarray())
cat_encoder.categories_<jupyter_output><empty_output><jupyter_text>3.3 Feature Scaling Input variables of ML algorithms almost always need to be scaled or transformed in order to make them useable. Here we are aiming to get values within (or near enough) the -1 to +1 range, and also we aim to transform each feature such that it either looks like a uniform or a normal distribution. The target features do not usually need to be changed.The two most common ways of scaling features are min-max scaling and standardisation (sometimes called normalisation). Standardisation is simply subtracting the mean of the data and dividing by the standard deviation, and is therefore much less affected by outliers than min-max scaling. Since this dataset has many outliers, we shall use standardisation. Note that it is important to only fit these transformations to the training data and not the test data.<jupyter_code>from sklearn.preprocessing import StandardScaler<jupyter_output><empty_output><jupyter_text>3.4 Pipelining A nice way of connecting all of these things together is to use sklearn's pipelining functionality. Pipelining allows us to sequence transformations in a way such that reusing and editing these transformations is simple. Additionally sklearn provides a class to interact with pandas easily called ColumnTransformer.<jupyter_code>from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer<jupyter_output><empty_output><jupyter_text>We need to define a pipeline for the numerical attributes and the categorical attributes seperately, and then combine them.<jupyter_code>num_steps = [
    # ('data_cleaner', DataCleaner()),   # some list of Estimators and Transformers
    ('std_scaler', StandardScaler()),  # all but the last one must have a .fit_transform() method
]
num_pipeline = Pipeline(num_steps)
attributes = train_inputs.columns

num_attributes = list(train_inputs.drop(['location', 'timestamp'], axis=1).columns)
cat_attributes = ['location']
print(num_attributes, cat_attributes)

full_pipeline = ColumnTransformer([
    ('num', num_pipeline, num_attributes),
    ('cat', OneHotEncoder(categories='auto'), cat_attributes),
])
train_inputs_prepared = full_pipeline.fit_transform(train_inputs)
train_inputs_prepared.shape
prepared = pd.DataFrame(train_inputs_prepared, columns=['family_member', 'age', 'male_members',
       'female_members', 'building_age', 'land_size', 'outside_room', 'total_family_members', 'l1', 'l2', 'l3', 'l4', 'l5'])
prepared.head()
prepared.describe()<jupyter_output><empty_output><jupyter_text>4. Train a Model [Back to top](top) There are loads of different types of model provided by sklearn, below a number of them are tried and evaluated.<jupyter_code>from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_val_score, GridSearchCV, RandomizedSearchCV<jupyter_output><empty_output><jupyter_text>4.1 Basic Models<jupyter_code>lin_reg = LinearRegression()
lin_reg.fit(train_inputs_prepared, train_targets)<jupyter_output><empty_output><jupyter_text>Fitting a model can be as simple as that. Now all that is needed is to see how it did. We can do this by just doing some predictions on a subset of the training data, or by looking at a performance metric, such as RMSE (root mean squared error).<jupyter_code>some_data = train_inputs.iloc[:5]
some_targets = train_targets.iloc[:5]
some_data_prepared = full_pipeline.transform(some_data)
print('predictions: ', lin_reg.predict(some_data_prepared))
print('Labels: ', some_targets)
lin_predictions = lin_reg.predict(train_inputs_prepared)
lin_mse = mean_squared_error(train_targets, lin_predictions)
lin_rmse = np.sqrt(lin_mse)
print(lin_rmse)<jupyter_output>82.2232976197475<jupyter_text>Given that the target variables have the following statistics, this isn't very good.| data  | `male_guest_room` 	| `female_guest_room` 	| `living_room` 	||------	|-------------------	|:-------------------:	|---------------	|| Mean 	|       39.86       	|        33.70        	|     323.62    	|| Std  	|       79.10       	|        70.70        	|     96.48     	|We can try a different type of model to solve this.<jupyter_code>tree_reg = DecisionTreeRegressor()
tree_reg.fit(train_inputs_prepared, train_targets)
tree_predictions = tree_reg.predict(train_inputs_prepared)
tree_mse = mean_squared_error(train_targets, tree_predictions)
tree_rmse = np.sqrt(tree_mse)
print(tree_rmse)<jupyter_output>0.3465210660024545<jupyter_text>This is tiny! Clearly this is overfitting.  4.2 Cross-Validation and Ensemble Methods Cross validation splits the dataset into k *folds*, or subsets. The model is then trained k times, picking a different fold as the test set each time. This allows for better evaluation of the model in general. Note that sklearn's cross-validation requires a utility function rather than a cost function so we just negate the mean squared error.<jupyter_code>tree_scores = cross_val_score(tree_reg, train_inputs_prepared, train_targets, scoring='neg_mean_squared_error', cv=10)
tree_rmses = np.sqrt(-tree_scores)
lin_scores = cross_val_score(lin_reg, train_inputs_prepared, train_targets, scoring='neg_mean_squared_error', cv=10)
lin_rmses = np.sqrt(-lin_scores)

print("Tree: ", tree_rmses, tree_rmses.mean(), tree_rmses.std())
print("Lin: ", lin_rmses, lin_rmses.mean(), lin_rmses.std())<jupyter_output>Tree:  [105.76631924 120.89336112 124.06803049 119.82771362 111.24008939
 118.08328347 112.23885417 106.10886762 116.18629884 140.72146745] 117.5134285423917 9.673035975309215
Lin:  [101.40300385  82.53788236  96.06284975  92.64017876  84.55963843
  78.74102335  76.78144041  77.41080709  74.83302167  92.12868106] 85.70985267236362 8.772316945803547<jupyter_text>The tree regressor is now even worse than the linear regressor! Note that cross-validation also allows us to measure how precise the estimate of the error is, with the standard deviation.Often we can get very good performance using one of these models that are poor on their own by using an **Ensemble Method**. Random forests are an ensemble method for decision trees that work by training many decision trees on random subsets of data and averageing out their predictions.<jupyter_code>forest_reg = RandomForestRegressor()
forest_reg.fit(train_inputs_prepared, train_targets)

forest_scores = cross_val_score(forest_reg, train_inputs_prepared, train_targets, scoring='neg_mean_squared_error', cv=10)
forest_rmses = np.sqrt(-forest_scores)
print("Forest: ", forest_rmses, forest_rmses.mean(), forest_rmses.std())<jupyter_output>C:\Users\Alex\Miniconda3\lib\site-packages\sklearn\ensemble\forest.py:248: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)<jupyter_text>This actually performed even worse than the regular decision tree regressor. However, the result is more general, which suggests that these models are not complex enough for the data at hand - or that the input data does not capture any relationship with the target data.  4.3 Hyperparameter TuningHyperparameters are paramters govern how the model is set up, e.g. the use of 10 folds in the cross_val above, or the squared exponent in the linear regressor's discriminator - they can be anything that you set before learning begins. Tuning them is a good way of improving your model.<jupyter_code>param_grid = [
    {'n_estimators': [3, 10, 30, 100], 'max_features': [2, 4, 6, 8]},  # try to go up in powers of 10 (note log(3) ~ 0.5 so these increase as 0.5, 1, 1.5, ... in the log scale)
    {'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [2, 3, 4]},    
]
forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5,
                          scoring='neg_mean_squared_error',
                          return_train_score=True)
grid_search.fit(train_inputs_prepared, train_targets)
print(grid_search.best_params_)  # sklearn stores its learnt parameters in attributes with an underscore suffix
print(grid_search.best_estimator_)
res = grid_search.cv_results_
for mean_score, params in zip(res['mean_test_score'], res['params']):
    print(np.sqrt(-mean_score), params)<jupyter_output>{'max_features': 2, 'n_estimators': 100}
RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
           max_features=2, max_leaf_nodes=None, min_impurity_decrease=0.0,
           min_impurity_split=None, min_samples_leaf=1,
           min_samples_split=2, min_weight_fraction_leaf=0.0,
           n_estimators=100, n_jobs=None, oob_score=False,
           random_state=None, verbose=0, warm_start=False)
96.93361790079484 {'max_features': 2, 'n_estimators': 3}
88.32521682120553 {'max_features': 2, 'n_estimators': 10}
85.60279610212307 {'max_features': 2, 'n_estimators': 30}
84.57066586505252 {'max_features': 2, 'n_estimators': 100}
98.32590945453207 {'max_features': 4, 'n_estimators': 3}
88.54517173187007 {'max_features': 4, 'n_estimators': 10}
86.4852381497045 {'max_features': 4, 'n_estimators': 30}
84.9688168943982 {'max_features': 4, 'n_estimators': 100}
98.44029389026545 {'max_features': 6, 'n_estimators': 3}
89.48603985604291 {'max_features': 6, 'n_estimators': 10}
[...]<jupyter_text>Since 300 is the max n_estimators tries, it is probably a good idea to gridsearch through higher values. Still not much improvement, but this was expected as I don't think the data is sufficient. However, this does show you how to conduct a grid search. Another (better if lots of hyperparams) way of doing this is through a randomised search `RandomizedSearchCV`.  4.4 Analyse the Model<jupyter_code>feature_importances = grid_search.best_estimator_.feature_importances_
#print(feature_importances)
cat_encoder = full_pipeline.named_transformers_['cat']
cat_one_hot_attributes = list(cat_encoder.categories_[0])
attributes = num_attributes + cat_one_hot_attributes

sorted(zip(feature_importances, attributes), reverse=True)<jupyter_output><empty_output><jupyter_text>From this we can see that `building_age`, `age` and `land_size` are the most important catagories. The location category is of almost no importance. However, take this with a grain of salt, as the model has not performed better than by just taking the mean of the target values.  5. Evaluate [Back to top](top)Finally we evaluate the model on the test set. To do this we call `.transform()`, not `.fit_transform()`.<jupyter_code>final_model = grid_search.best_estimator_
test_inputs = strat_test_set.drop(['male_guest_room', 'female_guest_room', 'living_room'], axis=1)
test_targets = strat_test_set[['male_guest_room', 'female_guest_room', 'living_room']].copy()

test_prepared = full_pipeline.transform(test_inputs)
final_predictions = final_model.predict(test_prepared)
final_mse = mean_squared_error(test_targets, final_predictions)
final_rmse = np.sqrt(final_mse)
print(final_rmse)<jupyter_output>81.76040767873596
{'fraction_non_alphanumeric': 0.08824849993877301, 'fraction_numerical': 0.03249112208661578, 'mean_word_length': 5.5718884120171674, 'pattern_counts': {'":': 0, '<': 101, '<?xml version=': 0, '>': 100, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '70883', 'license': 'MIT', 'n_tokens_mistral': 7911, 'n_tokens_neox': 7137, 'n_words': 2678, 'path': 'SaudiArchitecture/SaudiArchitecure.ipynb', 'repo_name': 'alexdarch/PythonProjects', 'size': '24499'}