{"cells":[{"metadata":{"_uuid":"233ebcdcfdd40d5fcf16cf0c10a440f791a49511"},"cell_type":"markdown","source":"# Housing Prices prediction with real life examples\n\nThis is a solution for the Housing Prices challenge on Kaggle with real life examples.\nhttps://www.kaggle.com/c/home-data-for-ml-course"},{"metadata":{"trusted":true,"_uuid":"8bf8a648dee8922c6958f88e03c4412775ef0633"},"cell_type":"code","source":"from IPython.display import HTML\nHTML('')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b1cc1de7925838ebcc17ca5aadc047a952f0d88f"},"cell_type":"markdown","source":"### Imported Python packages"},{"metadata":{"trusted":false,"_uuid":"06a1e61264a6505326c0df13e3e79a7c953cb6e8"},"cell_type":"code","source":"# For simple vectorized calculations\nimport numpy as np\n\n# Mainly data handling and representation\nimport pandas as pd\n\n# For statisctics\nimport scipy\nfrom scipy import stats\n\n# Models\nfrom sklearn.ensemble import GradientBoostingRegressor\nfrom xgboost import XGBRegressor\n\n# Data preparation\nfrom sklearn.preprocessing import MinMaxScaler, PowerTransformer, OneHotEncoder\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.model_selection import train_test_split\n\n# Model validation, scoreing\nfrom sklearn.metrics import mean_absolute_error, make_scorer\nfrom sklearn.model_selection import cross_val_score\n\n# Math helper\nfrom math import isnan\n\n# Plotting and display\nfrom IPython.display import display\nfrom matplotlib import pyplot as plt\nimport seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":false,"_uuid":"27687e25af685a12cf554950fd4eef3ae6e05adb"},"cell_type":"code","source":"# Path of the file to read.\ntrain_file_path = '../input/train.csv'\n\n# Read the file\nhome_data = pd.read_csv(train_file_path)\n\n# The shape of the data\nhome_data.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"829fedf5b9652cc7a3f0979b2fc7f6633ac10722"},"cell_type":"markdown","source":"There are 1460 measurements of the dataset and 81 features were measures. One of the features is the house price which we would like to predict. So it leaves 80 features to estimate the prices."},{"metadata":{"_uuid":"08c2f7ea59d9fe408a93f848ee91e4738628c6fe"},"cell_type":"markdown","source":"## Data inspection and preprocessing\n### Inspecting and preprocessing y"},{"metadata":{"trusted":false,"_uuid":"690278f9cadb90dfcba6246c5ae21f0716bc5797"},"cell_type":"code","source":"# Create target object and call it y_orig, because we will format this value later\ny_orig = pd.DataFrame(home_data.SalePrice)\n# Let's see the distribution for y\nsns.distplot(y_orig)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1438d6c4922a0191b3ce74a41d2943a7c1e1dd2a"},"cell_type":"markdown","source":"The distribution is skewed because the number of houses sold on a given price shows logarithmic like behaviur. There are more expensive houses than cheap houses.\n\nWe should correct the skewness because the learning algorithm estimates easier/better and it is easier to find and cut the outliers.\n\nWe could simply use a logarithm here but the there are multiple options. For example the PowerTransformer from the Scikit Learn preprocessing toolkit. Here is a great article explaining it: https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html\n\nNote that we need to inverse transform the y later to get back to the original scale of the housing prices.\n\n\"PowerTransformer applies a power transformation to each feature to make the data more Gaussian-like. ... The power transform finds the optimal scaling factor to stabilize variance and mimimize skewness through maximum likelihood estimation. ...\"\n\nWe can compare the log of the y and the powertransformer's performance by measuring the skew and kurtosis."},{"metadata":{"trusted":false,"_uuid":"5b6008a0ee9f107e7f8dd6af122ce7fd42e180e5"},"cell_type":"code","source":"# Calculate log(y)\ny_log = np.log(y_orig)\n\n# Calculate y_power by PowerTransformer\ny_power_scaler = PowerTransformer()\ny_power_orig = y_power_scaler.fit_transform(y_orig)\n\n# Calculating skews\norig_skew = scipy.stats.skew(y_orig)[0]\nlog_skew = scipy.stats.skew(y_log)[0]\npower_skew = scipy.stats.skew(y_power_orig)[0]\n\n# Calculating kurtosises\noriginal_kurtosis = scipy.stats.kurtosis(y_orig)[0]\nlog_kurtosis = scipy.stats.kurtosis(y_log)[0]\npower_kurtosis = scipy.stats.kurtosis(y_power_orig)[0]\n\ncolumns = [\"skew\", \"kurtosis\"]\n\nkurtosis_skew_comparison = pd.DataFrame([[orig_skew, original_kurtosis],\n [log_skew, log_kurtosis],\n [power_skew, power_kurtosis],],\n index=[\"original\", \"logarithm\", \"power\"],\n columns=columns)\nkurtosis_skew_comparison","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"99198a9fb2f13311d25af5b0595535cf93b067d3"},"cell_type":"code","source":"# Plotting distribution of log_y\nfig = plt.figure(1)\nfig.set_size_inches(18.5, 6)\nplt.subplot(121)\nplt.title('Log(y)')\nplt.xlabel('log(y)')\nplt.ylabel('distribution')\n\nsns.distplot(y_log)\n\n# Plotting distribution of power_y\nplt.subplot(122)\nplt.title('power_y')\nplt.xlabel('power_y')\nplt.ylabel('distribution')\nsns.distplot(y_power_orig)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"228194b508e3ab8094f830a0bfb63459ff1e37c1"},"cell_type":"code","source":"# Dropping y outliers - the model training showed that there is no need to drop any outliers\nfrom scipy.stats import norm\n\nplt.plot(np.sort(norm.pdf(y_power_orig), axis=0))\nplt.title('Normal distribution probabilities ordered.')\nplt.xlabel('Datapoints')\nplt.ylabel('Probability')\n\nplt.show()\nthreshold = 0.00\ny_filter = [False if x < threshold else True for x in norm.pdf(y_power_orig)]\n\ny_power_pd = pd.DataFrame(y_power_orig)\nprint(\"New shape of the y: {}\".format(y_power_pd.shape))\n\ny_power = y_power_pd.iloc[y_filter, 0].values.reshape(-1,1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"366c67188f124aba2bb450a27d0925a6d65f3e8d"},"cell_type":"markdown","source":"We can colnclude that the best choice is power_y because of the lowest skewness."},{"metadata":{"_uuid":"feccf02994cfd4540cfaeba6fa47ed89c1395457"},"cell_type":"markdown","source":"## Inspecting and transforming X\n\nFirst the Padas can show us a good summary of the dataset:"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"31ac6ed31518ddeb3e376630e5ad661dca70b039"},"cell_type":"code","source":"home_data.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f752083e7d29ae0279677b3d8cc32cfad9e3b1f4"},"cell_type":"markdown","source":"We can see that there are 38 columns but shouldn't be 80? If we inspect the datatypes of the data we can see that there are object type columns also. This means those are categorical values represented by strings. Let's see:"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"9f05f3ac3510204b19011307c19e7c37d2b25bec"},"cell_type":"code","source":"home_data.info()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"79093cffa9a18d0438920d22d6dbf5dbe5a51a91"},"cell_type":"markdown","source":"At the bottom you can see which and how many types there are. float64(3), int64(35), object(43)\nAnd you can see that the number of measurement by feature is not always 1460 so there are missing values. We address this problem later.\nAlso we definitely doesn't need the Id column because the house price does not depend on it."},{"metadata":{"trusted":false,"_uuid":"fd8a2312a645d626ad2c8a3a2c79273d98b0d07d"},"cell_type":"code","source":"# Let's copy the original data.\nhome_data_copy = home_data.copy()\nif \"Id\" in home_data_copy.columns:\n home_data_copy = home_data_copy.drop([\"SalePrice\", \"Id\"], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7cc63cd7d18d661b1ad50202a626c663f63a911d"},"cell_type":"markdown","source":"The categorical data can be separated by the other features. and the learning algorithm can use it using one-hot encoding. Let's define a function which separates the categorical features from the other features."},{"metadata":{"_uuid":"44f3d9d87d9900954a76a02c3f7ea8a3eb6d6011","trusted":false},"cell_type":"code","source":"# Separate the categorical values from the numerical values\ndef separate_X_numerical_categorical(X):\n \n X_numericals = X.copy()\n X_categoricals = X.copy()\n \n # Loop over the columns\n for column in X.columns[1:]:\n if str(X[column].dtype) == \"object\":\n X_numericals.drop(column, axis=1, inplace=True)\n else:\n X_categoricals.drop(column, axis=1, inplace=True)\n \n return X_numericals, X_categoricals","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"2b7e14470450ddf9d4fc19f56e59918a61b0c77b"},"cell_type":"code","source":"# Fill nan values with \"nan\" strings\ndef fill_nan_X_categoricals(X_categoricals):\n \n X_categoricals_filled = X_categoricals.fillna(value=\"nan\")\n \n return X_categoricals_filled","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"9f6b2224a0304841eccedb6ee04720c51105275f"},"cell_type":"code","source":"# One-hot encode the categorical values\ndef one_hot_encode_categories(X_categoricals):\n \n encoder = OneHotEncoder(handle_unknown='ignore', sparse=False)\n\n X_cat_one_hot = pd.DataFrame(encoder.fit_transform(X_categoricals), columns=encoder.get_feature_names())\n \n return X_cat_one_hot","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"29226548c4eb5bab9f4aaf8f06e054b37b7f8211"},"cell_type":"code","source":"def fill_nan_X_numericals(X_numericals, my_imputer=None):\n \n if my_imputer is None:\n # Imputation\n my_imputer = SimpleImputer()\n \n X_numericals_filled = my_imputer.fit_transform(X_numericals)\n else:\n X_numericals_filled = my_imputer.transform(X_numericals)\n \n return pd.DataFrame(X_numericals_filled, columns=X_numericals.columns), my_imputer","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"ca8f731d60740bfc07a2fd7f34f73753daa6e926"},"cell_type":"code","source":"def X_numericals_transform(X_numericals, X_scaler=None):\n \n if X_scaler == None:\n # Create transformer\n X_scaler = PowerTransformer()\n \n # Fit and transform\n X_numericals_scaled = X_scaler.fit_transform(X_numericals)\n \n else:\n # Only transform\n X_numericals_scaled = X_scaler.transform(X_numericals)\n \n return pd.DataFrame(X_numericals_scaled, columns=X_numericals.columns), X_scaler","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"1476809890c2c38aa41ee87aa55f9f1c8d98e587"},"cell_type":"code","source":"from sklearn.decomposition import PCA\n\ndef PCA_transform(X_numericals, PCA_transformer=None):\n if PCA_transformer == None:\n # Create transformer\n PCA_transformer = PCA()\n \n # Fit and transform\n X_numericals_scaled = PCA_transformer.fit_transform(X_numericals)\n \n else:\n # Only transform\n X_numericals_scaled = PCA_transformer.transform(X_numericals)\n \n return pd.DataFrame(X_numericals_scaled, columns=X_numericals.columns), PCA_transformer","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2ca55624fdf696236057c718f6dcef8e68ec501c"},"cell_type":"markdown","source":"And add the previously seen but not present categories to the test features."},{"metadata":{"trusted":false,"_uuid":"1e2c1a924c362bf2b6dda00198919b49a6ce9eb7"},"cell_type":"code","source":"def drop_unknown_categories(X_categoricals_test, X_categoricals):\n \n for column in X_categoricals_test.columns:\n if column not in X_categoricals.columns:\n X_categoricals_test.drop(column, axis=1, inplace=True)\n \n return X_categoricals_test","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8cfbf6713536dd7af2685f2115bce3e40b89f83f"},"cell_type":"markdown","source":"To make the new categories one-hot vectos in sync with the new data drop the previously unseen categories."},{"metadata":{"trusted":false,"_uuid":"166118c4da39c4c1d1d638fb7fe352dd4fb64bde"},"cell_type":"code","source":"def add_known_categories(X_categoricals_test, X_categoricals):\n \n for column in X_categoricals.columns:\n if column not in X_categoricals_test.columns:\n X_categoricals_test = pd.concat([X_categoricals_test,\n pd.DataFrame(np.zeros((X_categoricals_test.shape[0])).reshape(-1,1), columns=[column])],\n axis=1).reset_index(drop=True)\n\n return X_categoricals_test","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"225260b00ed39c7c1beb3196425e197e8f9de549","trusted":false},"cell_type":"code","source":"# Separate the categoricals from the numericals\nX_numericals, X_categoricals = separate_X_numerical_categorical(home_data_copy.iloc[y_filter])\n\n##############\n# Fill X_categoricals nan values with \"nan\"\nX_cat_filled = fill_nan_X_categoricals(X_categoricals)\n\n# On hot encode the categoricals\nX_cat_one_hot = one_hot_encode_categories(X_cat_filled)\n\n\n##############\n# Fill X_numericals nans with imputing\nX_numericals_filled, my_imputer = fill_nan_X_numericals(X_numericals)\n\n# Transform also the X_numericals values by power transformer\nX_numericals_scaled, X_scaler = X_numericals_transform(X_numericals_filled)\n\n# The principal component analysis did not help the model\n#X_numericals_PCA, X_PCA_transformer = PCA_transform(X_numericals_scaled)\n\n#X_numericals_PCA.describe()\nX_numericals_outlier = X_numericals_scaled","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"224f7c70959c7dcfb4e6f821100a69250de01c8f"},"cell_type":"markdown","source":"### Calculate the multivariate normal distribution for outlier detection"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"060bcb866db7e6236913809bb66def3c4c1c90ac"},"cell_type":"code","source":"# Calculate the mean of X_numericals_scaled\nX_mean = np.mean(X_numericals_outlier.values, axis=0, keepdims=True)[0]\n\n# Calculate the covariance matrix of X_numericals_scaled\nX_covariance_matrix = np.cov(X_numericals_outlier.values.T)\n\n# Calculate the multivariate normal distribution, which is a matrix filled with measurement probabilities being not outlier.\nMGD = stats.multivariate_normal.pdf(X_numericals_outlier, mean=X_mean, cov=X_covariance_matrix)\n\n# Use the log for plotting so the outliers can be visualised better\nplt.plot(np.log(np.sort(MGD)))\nplt.title('Multivariate normal distribution probabilities ordered.')\nplt.xlabel('Datapoints')\nplt.ylabel('Probability')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2a7b44609a7d81fbd4d24a15dff4a7b04586f6dd"},"cell_type":"markdown","source":"### Cut the outliers from X by Multivariate normal distribution"},{"metadata":{"_uuid":"765a01153323b05b28d81386473abafe7fb39ea4","scrolled":true,"trusted":false},"cell_type":"code","source":"# Select a threshold\nthreshold = -150\n\nfilter_array = [True if np.log(meas) > threshold else False for meas in MGD]\n\nX_numericals_outlier_scaled, X_scaler2 = X_numericals_transform(X_numericals_outlier)\n\nprint(\"Original number of measurement points: {}\".format(X_cat_one_hot.shape[0]))\n\n# Remove the outliers from the numerical features\nX_numericals_outlierless = X_numericals_outlier_scaled.iloc[filter_array]\n\n# Remove the outliers from the categorical features\nX_categoricals_outlierless = X_cat_one_hot.iloc[filter_array]\n\n# Concatenate the numerical and categorical features\nX_outlierless = pd.concat([X_numericals_outlierless, X_categoricals_outlierless], axis=1)\n\n# Remove the outliers from the output\nY_outlierless = y_power[filter_array]\n\nprint(\"New number of measurement points: {}\".format(Y_outlierless.shape[0]))\nplt.plot(np.log(np.sort(MGD)))\nplt.plot(np.ones(len(MGD))*threshold)\nplt.title('Multivariate normal distribution probabilities ordered.')\nplt.xlabel('Datapoints')\nplt.ylabel('Probability')\nplt.show()\n\n# The outliers are cut below the orange line","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"7530850776be2ac6a542e13ad7c6e9841fd0d8d6"},"cell_type":"code","source":"# Concatenate the X_categoricals and X_numericals\nX = X_outlierless\ny = Y_outlierless","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"71218c92faa9e4e5f7d42e5c43ae5ccc2834cde1","scrolled":true,"trusted":false},"cell_type":"code","source":"print(\"Final shape of the features: {}\".format(X.shape))\nprint(\"Final shape of the features: {}\".format(y.shape))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"10fdb043f6f650ac543ada4b14054057d8725a98","scrolled":true,"trusted":false},"cell_type":"code","source":"# Define the score metrics to be the mean absolute error\ndef mae(predict, actual):\n score = mean_absolute_error(predict, actual)\n return score\n\nmae_score = make_scorer(mae)\n\n# Define the cross validation function\ndef score(model, X, y):\n score = cross_val_score(model, X, (y.reshape(-1)), cv=5, scoring=mae_score).mean()\n return score","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d9486fe36977a45657087c185c41fcd22c71fa77"},"cell_type":"markdown","source":"### Grid search using cross validation"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"15aac50e13ea436d10a281a96c45236ef083c983"},"cell_type":"code","source":"hyperparameters = pd.DataFrame()\n\nmodels = []\nscores = []\n\n# grid search\nfor i in range(10):\n \n # The hyperparameters:\n \n # The learning rate between 0.03 and 0.01\n learning_rate = 0.03 - 0.02 * (np.random.rand())\n \n # The max depth between 5-6\n max_depth = np.random.randint(3, 6)\n \n # The number of estimators between 1000 and 3000\n n_estimators = np.random.randint(1000, 3000)\n \n # Base score between 0.8 and 0.4\n base_score = 0.8 - 0.4 * (np.random.rand())\n \n # Subsample between 0.8 and .04\n subsample = 0.8 - 0.4 * (np.random.rand())\n \n # Create the model\n model = XGBRegressor(max_depth=max_depth,\n learning_rate=learning_rate,\n n_estimators=n_estimators,\n silent=True,\n objective='reg:linear',\n booster='gbtree',\n subsample=subsample,\n base_score=base_score,\n random_state=0,\n importance_type='gain')\n \n # Calculate the score\n mean_cross_validation = score(model, X, y)\n \n # Save the scores for later evaluation\n scores.append(mean_cross_validation)\n \n # Save the model for later evaluation\n models.append(model)\n \n # Append the hyperparameters for lates examination\n hyperparameters = hyperparameters.append(pd.DataFrame([[mean_cross_validation,\n learning_rate, \n max_depth,\n n_estimators,\n base_score,\n subsample]],\n columns=[\"mean_cross_validation\",\n \"learning_rate\",\n \"max_depth\",\n \"n_estimators\",\n \"base_score\",\n \"subsample\"]))\n display(hyperparameters)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c7d57b4cace092d0ca2770c93cc8ad6c2dfe9ae7"},"cell_type":"markdown","source":"### Let's see the importance of the hyperparameters by calculating the crossvalidation"},{"metadata":{"scrolled":false,"trusted":false,"_uuid":"1aed21bcedca4d853e7b1f2456e99d0be410e3d7"},"cell_type":"code","source":"corr = hyperparameters.corr()\ncorr.columns = hyperparameters.columns\n\ncorr[\"index\"] = hyperparameters.columns\n\ncorr.set_index(\"index\", inplace=True)\nprint(corr)\n\nsns.heatmap(data=corr, center=0)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fc5835f1740193d288598ff95c2185db04319332"},"cell_type":"markdown","source":"### Compare results using train-validation data splitting\n\nYou can see the difference in the learning on the resulting graphs. The early stopping causes not constant learning steps. The gap between the "},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"ed389ed3cbd780ba391100788ec27ce97710d8e3"},"cell_type":"code","source":"# Split into validation and training data\ntrain_X, val_X, train_y, val_y = train_test_split(X, y, random_state=1)\n\nfor model in models:\n eval_set = [(train_X, train_y), (val_X, val_y)]\n eval_metric = [\"mae\"]\n \n # Fit the model\n %time history = model.fit(train_X, train_y, eval_metric=eval_metric, eval_set=eval_set, verbose=False, early_stopping_rounds=150)\n\n # Make the predictions on the training set\n train_predictions = model.predict(train_X)\n\n # Rescale the output to original\n train_predictions_scaled = y_power_scaler.inverse_transform(train_predictions.reshape(-1,1))\n train_y_scaled = y_power_scaler.inverse_transform(train_y)\n \n # Calculate the Mean Average Error of the training set\n train_mae = mean_absolute_error(train_predictions_scaled.astype(\"float64\"), train_y_scaled.astype(\"float64\"))\n print(\"Mean average error of the training set: {}\".format(train_mae))\n\n # Make the predictions on the validation set\n val_predictions = model.predict(val_X)\n \n # Rescale the output to original\n val_predictions_scaled = y_power_scaler.inverse_transform(val_predictions.reshape(-1,1))\n val_y_scaled = y_power_scaler.inverse_transform(val_y)\n \n # Calculate the Mean Average Error of the validation set\n val_mae = mean_absolute_error(val_predictions_scaled.astype(\"float64\"), val_y_scaled.astype(\"float64\"))\n print(\"Mean average error of the validation set: {}\".format(val_mae))\n\n fig = plt.figure()\n ax = plt.subplot(111)\n ax.plot(history.evals_result()[\"validation_0\"][\"mae\"], label='training set')\n ax.plot(history.evals_result()[\"validation_1\"][\"mae\"], label='validation set')\n ax.legend()\n plt.title('Score during training.')\n plt.xlabel('Training step')\n plt.ylabel('Score (MAE)')\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a565a792f7a4d0786dfb04216e806f0f684ced38"},"cell_type":"markdown","source":"### Let's see the scores"},{"metadata":{"trusted":false,"_uuid":"3fcd544c6e715f2367a832fc787984005243fbab"},"cell_type":"code","source":"scores","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1ff5369d6b9c2a3459875263283ebd61861ed211"},"cell_type":"markdown","source":"### Chose the best model"},{"metadata":{"trusted":false,"_uuid":"8fd2bc0602cd5112fb577e216b033abd31a55b43"},"cell_type":"code","source":"model = models[np.argmin(np.array(scores))]\nprint(model)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c1aeba7c004d9005ae489f0e6068890bc44ea708"},"cell_type":"markdown","source":"### Fit the chosen model using all of the data"},{"metadata":{"trusted":false,"_uuid":"b42362f3780021d9b4ab4dec3b260888d0c1146c"},"cell_type":"code","source":"model.fit(X, y.reshape(-1,1))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0d4af3e68d6291721b2afd5ea8669323e49ec648"},"cell_type":"markdown","source":"## Process the test data"},{"metadata":{"trusted":false,"_uuid":"2cea64f7961733d1c82a21b1f783347c387eda03"},"cell_type":"code","source":"def convert_test_data(test_data, X_cat_one_hot, my_imputer, X_scaler, X_scaler2):\n\n ##############\n # Separate the test set categoricals from the numericals\n X_numericals_test, X_categoricals_test = separate_X_numerical_categorical(test_data)\n\n ##############\n # Fill X_categoricals_test nans with \"nan\"\n X_categorical_filled_test = fill_nan_X_categoricals(X_categoricals_test)\n\n # One-hot encode X_categoricals_test\n X_cat_one_hot_test = one_hot_encode_categories(X_categorical_filled_test)\n\n # Add missing categories\n X_cat_one_hot_test_dropped = add_known_categories(X_cat_one_hot_test, X_cat_one_hot)\n\n # Remove unkown categories\n X_cat_one_hot_filled_test = drop_unknown_categories(X_cat_one_hot_test_dropped, X_cat_one_hot)\n\n ##############\n # Fill X_numericals_test nans with imputing\n X_numericals_filled_test, _ = fill_nan_X_numericals(X_numericals_test, my_imputer)\n\n # Transform also the X_numericals_test values by power transformer\n X_numericals_test_scaled, _ = X_numericals_transform(X_numericals_filled_test, X_scaler)\n\n X_numericals_test_scaled2, _ = X_numericals_transform(X_numericals_test_scaled, X_scaler2)\n\n #X_numericals_test_PCA, _ = PCA_transform(X_numericals_test_scaled, X_PCA_transformer)\n\n X_test = pd.concat([X_numericals_test_scaled2, X_cat_one_hot_filled_test], axis=1)\n print(\"Test data dimensions: {}\".format(X_test.shape))\n \n return X_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"45ec36c46a1943a6ac4b7b4c714c81d69a5d0db9"},"cell_type":"code","source":"# path to file you will use for predictions\ntest_data_path = '../input/test.csv'\n\n# read test data file using pandas\ntest_data_orig = pd.read_csv(test_data_path)\n\n# Drop the index column\nif \"Id\" in test_data_orig.columns:\n test_data = test_data_orig.drop([\"Id\"], axis=1)\n\n# Preprocess the test data\nX_test = convert_test_data(test_data, X_cat_one_hot, my_imputer, X_scaler, X_scaler2)\n\n# Sort the colums like the original in case of misalignements\nX_test = X_test[X.columns]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8bb348e6e4ec800975e9cd567cfe97ec0e49eae6"},"cell_type":"markdown","source":"## Save the output"},{"metadata":{"_uuid":"8c9dc12671da530ad1139ebbc7f7ffef01586073","trusted":false},"cell_type":"code","source":"# Make predictions which we will submit.\ntest_preds_unscaled = model.predict(X_test).reshape(-1, 1)\n\n# Inverse transform the predictions to the original scale\ntest_preds = y_power_scaler.inverse_transform(test_preds_unscaled)[:,0]\n\n# Save the predictions\noutput = pd.DataFrame({'Id': test_data_orig['Id'],\n 'SalePrice': test_preds})\n\noutput.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"28407a73956cae574cca1af863040071e7401e63"},"cell_type":"markdown","source":"Remark: An improvement can be made if we would first include the test data features in our features before normalization because the distributions would fit better the real life and outliers could be found more easily. This means in practice that we need to retrain our network for every new measurement or at least time-to-time. There are application where the additional time cost does not allow this."},{"metadata":{"_uuid":"819f892188daccb69fd3bd7962c3ae1708b66f4b"},"cell_type":"markdown","source":"# Real life example - Houses in Ames on Craigslist\n\n## Webscrape real life data"},{"metadata":{"trusted":false,"_uuid":"a929dde0594b21c879d6f80000151dcaee96a914"},"cell_type":"code","source":"from requests import get\nfrom requests.exceptions import RequestException\nfrom contextlib import closing\nfrom bs4 import BeautifulSoup\nimport pandas as pd\n\ndef get_data(url):\n raw_html = simple_get(url)\n html = BeautifulSoup(raw_html, 'html.parser')\n return parse_html(html)\n\ndef parse_html(html):\n SalePrice = 0 # Done\n YearSold = 2019 # Done\n FullBathroom = 0 # Done\n HalfBathroom = 0 # Done\n Bedroom = 0 # Done\n Utilities = \"AllPub\" # Done\n Garage = 0\n Fence = 0\n\n SalePrice = int(html.find_all('span', class_='price')[0].text[1:])\n brbaline = html.find_all('span', class_='shared-line-bubble')[0].text\n Bedroom = int(brbaline[:brbaline.index(\"BR\")])\n bath = float(brbaline[brbaline.index(\"/\")+1:brbaline.index(\"Ba\")])\n FullBathroom = int(bath)\n HalfBathroom = int((bath-FullBathroom)>0)\n\n Garage = len([1 for x in html.select('p.attrgroup span') if \"garage\" in x.text.lower()])>0\n\n for i, sec in enumerate(html.select('section')):\n if sec.get(\"id\") is not None:\n if \"postingbody\" in sec.get(\"id\"):\n if \"fence\" in sec.text.lower():\n Fence = 1\n if \"garage\" in sec.text.lower():\n Garage = 1\n\n output = pd.DataFrame([[SalePrice,YearSold,FullBathroom,HalfBathroom,Bedroom,Utilities,Garage,Fence]],\n columns=[\"SalePrice\", \"YrSold\", \"FullBath\", \"HalfBath\", \"BedroomAbvGr\", \"Utilities\", \"GarageCars\", \"Fence\"])\n return output\n\ndef simple_get(url):\n \"\"\"\n Attempts to get the content at `url` by making an HTTP GET request.\n If the content-type of response is some kind of HTML/XML, return the\n text content, otherwise return None.\n \"\"\"\n try:\n with closing(get(url, stream=True)) as resp:\n if is_good_response(resp):\n return resp.content\n else:\n return None\n\n except RequestException as e:\n log_error('Error during requests to {0} : {1}'.format(url, str(e)))\n return None\n\n\ndef is_good_response(resp):\n \"\"\"\n Returns True if the response seems to be HTML, False otherwise.\n \"\"\"\n content_type = resp.headers['Content-Type'].lower()\n return (resp.status_code == 200\n and content_type is not None\n and content_type.find('html') > -1)\n\n\ndef log_error(e):\n \"\"\"\n It is always a good idea to log errors.\n This function just prints them, but you can\n make it do anything.\n \"\"\"\n print(e)\n\n\nweburl = 'https://ames.craigslist.org/reb/d/ames-stop-renting-when-you-can-own/6796373030.html'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"354117cf71f949cda23489abf83394afbeb96269"},"cell_type":"markdown","source":"(The links are maybe outdated.)"},{"metadata":{"trusted":false,"_uuid":"e4713c219581c5a0170dd364f8a46521a0c58e4e"},"cell_type":"code","source":"urls = [\"https://ames.craigslist.org/reb/d/ames-stop-renting-when-you-can-own/6796373030.html\",\n \"https://ames.craigslist.org/reo/d/kelloggjasperth-ave-kellogg-ia/6796762804.html\",\n \"https://ames.craigslist.org/reo/d/story-city-commercial-light-industrial/6799601680.html\",\n \"https://ames.craigslist.org/reb/d/ankeny-split-foyer/6800347988.html\",\n \"https://ames.craigslist.org/reo/d/polk-city-brand-new-villa-in-wolf-creek/6777063925.html\",\n \"https://ames.craigslist.org/reo/d/omaha-walk-out-ranch/6797562631.html\",\n \"https://ames.craigslist.org/reb/d/ames-open-sat-1-12-at-130-to-3-pm/6791622440.html\",\n \"https://ames.craigslist.org/reb/d/waukee-modern-living/6784701226.html\",\n \"https://ames.craigslist.org/reo/d/ames-ames-ia-why-pay-rent-buy-and-save/6767869407.html\",\n \"https://ames.craigslist.org/reo/d/collins-house-for-sale/6798248546.html\"]\n\nreal_data_list = []\n\nfor url in urls:\n print(url)\n real_data_list.append(get_data(url))\n\nreal_data = get_data(weburl)\nreal_data = pd.concat(real_data_list)\nreal_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"dbfe4b33d968f24dadaded79ad1a5c3c0d4f5ec0"},"cell_type":"code","source":"scraped_columns = [\"SalePrice\", \"YrSold\", \"FullBath\", \"HalfBath\", \"BedroomAbvGr\", \"Utilities\", \"GarageCars\", \"Fence\"]\nX_columns = [\"YrSold\", \"FullBath\", \"HalfBath\", \"BedroomAbvGr\", \"Utilities\", \"GarageCars\", \"Fence\"]\nreal_data = real_data[scraped_columns]\nX_real_data = real_data[X_columns]\nX_real_data.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"0fa856296d6ae9a452981cfbf3cdd5fdc8f18786"},"cell_type":"code","source":"def X_preprocess(home_data_copy):\n # Separate the categoricals from the numericals\n X_numericals, X_categoricals = separate_X_numerical_categorical(home_data_copy)\n\n ##############\n # Fill X_categoricals nan values with \"nan\"\n X_cat_filled = fill_nan_X_categoricals(X_categoricals)\n\n # On hot encode the categoricals\n X_cat_one_hot = one_hot_encode_categories(X_cat_filled)\n\n\n ##############\n # Fill X_numericals nans with imputing\n X_numericals_filled, my_imputer = fill_nan_X_numericals(X_numericals)\n\n # Transform also the X_numericals values by power transformer\n X_numericals_scaled, X_scaler = X_numericals_transform(X_numericals_filled)\n\n #X_numericals_PCA, X_PCA_transformer = PCA_transform(X_numericals_scaled)\n\n #X_numericals_PCA.describe()\n X_numericals_outlier = X_numericals_scaled\n\n # Calculate the mean of X_numericals_scaled\n X_mean = np.mean(X_numericals_outlier.values, axis=0, keepdims=True)[0]\n\n # Calculate the covariance matrix of X_numericals_scaled\n X_covariance_matrix = np.cov(X_numericals_outlier.values.T)\n\n # Calculate the multivariate normal distribution, which is a matrix filled with measurement probabilities being not outlier.\n MGD = stats.multivariate_normal.pdf(X_numericals_outlier, mean=X_mean, cov=X_covariance_matrix)\n\n # Select a threshold\n threshold = -150\n\n # Calculate a filter array for removing the outliers from the dataset\n filter_array = [True if np.log(meas) > threshold else False for meas in MGD]\n\n X_numericals_outlier_scaled, X_scaler2 = X_numericals_transform(X_numericals_outlier)\n \n # Remove the outliers from the categoricals and numericals\n X_numericals_outlierless = X_numericals_outlier_scaled.iloc[filter_array]\n X_categoricals_outlierless = X_cat_one_hot.iloc[filter_array]\n\n # Concatenate the X_numericals_scaled and X_categoricals\n X_outlierless = pd.concat([X_numericals_outlierless, X_categoricals_outlierless], axis=1)\n print(\"The final shape of the features: {}\".format(X_outlierless.shape))\n \n \n Y_outlierless = y_power[filter_array]\n print(\"The final shape of the output: {}\".format(Y_outlierless.shape))\n\n plt.plot(np.log(np.sort(MGD)))\n plt.plot(np.ones(len(MGD)) * threshold)\n plt.title('Multivariate normal distribution probabilities ordered.')\n plt.xlabel('Datapoints')\n plt.ylabel('Probability')\n plt.show()\n \n return X_outlierless, Y_outlierless, X_cat_one_hot, my_imputer, X_scaler, X_scaler2\n\n\ndef grid_search(iterations):\n hyperparameters = pd.DataFrame()\n\n models = []\n\n # grid search\n for i in range(iterations):\n\n # The hyperparameters:\n\n # The learning rate between 0.333 and 0.00333\n learning_rate = 10 ** (-2 * np.random.rand()) / 3\n\n # The max_depth between 3-12\n max_depth = np.random.randint(3, 12)\n\n # The number of estimators between 1000 and 10000\n n_estimators = np.random.randint(1000, 10000)\n\n # Base score between 1 and 0.5\n base_score = 1 - 0.5 * (np.random.rand())\n\n # Subsample between 1 and 0.6\n subsample = 1 - 0.4 * (np.random.rand())\n\n model = XGBRegressor(max_depth=max_depth,\n learning_rate=learning_rate,\n n_estimators=n_estimators,\n subsample=subsample,\n base_score=base_score,\n random_state=0)\n\n # Calculate the score\n mean_cross_validation = score(model, X, y)\n\n models.append(model)\n\n # Append the hyperparameters for lates examination\n hyperparameters = hyperparameters.append(pd.DataFrame([[mean_cross_validation,\n learning_rate, \n max_depth,\n n_estimators,\n base_score,\n subsample]],\n columns=[\n \"mean_cross_validation\",\n \"learning_rate\",\n \"max_depth\",\n \"n_estimators\",\n \"base_score\",\n \"subsample\"\n ]))\n\n return models, hyperparameters","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"331b717bdca616c2ecb2af89bc24039281f16a25"},"cell_type":"code","source":"X_real_data.info()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"c0133cbdf60aa4133b3ca8e442d055dcc226f482"},"cell_type":"code","source":"home_data[scraped_columns].info()\nhome_data_fenced = home_data_copy.copy()\n\n# The fence feature of the scraped data is only categorical so concentrate the train data\nhome_data_fenced[\"Fence\"] = home_data_copy[\"Fence\"].apply(lambda x: 1 if isinstance(x, str) else 0)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"4cf0bf19623d1fda2ffc88385947ab9ba06d4255"},"cell_type":"code","source":"X, y, X_cat_one_hot, my_imputer, X_scaler, X_scaler2 = X_preprocess(home_data_fenced[X_columns])\nX.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"286a46451c3d9935a53c35fcdd51eaf892ed8e89"},"cell_type":"code","source":"# Teach previous model on new columns\nmodels, hyperparameters = grid_search(20)\nscores = hyperparameters[\"mean_cross_validation\"]\ndisplay(hyperparameters)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"0fc7c77971cdcf529d1421a2a8f84135218a1197"},"cell_type":"code","source":"# Try the default GradientBoostingRegressor\nmodel = GradientBoostingRegressor()\nscores.append(models.append(score(model, X, y)))\nmodels.append(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"644dea726c5beb440569fe866546c18a78021bb8"},"cell_type":"code","source":"# Try the default XGBRegressor\nmodel = XGBRegressor()\nscores.append(models.append(score(model, X, y)))\nmodels.append(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"425f756d47fe1e1321fa3478a3e56042e6e22320"},"cell_type":"code","source":"model = models[np.argmin(scores.values)]\nprint(model)\nmodel.fit(X, y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"89122cf8474b05ea0d09258a76fb65f25b9a8f7b"},"cell_type":"code","source":"X_test = convert_test_data(X_real_data, X_cat_one_hot, my_imputer, X_scaler, X_scaler2)\nX_test = X_test[X.columns]\nX_test.shape, X_real_data.shape\nX_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"fb884f78e05c1a71a2d8e0454fd9bda2ab5273b7"},"cell_type":"code","source":"assert 0 == np.sum(~np.isfinite(X_test.values)), \"There are infinite values in the X_test.\"\nassert 0 == np.sum(np.isnan(X_test.values)), \"There are nan values in the X_test.\"\n\n# Make predictions which we will submit.\ntest_preds_unscaled = model.predict(X_test).reshape(-1, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"e9502dddc8baa8c145a0aead4c312bbb7f2779c5"},"cell_type":"code","source":"# Inverse transform the predictions to the original scale\npredicted_prices = list(y_power_scaler.inverse_transform(test_preds_unscaled)[:,0].astype('int'))\nactual_prices = list(real_data[\"SalePrice\"])\n\npd.DataFrame(np.array([predicted_prices, actual_prices]).T, columns=[\"Predicted prices (USD)\",\n \"Actual prices (USD)\"]).transpose()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"c9413a57de521a6f9f1ae8612d37842158e09a28"},"cell_type":"code","source":"prices_comparison = pd.DataFrame([predicted_prices, actual_prices]).transpose().sort_values(1)\nprices_comparison.columns = [\"predicted\", \"actual\"]\n\nax = prices_comparison.plot(kind=\"bar\", title=\"Comparison real vs predicted house prices.\")\n\nax.set_xlabel('House Id')\nax.set_ylabel('Price (USD)')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"b11390aa6ddce77031f888196f47a570d2a58b94"},"cell_type":"code","source":"prices_comparison.corr()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"468a70b386e20269d28a0f72975a14853ffbec39"},"cell_type":"code","source":"print(\"The final MAE score of the real life data: {}\".format(mean_absolute_error(prices_comparison[\"predicted\"],\n prices_comparison[\"actual\"])))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2c11388a53810d76f5a7ff0a9dff28e557fbdc2e"},"cell_type":"markdown","source":"## Final conclusion\n\nThe final real life guesses are correlating with the actual values and in the same range but there are some things to mention:\nThe houses were not sold yet so the price is actually not the ground truth.\nWe used much less features in the real life example.\nThe real life example and the dataset is close but not from the same distribution meaning not from the same source. We could include part of the real life examples in the train data which would help. Also collecting more real life data would help.\n\nThank you for your attention! I hope you got some good insights. Happy machine learning!\n\nCreated by Andor Hofecker and Tamás Kardos"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}