123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442 |
- # -*- coding: utf-8 -*-
- # This code is initially based on the Kaggle kernel from Sergei Neviadomski, which can be found in the following link
- # https://www.kaggle.com/neviadomski/how-to-get-to-top-25-with-simple-model-sklearn/notebook
- # and the Kaggle kernel from Pedro Marcelino, which can be found in the link below
- # https://www.kaggle.com/pmarcelino/comprehensive-data-exploration-with-python/notebook
- # Also, part of the preprocessing and modelling has been inspired by this kernel from Serigne
- # https://www.kaggle.com/serigne/stacked-regressions-top-4-on-leaderboard
- # And this kernel from juliencs has been pretty helpful too!
- # https://www.kaggle.com/juliencs/a-study-on-regression-applied-to-the-ames-dataset
- # Adding needed libraries and reading data
- import pandas as pd
- import numpy as np
- import matplotlib.pyplot as plt
- import seaborn as sns
- from sklearn import ensemble, tree, linear_model, preprocessing
- from sklearn.preprocessing import LabelEncoder, RobustScaler
- from sklearn.linear_model import ElasticNet, Lasso
- from sklearn.kernel_ridge import KernelRidge
- from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone
- from sklearn.model_selection import KFold, train_test_split, cross_val_score
- from sklearn.metrics import r2_score, mean_squared_error
- from sklearn.pipeline import make_pipeline
- from sklearn.utils import shuffle
- from scipy import stats
- from scipy.stats import norm, skew, boxcox
- from scipy.special import boxcox1p
- import xgboost as xgb
- import lightgbm as lgb
- import warnings
- warnings.filterwarnings('ignore')
- # Class AveragingModels
- # This class is Serigne's simplest way of stacking the prediction models, by
- # averaging them. We are going to use it as it represents the same that we have
- # been using in the late submissions, but this applies perfectly to rmsle_cv function.
- class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):
- def __init__(self, models):
- self.models = models
- # we define clones of the original models to fit the data in
- def fit(self, X, y):
- self.models_ = [clone(x) for x in self.models]
- # Train cloned base models
- for model in self.models_:
- model.fit(X, y)
- return self
- #Now we do the predictions for cloned models and average them
- def predict(self, X):
- predictions = np.column_stack([
- model.predict(X) for model in self.models_
- ])
- return np.mean(predictions, axis=1)
- train = pd.read_csv("../../train.csv")
- test = pd.read_csv("../../test.csv")
- #Save the 'Id' column
- train_ID = train['Id']
- test_ID = test['Id']
- train.drop('Id', axis=1, inplace=True)
- test.drop('Id', axis=1, inplace=True)
- # Visualizing outliers
- fig, ax = plt.subplots()
- ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])
- plt.ylabel('SalePrice', fontsize=13)
- plt.xlabel('GrLivArea', fontsize=13)
- #plt.show()
- # Now the outliers can be deleted
- train = train.drop(train[(train['GrLivArea'] > 4000) & (train['SalePrice'] < 300000)].index)
- #Check the graphic again, making sure there are no outliers left
- fig, ax = plt.subplots()
- ax.scatter(train['GrLivArea'], train['SalePrice'])
- plt.ylabel('SalePrice', fontsize=13)
- plt.xlabel('GrLivArea', fontsize=13)
- #plt.show()
- #We use the numpy fuction log1p which applies log(1+x) to all elements of the column
- train["SalePrice"] = np.log1p(train["SalePrice"])
- #Check the new distribution
- sns.distplot(train['SalePrice'] , fit=norm);
- # Get the fitted parameters used by the function
- (mu, sigma) = norm.fit(train['SalePrice'])
- print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))
- #Now plot the distribution
- plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
- loc='best')
- plt.ylabel('Frequency')
- plt.title('SalePrice distribution')
- #Get also the QQ-plot
- fig = plt.figure()
- res = stats.probplot(train['SalePrice'], plot=plt)
- #plt.show()
- # Splitting to features and labels
- train_labels = train.pop('SalePrice')
- # Test set does not even have a 'SalePrice' column, so both sets can be concatenated
- features = pd.concat([train, test], keys=['train', 'test'])
- # Checking for missing data, showing every variable with at least one missing value in train set
- total_missing_data = features.isnull().sum().sort_values(ascending=False)
- missing_data_percent = (features.isnull().sum()/features.isnull().count()).sort_values(ascending=False)
- missing_data = pd.concat([total_missing_data, missing_data_percent], axis=1, keys=['Total', 'Percent'])
- print(missing_data[missing_data['Percent']> 0])
- # Deleting non-interesting variables for this case study
- features.drop(['Utilities'], axis=1, inplace=True)
- # Imputing missing values and transforming certain columns
- # Converting OverallCond to str
- features.OverallCond = features.OverallCond.astype(str)
- # MSSubClass as str
- features['MSSubClass'] = features['MSSubClass'].astype(str)
- # MSZoning NA in pred. filling with most popular values
- features['MSZoning'] = features['MSZoning'].fillna(features['MSZoning'].mode()[0])
- # LotFrontage NA filling with median according to its OverallQual value
- median = features.groupby('OverallQual')['LotFrontage'].transform('median')
- features['LotFrontage'] = features['LotFrontage'].fillna(median)
- # Alley NA in all. NA means no access
- features['Alley'] = features['Alley'].fillna('NoAccess')
- # MasVnrArea NA filling with median according to its OverallQual value
- median = features.groupby('OverallQual')['MasVnrArea'].transform('median')
- features['MasVnrArea'] = features['MasVnrArea'].fillna(median)
- # MasVnrType NA in all. filling with most popular values
- features['MasVnrType'] = features['MasVnrType'].fillna(features['MasVnrType'].mode()[0])
- # BsmtQual, BsmtCond, BsmtExposure, BsmtFinType1, BsmtFinType2
- # NA in all. NA means No basement
- for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
- features[col] = features[col].fillna('NoBsmt')
- # TotalBsmtSF NA in pred. I suppose NA means 0
- features['TotalBsmtSF'] = features['TotalBsmtSF'].fillna(0)
- # Electrical NA in pred. filling with most popular values
- features['Electrical'] = features['Electrical'].fillna(features['Electrical'].mode()[0])
- # KitchenAbvGr to categorical
- features['KitchenAbvGr'] = features['KitchenAbvGr'].astype(str)
- # KitchenQual NA in pred. filling with most popular values
- features['KitchenQual'] = features['KitchenQual'].fillna(features['KitchenQual'].mode()[0])
- # FireplaceQu NA in all. NA means No Fireplace
- features['FireplaceQu'] = features['FireplaceQu'].fillna('NoFp')
- # Garage-like features NA in all. NA means No Garage
- for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageYrBlt', 'GarageCond'):
- features[col] = features[col].fillna('NoGrg')
- # GarageCars and GarageArea NA in pred. I suppose NA means 0
- for col in ('GarageCars', 'GarageArea'):
- features[col] = features[col].fillna(0.0)
- # SaleType NA in pred. filling with most popular values
- features['SaleType'] = features['SaleType'].fillna(features['SaleType'].mode()[0])
- # PoolQC NA in all. NA means No Pool
- features['PoolQC'] = features['PoolQC'].fillna('NoPool')
- # MiscFeature NA in all. NA means None
- features['MiscFeature'] = features['MiscFeature'].fillna('None')
- # Fence NA in all. NA means no fence
- features['Fence'] = features['Fence'].fillna('NoFence')
- # BsmtHalfBath and BsmtFullBath NA means 0
- for col in ('BsmtHalfBath', 'BsmtFullBath'):
- features[col] = features[col].fillna(0)
- # Functional NA means Typ
- features['Functional'] = features['Functional'].fillna('Typ')
- # NA in Bsmt SF variables means not that type of basement, 0 square feet
- for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF', 'TotalBsmtSF'):
- features[col] = features[col].fillna(0)
- # NA in Exterior1st filled with the most common value
- features['Exterior1st'] = features['Exterior1st'].fillna(features['Exterior1st'].mode()[0])
- # NA in Exterior2nd means No 2nd material
- features['Exterior2nd'] = features['Exterior2nd'].fillna('NoExt2nd')
- # Year and Month to categorical
- features['YrSold'] = features['YrSold'].astype(str)
- features['MoSold'] = features['MoSold'].astype(str)
- # Adding total sqfootage feature and removing Basement, 1st and 2nd floor features
- features['TotalSF'] = features['TotalBsmtSF'] + features['1stFlrSF'] + features['2ndFlrSF']
- #features.drop(['TotalBsmtSF', '1stFlrSF', '2ndFlrSF'], axis=1, inplace=True)
- # Let's rank those categorical features that can be understood to have an order
- # Criterion: give higher ranking to better feature values
- features = features.replace({'Street' : {'Grvl':1, 'Pave':2},
- 'Alley' : {'NoAccess':0, 'Grvl':1, 'Pave':2},
- 'LotShape' : {'I33':1, 'IR2':2, 'IR1':3, 'Reg':4},
- 'LandContour' : {'Low':1, 'HLS':2, 'Bnk':3, 'Lvl':4},
- 'LotConfig' : {'FR3':1, 'FR2':2, 'CulDSac':3, 'Corner':4, 'Inside':5},
- 'LandSlope' : {'Gtl':1, 'Mod':2, 'Sev':3},
- 'HouseStyle' : {'1Story':1, '1.5Fin':2, '1.5Unf':3, '2Story':4, '2.5Fin':5, '2.5Unf':6, 'SFoyer':7, 'SLvl':8},
- 'ExterQual' : {'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'ExterCond' : {'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'BsmtQual' : {'NoBsmt':0, 'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'BsmtCond' : {'NoBsmt':0, 'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'BsmtExposure' : {'NoBsmt':0, 'No':1, 'Mn':2, 'Av':3, 'Gd':4},
- 'BsmtFinType1' : {'NoBsmt':0, 'Unf':1, 'LwQ':2, 'BLQ':3, 'Rec':4, 'ALQ':5, 'GLQ':6},
- 'BsmtFinType2' : {'NoBsmt':0, 'Unf':1, 'LwQ':2, 'BLQ':3, 'Rec':4, 'ALQ':5, 'GLQ':6},
- 'HeatingQC' : {'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'CentralAir' : {'N':0, 'Y':1},
- 'KitchenQual' : {'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'Functional' : {'Sal':0, 'Sev':1, 'Maj2':2, 'Maj1':3, 'Mod':4, 'Min2':5, 'Min1':6, 'Typ':7},
- 'FireplaceQu' : {'NoFp':0, 'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'GarageType' : {'NoGrg':0, 'Detchd':1, 'CarPort':2, 'BuiltIn':3, 'Basment':4, 'Attchd':5, '2Types':6},
- 'GarageFinish' : {'NoGrg':0, 'Unf':1, 'RFn':2, 'Fin':3},
- 'GarageQual' : {'NoGrg':0, 'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'GarageCond' : {'NoGrg':0, 'Po':1, 'Fa':2, 'TA':3, 'Gd':4, 'Ex':5},
- 'PavedDrive' : {'N':0, 'P':1, 'Y':2},
- 'PoolQC' : {'NoPool':0, 'Fa':1, 'TA':2, 'Gd':3, 'Ex':4},
- 'Fence' : {'NoFence':0, 'MnWw':1, 'MnPrv':2, 'GdWo':3, 'GdPrv':4}
- })
- ###################################################################################################
- # Now we try to simplify some of the ranked features, reducing its number of values
- features['SimplifiedOverallCond'] = features.OverallCond.replace({1:1, 2:1, 3:1, # bad
- 4:2, 5:2, 6:2, # average
- 7:3, 8:3, # good
- 9:4, 10:4 # excellent
- })
- features['SimplifiedHouseStyle'] = features.HouseStyle.replace({1:1, # 1 storey houses
- 2:2, 3:2, # 1.5 storey houses
- 4:3, # 2 storey houses
- 5:4, 6:4, # 2.5 storey houses
- 7:5, 8:5 # splitted houses
- })
- features['SimplifiedExterQual'] = features.ExterQual.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedExterCond'] = features.ExterCond.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedBsmtQual'] = features.BsmtQual.replace({1:1, 2:1, # bad, not necessary to check 0 value because will remain 0
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedBsmtCond'] = features.BsmtCond.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedBsmtExposure'] = features.BsmtExposure.replace({1:1, 2:1, # bad
- 3:2, 4:2 # good
- })
- features['SimplifiedBsmtFinType1'] = features.BsmtFinType1.replace({1:1, 2:1, 3:1, # bad
- 4:2, 5:2, # average
- 6:3 # good
- })
- features['SimplifiedBsmtFinType2'] = features.BsmtFinType2.replace({1:1, 2:1, 3:1, # bad
- 4:2, 5:2, # average
- 6:3 # good
- })
- features['SimplifiedHeatingQC'] = features.HeatingQC.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedKitchenQual'] = features.KitchenQual.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedFunctional'] = features.Functional.replace({0:0, 1:0, # really bad
- 2:1, 3:1, # quite bad
- 4:2, 5:2, 6:2, # small deductions
- 7:3 # working fine
- })
- features['SimplifiedFireplaceQu'] = features.FireplaceQu.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedGarageQual'] = features.GarageQual.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedGarageCond'] = features.GarageCond.replace({1:1, 2:1, # bad
- 3:2, 4:2, # good/average
- 5:3 # excellent
- })
- features['SimplifiedPoolQC'] = features.PoolQC.replace({1:1, 2:1, # average
- 3:2, 4:2 # good
- })
- features['SimplifiedFence'] = features.Fence.replace({1:1, 2:1, # bad
- 3:2, 4:2 # good
- })
- # Now, let's combine some features to get newer and cooler features
- # Overall Score of the house (and simplified)
- features['OverallScore'] = features['OverallQual'] * features['OverallCond']
- features['SimplifiedOverallScore'] = features['OverallQual'] * features['SimplifiedOverallCond']
- # Overall Score of the garage (and simplified garage)
- features['GarageScore'] = features['GarageQual'] * features['GarageCond']
- features['SimplifiedGarageScore'] = features['SimplifiedGarageQual'] * features['SimplifiedGarageCond']
- # Overall Score of the exterior (and simplified exterior)
- features['ExterScore'] = features['ExterQual'] * features['ExterCond']
- features['SimplifiedExterScore'] = features['SimplifiedExterQual'] * features['SimplifiedExterCond']
- # Overall Score of the pool (and simplified pool)
- features['PoolScore'] = features['PoolQC'] * features['PoolArea']
- features['SimplifiedPoolScore'] = features['SimplifiedPoolQC'] * features['PoolArea']
- # Overall Score of the kitchens (and simplified kitchens)
- features['KitchenScore'] = features['KitchenQual'] * features['KitchenAbvGr']
- features['SimplifiedKitchenScore'] = features['SimplifiedKitchenQual'] * features['KitchenAbvGr']
- # Overall Score of the fireplaces (and simplified fireplaces)
- features['FireplaceScore'] = features['FireplaceQu'] * features['Fireplaces']
- features['SimplifiedFireplaceScore'] = features['SimplifiedFireplaceQu'] * features['Fireplaces']
- # Total number of bathrooms
- features['TotalBaths'] = features['FullBath'] + (0.5*features['HalfBath']) + features['BsmtFullBath'] + (0.5*features['BsmtHalfBath'])
- ###################################################################################################
- # Box-cox transformation to most skewed features
- numeric_feats = features.dtypes[features.dtypes != 'object'].index
- # Check the skew of all numerical features
- skewed_feats = features[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
- print("\nSkew in numerical features:")
- skewness = pd.DataFrame({'Skew' :skewed_feats})
- skewness.head(10)
- # Box-cox
- skewness = skewness[abs(skewness) > 0.75]
- print("There are {} skewed numerical features to Box Cox transform\n".format(skewness.shape[0]))
- from scipy.special import boxcox1p
- skewed_features = skewness.index
- lam = 0.15
- for feat in skewed_features:
- features[feat] = boxcox1p(features[feat], lam)
- # Label encoding to some categorical features
- categorical_features = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond',
- 'ExterQual', 'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1',
- 'BsmtFinType2', 'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope',
- 'LotShape', 'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond',
- 'YrSold', 'MoSold')
- lbl = LabelEncoder()
- for col in categorical_features:
- lbl.fit(list(features[col].values))
- features[col] = lbl.transform(list(features[col].values))
- # Getting Dummies
- features = pd.get_dummies(features)
- # Splitting features
- train_features = features.loc['train'].select_dtypes(include=[np.number]).values
- test_features = features.loc['test'].select_dtypes(include=[np.number]).values
- # Validation function
- n_folds = 5
- def rmsle_cv(model):
- kf = KFold(n_folds, shuffle=True, random_state=101010).get_n_splits(train_features)
- rmse= np.sqrt(-cross_val_score(model, train_features, train_labels, scoring="neg_mean_squared_error", cv = kf))
- return(rmse)
- # Modelling
- enet_model = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0002, l1_ratio=.9, random_state=101010))
- print(enet_model)
- #score = rmsle_cv(enet_model)
- #print("\nRMSLE: {:.4f} (+/- {:.4f})\n".format(score.mean(), score.std()))
- gb_model = ensemble.GradientBoostingRegressor(n_estimators=3000, learning_rate=0.02,
- max_depth=3, max_features='sqrt',
- min_samples_leaf=15, min_samples_split=10,
- loss='huber', random_state =101010)
- print(gb_model)
- #score = rmsle_cv(gb_model)
- #print("\nRMSLE: {:.4f} (+/- {:.4f})\n".format(score.mean(), score.std()))
- lasso_model = make_pipeline(RobustScaler(), Lasso(alpha=0.0005, random_state=101010))
- print(lasso_model)
- #score = rmsle_cv(lasso_model)
- #print("\nRMSLE: {:.4f} (+/- {:.4f})\n".format(score.mean(), score.std()))
- # Now let's check how do the averaged models work
- averaged_models = AveragingModels(models = (gb_model, enet_model, lasso_model))
- print("AVERAGED MODELS")
- score = rmsle_cv(averaged_models)
- print("\nRMSLE: {:.4f} (+/- {:.4f})\n".format(score.mean(), score.std()))
- # Getting our SalePrice estimation
- averaged_models.fit(train_features, train_labels)
- final_labels = np.exp(averaged_models.predict(test_features))
- # Saving to CSV
- pd.DataFrame({'Id': test_ID, 'SalePrice': final_labels}).to_csv('submission-18.csv', index =False)
|