2 |Mathematical Concepts\n\n\n
2.1 | Linear Algebra
\n
\n
2.2 | Geometric Interpretation of Principal Components
\n
\n
2.3 | Calculating Principal Components
\n
\n
2.4 | Scaling of Variables
\n
\n
2.5 | Component reproducibility
\n
\n
2.6 | Outliers Influence
\n
\n
2.7 | Explained Variance Ratio
\n
\n
2.8 | Optimal Number of Principal Components
\n
\n
3.1 | Feature Engineering
\n
📌 PCA Best Practices. There are a few things to keep in mind when applying PCA:\n\n* PCA only works with numeric features, like continuous quantities or counts.\n* PCA is sensitive to scale. It's good practice to standardize your data before applying PCA, unless you know you have good reason not to.\n* Consider removing or constraining outliers, since they can an have an undue influence on the results. \n
","metadata":{"_uuid":"7746e782-6dde-46e3-a0d3-d6a9ceee5acb","_cell_guid":"6ef2ceea-f541-4f7f-9d57-1c0a3ac81a5d","trusted":true}},{"cell_type":"code","source":"# Setup feedback system\nfrom learntools.core import binder\nbinder.bind(globals())\nfrom learntools.feature_engineering_new.ex5 import *\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom sklearn.decomposition import PCA\nfrom sklearn.feature_selection import mutual_info_regression\nfrom sklearn.model_selection import cross_val_score\nfrom xgboost import XGBRegressor\n\n# Set Matplotlib defaults\nplt.style.use(\"seaborn-whitegrid\")\nplt.rc(\"figure\", autolayout=True)\nplt.rc(\n \"axes\",\n labelweight=\"bold\",\n labelsize=\"large\",\n titleweight=\"bold\",\n titlesize=14,\n titlepad=10,\n)\n\n\ndef apply_pca(X, standardize=True):\n # Standardize\n if standardize:\n X = (X - X.mean(axis=0)) / X.std(axis=0)\n # Create principal components\n pca = PCA()\n X_pca = pca.fit_transform(X)\n # Convert to dataframe\n component_names = [f\"PC{i+1}\" for i in range(X_pca.shape[1])]\n X_pca = pd.DataFrame(X_pca, columns=component_names)\n # Create loadings\n loadings = pd.DataFrame(\n pca.components_.T, # transpose the matrix of loadings\n columns=component_names, # so the columns are the principal components\n index=X.columns, # and the rows are the original features\n )\n return pca, X_pca, loadings\n\n\ndef plot_variance(pca, width=8, dpi=100):\n # Create figure\n fig, axs = plt.subplots(1, 2)\n n = pca.n_components_\n grid = np.arange(1, n + 1)\n # Explained variance\n evr = pca.explained_variance_ratio_\n axs[0].bar(grid, evr)\n axs[0].set(\n xlabel=\"Component\", title=\"% Explained Variance\", ylim=(0.0, 1.0)\n )\n # Cumulative Variance\n cv = np.cumsum(evr)\n axs[1].plot(np.r_[0, grid], np.r_[0, cv], \"o-\")\n axs[1].set(\n xlabel=\"Component\", title=\"% Cumulative Variance\", ylim=(0.0, 1.0)\n )\n # Set up figure\n fig.set(figwidth=8, dpi=100)\n return axs\n\n\ndef make_mi_scores(X, y):\n X = X.copy()\n for colname in X.select_dtypes([\"object\", \"category\"]):\n X[colname], _ = X[colname].factorize()\n # All discrete features should now have integer dtypes\n discrete_features = [pd.api.types.is_integer_dtype(t) for t in X.dtypes]\n mi_scores = mutual_info_regression(X, y, discrete_features=discrete_features, random_state=0)\n mi_scores = pd.Series(mi_scores, name=\"MI Scores\", index=X.columns)\n mi_scores = mi_scores.sort_values(ascending=False)\n return mi_scores\n\n\ndef score_dataset(X, y, model=XGBRegressor()):\n # Label encoding for categoricals\n for colname in X.select_dtypes([\"category\", \"object\"]):\n X[colname], _ = X[colname].factorize()\n # Metric for Housing competition is RMSLE (Root Mean Squared Log Error)\n score = cross_val_score(\n model, X, y, cv=5, scoring=\"neg_mean_squared_log_error\",\n )\n score = -1 * score.mean()\n score = np.sqrt(score)\n return score\n\n\ndf = pd.read_csv(\"../input/fe-course-data/ames.csv\")","metadata":{"_uuid":"35494212-a9cc-4cbb-b740-8e31bc4cbc92","_cell_guid":"b2dfb909-a4ca-4e6e-8498-b5e4d450f554","collapsed":false,"jupyter":{"outputs_hidden":false},"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-05T10:33:57.827751Z","iopub.execute_input":"2022-07-05T10:33:57.82823Z","iopub.status.idle":"2022-07-05T10:33:59.821492Z","shell.execute_reply.started":"2022-07-05T10:33:57.828146Z","shell.execute_reply":"2022-07-05T10:33:59.820706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's choose a few features that are highly correlated with our target, `SalePrice`.","metadata":{"_uuid":"e3f48a39-5170-4591-8c4b-d7e867e799cf","_cell_guid":"71ab37ba-712b-4b26-acfe-2142bbc9b86b","trusted":true}},{"cell_type":"code","source":"features = [\n \"GarageArea\",\n \"YearRemodAdd\",\n \"TotalBsmtSF\",\n \"GrLivArea\",\n]\n\nprint(\"Correlation with SalePrice:\\n\")\nprint(df[features].corrwith(df.SalePrice))","metadata":{"_uuid":"c7820d00-82ba-462a-a35b-5a6109e41804","_cell_guid":"e55fc760-6ab3-4684-926c-72f229ca5113","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:33:59.823114Z","iopub.execute_input":"2022-07-05T10:33:59.823633Z","iopub.status.idle":"2022-07-05T10:33:59.848053Z","shell.execute_reply.started":"2022-07-05T10:33:59.823597Z","shell.execute_reply":"2022-07-05T10:33:59.847401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We'll rely on PCA to untangle the correlational structure of these features and suggest relationships that might be usefully modeled with new features. Run this cell to apply PCA and extract the loadings.","metadata":{"_uuid":"2e4e1a44-ef6a-48fe-823f-bf630243ede2","_cell_guid":"f3348b07-7ccb-4ba1-b1a4-732e3c15bd82","trusted":true}},{"cell_type":"code","source":"X = df.copy()\ny = X.pop(\"SalePrice\")\nX = X.loc[:, features]\n\n# `apply_pca`, defined above, reproduces the code from the tutorial\npca, X_pca, loadings = apply_pca(X)\nprint(loadings)","metadata":{"_uuid":"3ae6474a-0a94-469f-bc0e-888221501544","_cell_guid":"a3df0ef6-c13f-4af5-8f5b-b07b4bd6e80c","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:33:59.849513Z","iopub.execute_input":"2022-07-05T10:33:59.849996Z","iopub.status.idle":"2022-07-05T10:33:59.895832Z","shell.execute_reply.started":"2022-07-05T10:33:59.849959Z","shell.execute_reply":"2022-07-05T10:33:59.894536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_variance(pca)","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:33:59.897828Z","iopub.execute_input":"2022-07-05T10:33:59.898173Z","iopub.status.idle":"2022-07-05T10:34:00.361192Z","shell.execute_reply.started":"2022-07-05T10:33:59.898122Z","shell.execute_reply":"2022-07-05T10:34:00.360108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.1.1 | Interpret Component Loadings\n\nLook at the loadings for components `PC1` and `PC3`. Can you think of a description of what kind of contrast each component has captured? After you've thought about it, run the next cell for a solution.","metadata":{"_uuid":"d3a5fc97-3d80-4e8a-ab18-f2e2fa2b9507","_cell_guid":"dbbba74d-2cd6-4348-8022-a02d0aa696fc","trusted":true}},{"cell_type":"code","source":"# View the solution (Run this cell to receive credit!)\nq_1.check()","metadata":{"_uuid":"a646c0a5-a099-4a4d-bade-03ab38ba7c42","_cell_guid":"270523d8-26bf-40fe-b0a7-d775e3c57497","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:34:00.36263Z","iopub.execute_input":"2022-07-05T10:34:00.363141Z","iopub.status.idle":"2022-07-05T10:34:00.374227Z","shell.execute_reply.started":"2022-07-05T10:34:00.363097Z","shell.execute_reply":"2022-07-05T10:34:00.373248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-------------------------------------------------------------------------------\n\nYour goal in this question is to use the results of PCA to discover one or more new features that improve the performance of your model. One option is to create features inspired by the loadings, like we did in the tutorial. Another option is to use the components themselves as features (that is, add one or more columns of `X_pca` to `X`).\n\n### 3.1.2 | Create New Features\n\nAdd one or more new features to the dataset `X`. For a correct solution, get a validation score below 0.140 RMSLE. (If you get stuck, feel free to use the `hint` below!). For this first option we'll create some new features based on what PCA tells about the relations between previous features. ","metadata":{"_uuid":"3ae129fb-a01d-4277-b437-82cf4121c270","_cell_guid":"a331c2e1-bfe3-414a-bd01-be6c184b6b68","trusted":true}},{"cell_type":"code","source":"X = df.copy()\ny = X.pop(\"SalePrice\")\n\n# YOUR CODE HERE: Add new features to X.\nX[\"Feature1\"] = X.GrLivArea + X.TotalBsmtSF\nX[\"Feature2\"] = X.YearRemodAdd * X.TotalBsmtSF\n\nscore = score_dataset(X, y)\nprint(f\"Your score: {score:.5f} RMSLE\")\n\n# Check your answer\nq_2.check()","metadata":{"_uuid":"f85c812e-dde4-4f9c-ba10-6c8a3c5da4bf","_cell_guid":"0533bf18-b228-4e87-8577-971e1cb4115f","collapsed":false,"lines_to_next_cell":0,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:34:00.37561Z","iopub.execute_input":"2022-07-05T10:34:00.376939Z","iopub.status.idle":"2022-07-05T10:34:08.03073Z","shell.execute_reply.started":"2022-07-05T10:34:00.376893Z","shell.execute_reply":"2022-07-05T10:34:08.029818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Lines below will give you a hint or solution code\n#q_2.hint()\n#q_2.solution()","metadata":{"_uuid":"216ef7db-b90b-4f53-9024-ffcc55373721","_cell_guid":"85eef27e-89b6-4679-9fdd-554a25da4675","collapsed":false,"lines_to_next_cell":0,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:34:08.031989Z","iopub.execute_input":"2022-07-05T10:34:08.032305Z","iopub.status.idle":"2022-07-05T10:34:08.035954Z","shell.execute_reply.started":"2022-07-05T10:34:08.032273Z","shell.execute_reply":"2022-07-05T10:34:08.034841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = df.copy()\nX[\"Feature1\"] = X.GrLivArea + X.TotalBsmtSF\nX[\"Feature2\"] = X.YearRemodAdd * X.TotalBsmtSF\nsns.regplot(x=\"Feature1\", y='SalePrice', data=X);","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:08.037343Z","iopub.execute_input":"2022-07-05T10:34:08.037702Z","iopub.status.idle":"2022-07-05T10:34:08.498295Z","shell.execute_reply.started":"2022-07-05T10:34:08.037666Z","shell.execute_reply":"2022-07-05T10:34:08.49743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X[\"Feature2\"] = X.YearRemodAdd * X.TotalBsmtSF\nsns.regplot(x=\"Feature2\", y='SalePrice', data=X);","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:08.499659Z","iopub.execute_input":"2022-07-05T10:34:08.499919Z","iopub.status.idle":"2022-07-05T10:34:08.933395Z","shell.execute_reply.started":"2022-07-05T10:34:08.499888Z","shell.execute_reply":"2022-07-05T10:34:08.932841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's go on with 2nd option. We'll use now principal components as features. ","metadata":{}},{"cell_type":"code","source":"# Solution 2: Uses components\nX = X.join(X_pca)\ny = X.pop(\"SalePrice\")\n\nscore = score_dataset(X, y)\nprint(f\"Your score: {score:.5f} RMSLE\")","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:08.935291Z","iopub.execute_input":"2022-07-05T10:34:08.935696Z","iopub.status.idle":"2022-07-05T10:34:14.349719Z","shell.execute_reply.started":"2022-07-05T10:34:08.93566Z","shell.execute_reply":"2022-07-05T10:34:14.349031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-------------------------------------------------------------------------------\n\nThe next question explores a way you can use PCA to detect outliers in the dataset (meaning, data points that are unusually extreme in some way). Outliers can have a detrimental effect on model performance, so it's good to be aware of them in case you need to take corrective action. PCA in particular can show you anomalous *variation* which might not be apparent from the original features: neither small houses nor houses with large basements are unusual, but it is unusual for small houses to have large basements. That's the kind of thing a principal component can show you.\n\nRun the next cell to show distribution plots for each of the principal components you created above.","metadata":{"_uuid":"c344d2a9-8afb-4be8-9eb1-4f96e19642ab","_cell_guid":"758cb23d-70c9-46da-8b55-6caf8c14545e","trusted":true}},{"cell_type":"code","source":"sns.catplot(\n y=\"value\",\n col=\"variable\",\n data=X_pca.melt(),\n kind='boxen',\n sharey=False,\n col_wrap=2,\n);","metadata":{"_uuid":"db73b492-bb95-4630-b154-44846c78e3de","_cell_guid":"f8196c57-a18b-41c2-a671-23dbf7847465","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:34:14.350798Z","iopub.execute_input":"2022-07-05T10:34:14.351053Z","iopub.status.idle":"2022-07-05T10:34:15.012505Z","shell.execute_reply.started":"2022-07-05T10:34:14.351023Z","shell.execute_reply":"2022-07-05T10:34:15.011521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see, in each of the components there are several points lying at the extreme ends of the distributions -- outliers, that is.\n\nNow run the next cell to see those houses that sit at the extremes of a component:","metadata":{"_uuid":"8bb51c20-6561-472a-a890-d1be986fe805","_cell_guid":"bfe1ac56-fdfd-48df-9f91-aa00dddd4aca","trusted":true}},{"cell_type":"code","source":"# You can change PC1 to PC2, PC3, or PC4\ncomponent = \"PC1\"\n\nidx = X_pca[component].sort_values(ascending=False).index\ndf.loc[idx, [\"SalePrice\", \"Neighborhood\", \"SaleCondition\"] + features]","metadata":{"_uuid":"b8ab8659-20e0-428f-b847-8a030aa63af1","_cell_guid":"c9fd2e4a-a4dc-4864-971a-b828a285a966","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:34:15.013682Z","iopub.execute_input":"2022-07-05T10:34:15.013911Z","iopub.status.idle":"2022-07-05T10:34:15.038621Z","shell.execute_reply.started":"2022-07-05T10:34:15.013879Z","shell.execute_reply":"2022-07-05T10:34:15.037709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.1.3 | Outlier Detection\n\nDo you notice any patterns in the extreme values? Does it seem like the outliers are coming from some special subset of the data? After you've thought about your answer, run the next cell for the solution and some discussion.","metadata":{"_uuid":"c0ed9f55-050f-4325-98f2-823c289a60f8","_cell_guid":"443893bd-d0b4-4f9b-84ae-3d3f43836ef7","trusted":true}},{"cell_type":"code","source":"# View the solution (Run this cell to receive credit!)\nq_3.check()","metadata":{"_uuid":"c1d8a8dd-5d64-4c8b-9992-76bfb9f6d350","_cell_guid":"9e062e1d-e6b7-4a32-baa6-31b3ab3f2c0f","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2022-07-05T10:34:15.039957Z","iopub.execute_input":"2022-07-05T10:34:15.040172Z","iopub.status.idle":"2022-07-05T10:34:15.047854Z","shell.execute_reply.started":"2022-07-05T10:34:15.040147Z","shell.execute_reply":"2022-07-05T10:34:15.047009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In addition, I'll include how I always try to detect outliers. As it's well known, an outlier can be of two types: Univariate and Multivariate. These outliers can be found when we look at distribution of a single variable. Multi-variate outliers are outliers in an n-dimensional space. We'll start by detecting whether there are univariate outliers in our dataset or not. \n\n### 3.1.4 | Univariate Outliers\n\n#### 3.1.4.1 | Grubbs Test\n\n$$\n\\begin{array}{l}{\\text { Grubbs' test is defined for the hypothesis: }} \\\\ {\\begin{array}{ll}{\\text { Ho: }} {\\text { There are no outliers in the data set }} \\\\ {\\mathrm{H}_{\\mathrm{1}} :} {\\text { There is exactly one outlier in the data set }}\\end{array}}\\end{array}\n$$\n$$\n\\begin{array}{l}{\\text {The Grubbs' test statistic is defined as: }} \\\\ {\\qquad G_{calculated}=\\frac{\\max \\left|X_{i}-\\overline{X}\\right|}{SD}} \\\\ {\\text { with } \\overline{X} \\text { and } SD \\text { denoting the sample mean and standard deviation, respectively. }} \\end{array}\n$$\n$$\nG_{critical}=\\frac{(N-1)}{\\sqrt{N}} \\sqrt{\\frac{\\left(t_{\\alpha /(2 N), N-2}\\right)^{2}}{N-2+\\left(t_{\\alpha /(2 N), N-2}\\right)^{2}}}\n$$\n$$\n\\begin{array}{l}{\\text { If the calculated value is greater than critical, you can reject the null hypothesis and conclude that one of the values is an outlier }}\\end{array}$$","metadata":{}},{"cell_type":"code","source":"import scipy.stats as stats\ndef grubbs_test(x, feature):\n n = len(x)\n mean_x = np.mean(x)\n sd_x = np.std(x)\n numerator = max(abs(x-mean_x))\n g_calculated = numerator/sd_x\n print(\"Feature:\", feature)\n print(\"Grubbs Calculated Value:\",g_calculated)\n t_value = stats.t.ppf(1 - 0.05 / (2 * n), n - 2)\n g_critical = ((n - 1) * np.sqrt(np.square(t_value))) / (np.sqrt(n) * np.sqrt(n - 2 + np.square(t_value)))\n print(\"Grubbs Critical Value:\",g_critical)\n if g_critical > g_calculated:\n print(\"From grubbs_test we observe that calculated value is lesser than critical value, Accept null hypothesis and conclude that there is no outliers\\n\")\n else:\n print(\"From grubbs_test we observe that calculated value is greater than critical value, Reject null hypothesis and conclude that there is an outliers\\n\")\n \n print(\"==============================================================================================================================================\")\n \ngrubbs_test(X['PC1'], 'PC1')\ngrubbs_test(X['PC2'], 'PC2')\ngrubbs_test(X['PC3'], 'PC3')\ngrubbs_test(X['PC4'], 'PC4')","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:15.04931Z","iopub.execute_input":"2022-07-05T10:34:15.049837Z","iopub.status.idle":"2022-07-05T10:34:15.201488Z","shell.execute_reply.started":"2022-07-05T10:34:15.049801Z","shell.execute_reply":"2022-07-05T10:34:15.200548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.1.4.2 | Z-score method\n\nUsing Z score method,we can find out how many standard deviations value away from the mean. \n\n\n\n Figure in the left shows area under normal curve and how much area that standard deviation covers.\n* 68% of the data points lie between + or - 1 standard deviation.\n* 95% of the data points lie between + or - 2 standard deviation\n* 99.7% of the data points lie between + or - 3 standard deviation\n\n$\\begin{array}{l} {R.Z.score=\\frac{0.6745*( X_{i} - Median)}{MAD}} \\end{array}$\n\nIf the z score of a data point is more than 3 (because it cover 99.7% of area), it indicates that the data value is quite different from the other values. It is taken as outliers.","metadata":{}},{"cell_type":"code","source":"out=[]\ndef Zscore_outlier(df):\n m = np.mean(df)\n sd = np.std(df)\n row = 0\n for i in df: \n z = (i-m)/sd\n if np.abs(z) > 3: \n out.append(row)\n row += 1\n return out\n\nZscore_outlier(X['PC1'])","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:15.202854Z","iopub.execute_input":"2022-07-05T10:34:15.203599Z","iopub.status.idle":"2022-07-05T10:34:15.21844Z","shell.execute_reply.started":"2022-07-05T10:34:15.203551Z","shell.execute_reply":"2022-07-05T10:34:15.217785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.1.4.3 | Isolation Forest\n\n[Isolation Forest](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.IsolationForest.html#sklearn.ensemble.IsolationForest) is implemented in scikit-learn. \nIsolation forest is an algorithm to detect outliers. It partitions the data using a set of trees and provides an anomaly score looking at how isolated the point is in the structure found. The anomaly score is then used to tell apart outliers from normal observations.\nAn important concept in this method is the isolation number. The isolation number is the number of splits needed to isolate a data point. This number of splits is ascertained by following these steps:\n\n* A point “a” to isolate is selected randomly.\n* A random data point “b” is selected that is between the minimum and maximum value and different from “a”.\n* If the value of “b” is lower than the value of “a”, the value of “b” becomes the new lower limit.\n* If the value of “b” is greater than the value of “a”, the value of “b” becomes the new upper limit.\n* This procedure is repeated as long as there are data points other than “a” between the upper and the lower limit.\n\nIt requires fewer splits to isolate an outlier than it does to isolate a non-outlier, i.e. an outlier has a lower isolation number in comparison to a non-outlier point. A data point is therefore defined as an outlier if its isolation number is lower than the threshold. The threshold is defined based on the estimated percentage of outliers in the data, which is the starting point of this outlier detection algorithm.\n\n","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import IsolationForest\nfig, axs = plt.subplots(2, 3, figsize=(22, 12), facecolor='w', edgecolor='k')\naxs = axs.ravel()\n\nfor i, column in enumerate(X[['PC1','PC2','PC3','PC3','GrLivArea','TotalBsmtSF']].columns):\n isolation_forest = IsolationForest(contamination='auto')\n isolation_forest.fit(X[column].values.reshape(-1,1))\n\n xx = np.linspace(X[column].min(), X[column].max(), len(X)).reshape(-1,1)\n anomaly_score = isolation_forest.decision_function(xx)\n outlier = isolation_forest.predict(xx)\n \n axs[i].plot(xx, anomaly_score, label='anomaly score')\n axs[i].fill_between(xx.T[0], np.min(anomaly_score), np.max(anomaly_score), \n where=outlier==-1, color='r', \n alpha=.4, label='outlier region')\n axs[i].legend()\n axs[i].set_title(column)","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:15.219752Z","iopub.execute_input":"2022-07-05T10:34:15.220704Z","iopub.status.idle":"2022-07-05T10:34:18.388541Z","shell.execute_reply.started":"2022-07-05T10:34:15.220667Z","shell.execute_reply":"2022-07-05T10:34:18.387331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's analyse for example how much ratio of outliers are in `TotalBsmtSF` feature. ","metadata":{}},{"cell_type":"code","source":"X[(X['TotalBsmtSF'] > 2000) | (X['TotalBsmtSF'] < 650)].shape[0] / X.shape[0]","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:18.389929Z","iopub.execute_input":"2022-07-05T10:34:18.390415Z","iopub.status.idle":"2022-07-05T10:34:18.399287Z","shell.execute_reply.started":"2022-07-05T10:34:18.390381Z","shell.execute_reply":"2022-07-05T10:34:18.39814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.1.5 | Multivariate Outliers\n\nI will show multivariate outlier detection using two scikit-learn methods: For normally distributed data, use EllipticEnvelope.\n\n#### 3.1.5.1 | EllipticEnvelope","metadata":{}},{"cell_type":"code","source":"from sklearn.covariance import EllipticEnvelope\nfrom sklearn.cluster import DBSCAN\nclf = EllipticEnvelope()\n\ncols = ['PC1','PC2']\nX2 = X[cols].copy()\n\nlegend = {}\nxx, yy = np.meshgrid(np.linspace(-35, 35, 500), np.linspace(-35, 35, 500))\nplt.figure(1, figsize=(10,10))\nclf.fit(X2.values[:,:2])\nZ = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])\nZ = Z.reshape(xx.shape)\nlegend['EllipticEnvelope'] = plt.contour(\n xx, yy, Z, levels=[0], linewidths=2, colors=['m'])\n\nlegend_values_list = list(legend.values())\nlegend_keys_list = list(legend.keys())\n\nplt.figure(1, figsize=(10,10))# two clusters\nplt.title(\"Outlier detection on first two columns of dummydf\")\nplt.scatter(X2.values[:, 0], X.values[:, 1], color='black')\nbbox_args = dict(boxstyle=\"round\", fc=\"0.8\")\narrow_args = dict(arrowstyle=\"->\")\n\nplt.xlim((xx.min(), xx.max()))\nplt.ylim((yy.min(), yy.max()))\nplt.legend(legend_values_list[0].collections, legend_keys_list,\n loc=\"upper center\");\nplt.ylabel(X2.columns[1]);\nplt.xlabel(X2.columns[0]);","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:18.400813Z","iopub.execute_input":"2022-07-05T10:34:18.401051Z","iopub.status.idle":"2022-07-05T10:34:19.986621Z","shell.execute_reply.started":"2022-07-05T10:34:18.401025Z","shell.execute_reply":"2022-07-05T10:34:19.985625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.1.5.2 | DBSCAN (Density-Based Spatial Clustering of Applications with Noise)\n\nThis is a clustering algorithm (an alternative to K-Means) that clusters points together and identifies any points not belonging to a cluster as outliers. It’s like K-means, except the number of clusters does not need to be specified in advance.\nI will show you an example of using DBScan but before we start, let’s cover some important concepts. DBScan has three important concepts:\n\n- Core Points: In order to understand the concept of the core points, we need to visit some of the hyperparameters used to define DBScan job. First hyperparameter (HP) is min_samples. This is simply the minimum number of core points needed in order to form a cluster. second important HP is 'eps'. 'eps' is the maximum distance between two samples for them to be considered as in the same cluster.\n- Border Points are in the same cluster as core points but much further away from the centre of the cluster.\n- Everything else is called Noise Points, those are data points that do not belong to any cluster. They can be anomalous or non-anomalous and they need further investigation. Now, let’s see some code.\n\n","metadata":{}},{"cell_type":"code","source":"from sklearn.cluster import DBSCAN\nfrom sklearn.preprocessing import StandardScaler\n# scale data first\nX = StandardScaler().fit_transform(X2[cols].values)\n\ndb = DBSCAN(eps=3.0, min_samples=10).fit(X)\nlabels = db.labels_\n\npd.Series(labels).value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:19.988081Z","iopub.execute_input":"2022-07-05T10:34:19.988405Z","iopub.status.idle":"2022-07-05T10:34:20.169746Z","shell.execute_reply.started":"2022-07-05T10:34:19.988364Z","shell.execute_reply":"2022-07-05T10:34:20.168947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Visualize the outlier samples in the context of the first two features, Rooms and Price","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\n\nunique_labels = set(labels)\ncolors = ['blue', 'red']\n\nfor color,label in zip(colors, unique_labels):\n sample_mask = [True if l == label else False for l in labels]\n plt.plot(X[:,0][sample_mask], X[:, 1][sample_mask], 'o', color=color);\nplt.xlabel(X2.columns[0]);\nplt.ylabel(X2.columns[1]);","metadata":{"execution":{"iopub.status.busy":"2022-07-05T10:34:20.170818Z","iopub.execute_input":"2022-07-05T10:34:20.171029Z","iopub.status.idle":"2022-07-05T10:34:20.44717Z","shell.execute_reply.started":"2022-07-05T10:34:20.171004Z","shell.execute_reply":"2022-07-05T10:34:20.446286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4 |References\n\n* [Análisis de Componentes Principales](https://www.cienciadedatos.net/documentos/35_principal_component_analysis)\n* [Kaggle PCA Module, from Feature Engineering Course](https://www.kaggle.com/learn/feature-engineering)\n* [Outlier Detection Techniques: Simplified](https://www.kaggle.com/code/rpsuraj/outlier-detection-techniques-simplified#1)-DBSCAN-(Density-Based-Spatial-Clustering-of-Applications-with-Noise):)\n* [Outlier Detection Practice: uni/multivariate](https://www.kaggle.com/code/kevinarvai/outlier-detection-practice-uni-multivariate#Table-of-Contents)","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}