<jupyter_start><jupyter_text>![](https://storage.googleapis.com/kaggle-datasets-images/1379/2485/a52db2794593657403a4235bfc4147aa/dataset-cover.jpg) Step2 - Regression Ames House Prices: Advanced Regression Techniques(https://www.kaggle.com/c/house-prices-advanced-regression-techniques/overview) Study Regression with 'Boston Housing' Dataset and 'Ames House  Prices' Dataset  Learning Plan 이 과정을 통해서 배울 수 있는 것은 아래와 같습니다.- 첫째, Regression을 이용한 데이터 분석의 전체 Workflow- 둘째, ML - Regression 알고리즘  Lesson 2-1: Overview - Regression Lesson 2-2: EDA Lesson 2-3: Preprocessing Lesson 2-4: Model Selection Lesson 2-5: Model Evaluation & Prediction  Reference - '파이썬 머신러닝 완벽가이드' by 권철민: https://wikibook.co.kr/ml-definitive-guide/- Stacked Regressions : Top 4% on LeaderBoard by [Serigne](https://www.kaggle.com/serigne) : https://www.kaggle.com/serigne/stacked-regressions-top-4-on-leaderboard<jupyter_code>from google.colab import drive
drive.mount('/content/drive')<jupyter_output><empty_output><jupyter_text>Lesson 2-1: Overview - Regression Boston Housing: Predicting Boston Housing Prices(https://www.kaggle.com/samratp/boston-housing-prices-evaluation-validation)<jupyter_code>from sklearn.datasets import load_boston

boston = load_boston()
import pandas as pd
data = pd.DataFrame(boston.data, columns=boston.feature_names)
data['MEDV'] = boston.target
from sklearn.model_selection import train_test_split

train_data, test_data = train_test_split(data, test_size=0.25)
RM_mean = train_data.RM.mean()
RM_std = train_data.RM.std()

MEDV_mean = train_data.MEDV.mean()
MEDV_std = train_data.MEDV.std()

processed_data = (train_data
                  .assign(RM_standard = lambda x: (x.RM - RM_mean)/RM_std)
                  .assign(MEDV_standard = lambda x: (x.MEDV - MEDV_mean)/MEDV_std)
                 )
processed_data.head()
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>2-1-1 : 선형 회귀 (Linear regression) 이 강의에서는 $\theta$를 직접 정하지 않고 머신 러닝 알고리즘을 통해 학습하는 방법을 알아볼 것이다. 이 경우 가장 많이 쓰이는 알고리즘은 선형 회귀(linear regression)이다.구체적인 한 알고리즘에 대해 설명하고 있지만 이 강의에서 설명하는 모델을 정의하고 학습하는 과정은 일반적인 머신 러닝 모델에 모두 적용된다.   Model specification 모델을 정의하는 단계이다.   선형 회귀의 전제설명변수(explanatory variable)와 목적변수(target variable) 사이에 선형 관계가 있다.- 독립변수(independent variable) = 설명변수(explanatory variable) = 입력값이나 원인을 나타낸다. - 종속변수(dependent variable) = 목적변수(target variable) = 결과값이나 효과를 나타낸다. > "All models are wrong; some models are useful." - <NAME>  가정RM(방의 갯수)와 MEDV(집값) 사이에 선형 관계가 있다 이 예시에서는 RM(방의 갯수)와 MEDV(집값) 사이에 선형 관계가 있다고 가정한다. 즉, RM이 올라갈 수록 MEDV도 일정 비율만큼 상승한다고 가정하는 것이다.  Notation * $X$ : input variable, feature* $y$ : output variable, target variable* $\theta$ : model parameters * n : 데이터의 갯수* m : feature의 갯수* $x^{(i)}$ : i번째 데이터* $x^{(i)} = (x^{(i)}_0, x^{(i)}_1, ..., x^{(i)}_m)$* $\theta = (\theta_0, \theta_1, ... )$  가설 (Hypothesis) 가설은 $y$의 예측값을 x의 함수로 표현하는 것이다. 즉, 가설은 다음과 같이 나타낼 수 있다. $$\hat{y} = h_\theta(x)$$ 선형 회귀에서 설명변수 하나와 목적변수 사이의 선형 관계를 가정하고, bias가 없다고 가정하면 다음과 같다. $$h_\theta (x) = \theta_1x$$  가설 (Hypothesis) 이해하기<jupyter_code>def plot_h(theta_1):
    x = np.linspace(-3, 3)
    y_hat = theta_1 * x
    plt.scatter(x=processed_data.RM_standard, y=processed_data.MEDV_standard)
    plt.plot(x, y_hat, 'r');
plot_h(1.2)
plot_h(0.7)<jupyter_output><empty_output><jupyter_text>Cost function 좋은 $\theta$를 찾기 위해서는, 어떤 $\theta$가 좋은 $\theta$인지 평가할 수 있어야 한다.  Cost function은 모델이 얼마나 틀렸는지를 측정하는 함수이다. Loss function, Object function이라고 불리며 $J(\theta)$로 표기한다. $$Cost = J(\theta)$$  선형회귀 문제에서 Cost Function 앞 강의에서 보았듯이 회귀 모델을 평가할 때 널리 쓰이는 방법은 잔차 제곱의 합(RSS : Residual Sum of Squares)이다.  가설 (Hypothesis) $$\hat{y}^{(i)} = h_\theta(x^{(i)}) = \theta_1x^{(i)}$$  매개변수 (Parameter)  $$\theta_1$$  편차 $$\hat{y}^{(i)} - y^{(i)}$$  Cost function (MSE) $$J(\theta) = \frac{1}{2n}\sum^n_{i=1} (h_\theta(x^{(i)}) - y^{(i)})^2$$ 예컨대, $\theta_1$이 1.2와 0.7일 때 RSS는 다음과 같다.<jupyter_code>def RSS(theta_1):
    y_hat = theta_1 * processed_data.RM_standard
    residual = processed_data.MEDV_standard - y_hat
    n = len(residual)
    rss = 1/(2*n) * (residual ** 2).sum()
    return rss
RSS(1.2)
RSS(0.7)<jupyter_output><empty_output><jupyter_text>즉, $\theta_1 = 0.7$일 때 $\theta_1 = 1.2$일 때 보다 조금 더 낫다는 것을 알 수 있다.  $h_\theta(x)$와 $J(\theta)$의 차이$h_\theta(x)$의 함수값은 $\theta$가 고정되어 있고 $x$ 따라 변한다. 이에 비해 $J(\theta)$는 $x$와 $y$가 주어진 데이터로 고정되어 있고 $\theta$에 따라 값이 변한다.   $J(\theta)$의 그래프 $\theta$에 따라 변하는 $J(\theta)$의 그래프를 그려보면 다음과 같다.<jupyter_code>thetas = np.linspace(-3, 3)
plt.plot(thetas, [RSS(theta) for theta in thetas]);<jupyter_output><empty_output><jupyter_text>0.7 ~ 0.8 사이에서 최소값이 있다는 것을 직관적으로 알 수 있다.  Training  경사 하강법 (Gradient Descent) $J(\theta)$를 극소화하는 $\theta$를 자동으로 찾아주는 알고리즘이 경사 하강법 (Gradient Descent)이다.  $$\theta_j := \theta_j - \alpha{\partial\over\partial\theta_j}J(\theta)$$* $\alpha$ : 학습 속도 (learning rate)* ${\partial\over\partial\theta_j}J(\theta)$ : 기울기 (derivative)  $\alpha$의 의미 $\alpha$는 학습 속도 (learning rate)라고 불리며 경사 하강법의 하강 속도를 조절하는 역할을 한다. 이 값은 분석가가 직접 지정해주는 hyperparameter이다.* $\alpha$가 너무 클 때 : 경사 하강법이 수렴하지 못하고 발산하거나 진동한다.* $\alpha$가 너무 작을 때 : 경사 하강법이 수렴하는 데 너무 긴 시간이 걸린다.  극소점 (Local minimum)에 빠질 위험경사 하강법은 최솟값을 찾는 절대적인 알고리즘은 아니다. 최솟값이 아닌 극소점 (local minimum), 즉 주변보다는 낮지만 가장 낮은 값은 아닌 지점에 빠질 위험이 있다.![local minimum](https://files.slack.com/files-pri/T25783BPY-F5YFG1S8N/screenshot.png?pub_secret=7d691ade3a)  In code<jupyter_code>x = processed_data.RM_standard
y = processed_data.MEDV_standard
n = len(x)
alpha = 0.3
theta_1 = initial_theta_1 = -0.2
plot_h(theta_1)
costs = []
iterations = []
for i in range(20):
    y_hat = theta_1 * x
    residual = y_hat - y

    cost = (residual ** 2).sum()/(2*n)

    gradient = 1/n * (residual*x).sum()

    theta_1 -= alpha * gradient
    
    print('theta_1 :', theta_1, 'Cost :', cost, 'gradient :', gradient)
    plot_h(theta_1)
    plt.show()
    
    iterations.append(i)
    costs.append(cost)
plt.plot(iterations, costs)<jupyter_output><empty_output><jupyter_text>2-1-2 : Regularization학습이 과대적합(Overfitting) 되는 것을 방지하기 위해 Cost Function에 **penalty**를 부여하는 방법**L1 Regularization*** 가중치의 제곱의 합이 아닌 **가중치의 합**을 더한 값에 λ를 곱한다. * 일부 가중치(Weight)는 0이 되기도 하는데... 가중치가 0이 된다는 것은... 모델 영향력으로 부터 완전히 제외되는 Feature가 된다는 뜻이다.**L2 Regularization*** 각 **가중치(Weight)의 제곱의 합**에 λ(Regularization Strength)를 곱한다. * λ를 크게 하면 가중치가 더 많이 감소되고, λ를 작게 하면 가중치가 증가한다.**L2 Regularization가 L1 Regularization에 비해 더 안정적이어서 일반적으로 L2 Regularization가 더 많이 사용된다.** **라쏘(Lasso) - L1 Regularization**  $Cost(=Error) = MSE + λ|w|$**릿지(Ridge) - L2 Regularization**$Cost(=Error) = MSE + λw^2$   Lesson 2-2: EDA  2-2-1 : 문제 정의하기  아이오와 주 에임스(Ames)에 있는 주거용 주택을 설명하는 79가지 변수로 각 주택의 최종 가격을 예측합니다.  2-2-2 : 훈련, 테스트 데이터 확인하기<jupyter_code>#dirname = '/content/drive/My Drive/da-structured-data//data/ames_house_prices'
dirname = './data/ames_house_prices'
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
%matplotlib inline
import matplotlib.pyplot as plt  # Matlab-style plotting
import seaborn as sns
color = sns.color_palette()
sns.set_style('darkgrid')
import warnings
def ignore_warn(*args, **kwargs):
    pass
warnings.warn = ignore_warn #ignore annoying warning (from sklearn and seaborn)


from scipy import stats
from scipy.stats import norm, skew #for some statistics


pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x)) #Limiting floats output to 3 decimal points


from subprocess import check_output
print(check_output(["ls", dirname]).decode("utf8")) #check the files available in the directory
#Now let's import and put the train and test datasets in  pandas dataframe

train = pd.read_csv(dirname + '/train.csv')
test = pd.read_csv(dirname + '/test.csv')
##display the first five rows of the train dataset.
train.head(5)
##display the first five rows of the test dataset.
test.head(5)
#explore data
train.info()
#check the numbers of samples and features
print("The train data size before dropping Id feature is : {} ".format(train.shape))
print("The test data size before dropping Id feature is : {} ".format(test.shape))

#Save the 'Id' column
train_ID = train['Id']
test_ID = test['Id']

#Now drop the  'Id' colum since it's unnecessary for  the prediction process.
train.drop("Id", axis = 1, inplace = True)
test.drop("Id", axis = 1, inplace = True)

#check again the data size after dropping the 'Id' variable
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape)) 
print("The test data size after dropping Id feature is : {} ".format(test.shape))
train.head()<jupyter_output><empty_output><jupyter_text>Lesson 2-3: Preprocessing  2-3-1 : Outliers<jupyter_code>fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()<jupyter_output><empty_output><jupyter_text>Training 데이터의 GrLivArea 중 오른쪽 하단에 이상치(outlier) 2개를 볼 수 있습니다. 삭제합니다.Training 데이터의 이상치를 삭제하는 것은 안전합니다. 하지만 Test 데이터의 이상치를 삭제하는 것은 모델에 악영향을 줄 수 있습니다.<jupyter_code>#Deleting outliers
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

#Check the graphic again
fig, ax = plt.subplots()
ax.scatter(train['GrLivArea'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()<jupyter_output><empty_output><jupyter_text>2-3-2 : Target Variable**SalePrice** 는 예측해야할 변수입니다.- [Reference - 정규분포, QQ-Plot](https://datascienceschool.net/view-notebook/76acc92d28354e86940001f9fe85c50f/)<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>**SalePrice** 가 우편향되어 있습니다. 일반적으로 모델은 편향되지 않은 Normal Distribution이 좋습니다. **SalePrice** 분포를 변형합니다.- [Reference - Seaborn](https://datascienceschool.net/view-notebook/4c2d5ff1caab4b21a708cc662137bc65/)- [Reference - Skewed Distributions](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/skewed-distribution/)  Log-transformation of the target variable<jupyter_code>#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()<jupyter_output><empty_output><jupyter_text>2-3-3 : Features Engineering  2-3-3-1 Missing Data 먼저 동일한 DataFrame에 Training 데이터와 Test 데이터를 Concatenate 합니다.<jupyter_code>ntrain = train.shape[0]
ntest = test.shape[0]
y_train = train.SalePrice.values
all_data = pd.concat((train, test)).reset_index(drop=True)
all_data.drop(['SalePrice'], axis=1, inplace=True)
print("all_data size is : {}".format(all_data.shape))<jupyter_output>all_data size is : (2917, 79)<jupyter_text>Missing Data<jupyter_code>all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(20)
f, ax = plt.subplots(figsize=(15, 12))
plt.xticks(rotation='90')
sns.barplot(x=all_data_na.index, y=all_data_na)
plt.xlabel('Features', fontsize=15)
plt.ylabel('Percent of missing values', fontsize=15)
plt.title('Percent missing data by feature', fontsize=15)<jupyter_output><empty_output><jupyter_text>Data Correlation<jupyter_code>#Correlation map to see how features are correlated with SalePrice
corrmat = train.corr()
plt.subplots(figsize=(12,9))
sns.heatmap(corrmat, vmax=0.9, square=True)<jupyter_output><empty_output><jupyter_text>missing values 처리하기  PoolQCdata description says NA means "No Pool". That make sense, given the huge ratio of missing value (+99%) and majority of houses have no Pool at all in general.<jupyter_code>all_data["PoolQC"].value_counts()
all_data["PoolQC"] = all_data["PoolQC"].fillna("None")
all_data["PoolQC"].value_counts()<jupyter_output><empty_output><jupyter_text>MiscFeaturedata description says NA means "no misc feature"<jupyter_code>all_data["MiscFeature"] = all_data["MiscFeature"].fillna("None")<jupyter_output><empty_output><jupyter_text>Alley data description says NA means "no alley access"<jupyter_code>all_data["Alley"] = all_data["Alley"].fillna("None")<jupyter_output><empty_output><jupyter_text>Fencedata description says NA means "no fence"<jupyter_code>all_data["Fence"] = all_data["Fence"].fillna("None")<jupyter_output><empty_output><jupyter_text>FireplaceQudata description says NA means "no fireplace"<jupyter_code>all_data["FireplaceQu"] = all_data["FireplaceQu"].fillna("None")<jupyter_output><empty_output><jupyter_text>LotFrontageSince the area of each street connected to the house property most likely have a similar area to other houses in its neighborhood , we can fill in missing values by the median LotFrontage of the neighborhood.<jupyter_code>all_data["LotFrontage"].describe()
all_data["LotFrontage"].fillna(all_data.groupby("Neighborhood")["LotFrontage"].transform("median"), inplace=True)<jupyter_output><empty_output><jupyter_text>GarageType, GarageFinish, GarageQual and GarageCondReplacing missing data with None<jupyter_code>all_data[['GarageYrBlt','GarageType','GarageCars','GarageArea','GarageQual','GarageCond','GarageFinish']].head() 
for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
    all_data[col] = all_data[col].fillna('None')<jupyter_output><empty_output><jupyter_text>GarageYrBlt, GarageArea and GarageCarsReplacing missing data with 0 (Since No garage = no cars in such garage.)<jupyter_code>for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    all_data[col] = all_data[col].fillna(0)<jupyter_output><empty_output><jupyter_text>BsmtFinSF1, BsmtFinSF2, BsmtUnfSF, TotalBsmtSF, BsmtFullBath and BsmtHalfBathmissing values are likely zero for having no basement<jupyter_code>all_data[['BsmtFinType2','BsmtExposure','BsmtFinType1','BsmtCond','BsmtQual','BsmtFullBath','BsmtHalfBath','TotalBsmtSF','BsmtFinSF1','BsmtFinSF2','BsmtUnfSF']].head() 
for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    all_data[col] = all_data[col].fillna(0)<jupyter_output><empty_output><jupyter_text>BsmtQual, BsmtCond, BsmtExposure, BsmtFinType1 and BsmtFinType2For all these categorical basement-related features, NaN means that there is no basement.<jupyter_code>for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    all_data[col] = all_data[col].fillna('None')<jupyter_output><empty_output><jupyter_text>MasVnrArea and MasVnrTypeNA most likely means no masonry veneer for these houses. We can fill 0 for the area and None for the type.<jupyter_code>all_data["MasVnrType"] = all_data["MasVnrType"].fillna("None")
all_data["MasVnrArea"] = all_data["MasVnrArea"].fillna(0)<jupyter_output><empty_output><jupyter_text>MSZoning (The general zoning classification) 'RL' is by far the most common value. So we can fill in missing values with 'RL'<jupyter_code>all_data['MSZoning'].value_counts()
all_data['MSZoning'] = all_data['MSZoning'].fillna("RL")<jupyter_output><empty_output><jupyter_text>UtilitiesFor this categorical feature all records are "AllPub", except for one "NoSeWa" and 2 NA . Since the house with 'NoSewa' is in the training set, this feature won't help in predictive modelling. We can then safely remove it.<jupyter_code>all_data['Utilities'].value_counts()
all_data = all_data.drop(['Utilities'], axis=1)<jupyter_output><empty_output><jupyter_text>Functionaldata description says NA means typical<jupyter_code>all_data["Functional"] = all_data["Functional"].fillna("Typ")<jupyter_output><empty_output><jupyter_text>ElectricalIt has one NA value. Since this feature has mostly 'SBrkr', we can set that for the missing value.<jupyter_code>all_data['Electrical'].value_counts()
all_data['Electrical'] = all_data['Electrical'].fillna(all_data['Electrical'].mode()[0])<jupyter_output><empty_output><jupyter_text>KitchenQualOnly one NA value, and same as Electrical, we set 'TA' (which is the most frequent) for the missing value in KitchenQual.<jupyter_code>all_data['KitchenQual'].value_counts()
all_data['KitchenQual'].mode()[0]
all_data['KitchenQual'] = all_data['KitchenQual'].fillna(all_data['KitchenQual'].mode()[0])<jupyter_output><empty_output><jupyter_text>Exterior1st and Exterior2ndAgain Both Exterior 1 & 2 have only one missing value. We will just substitute in the most common string<jupyter_code>all_data['Exterior1st'] = all_data['Exterior1st'].fillna(all_data['Exterior1st'].mode()[0])
all_data['Exterior2nd'] = all_data['Exterior2nd'].fillna(all_data['Exterior2nd'].mode()[0])<jupyter_output><empty_output><jupyter_text>SaleTypeFill in again with most frequent which is "WD"<jupyter_code>all_data['SaleType'] = all_data['SaleType'].fillna(all_data['SaleType'].mode()[0])<jupyter_output><empty_output><jupyter_text>MSSubClassNa most likely means No building class. We can replace missing values with None<jupyter_code>all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")<jupyter_output><empty_output><jupyter_text>Check remaining missing values if any<jupyter_code>#Check remaining missing values if any 
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head()
all_data.head(10)
all_data.describe()<jupyter_output><empty_output><jupyter_text>2-3-3-2 Transforming some numerical variables to categoricalRemember...categorical data has no order.<jupyter_code>#MSSubClass=The building class
all_data['MSSubClass'] = all_data['MSSubClass'].apply(str)


#Changing OverallCond into a categorical variable
all_data['OverallCond'] = all_data['OverallCond'].astype(str)


#Year and month sold are transformed into categorical features.
all_data['YrSold'] = all_data['YrSold'].astype(str)
all_data['MoSold'] = all_data['MoSold'].astype(str)<jupyter_output><empty_output><jupyter_text>2-3-3-3 Transforming some categorical variables to numericalLabel Encoding some categorical variables that may contain information in their ordering setThe danger in label encoding is that your machine learning algorithm may learn to favor a over b.The common solution for encoding nominal data is one-hot encoding.(dummy val)<jupyter_code>from sklearn.preprocessing import LabelEncoder
cols = ('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')
# process columns, apply LabelEncoder to categorical features
for c in cols:
    lbl = LabelEncoder() 
    lbl.fit(list(all_data[c].values)) 
    all_data[c] = lbl.transform(list(all_data[c].values))

# shape        
print('Shape all_data: {}'.format(all_data.shape))<jupyter_output>Shape all_data: (2917, 78)<jupyter_text>2-3-3-4 Normalizeing skewd features  Check skewed features<jupyter_code>numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index
numeric_feats
# Check the skew of all numerical features
skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
print("\nSkew in numerical features: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})
skewness.head(10)<jupyter_output>Skew in numerical features:<jupyter_text>Box-Cox Transformation of (highly) skewed featuresWe use the scipy function boxcox1p which computes the Box-Cox transformation of 1+xNote that setting  λ=0 is equivalent to log1p used above for the target variable.See this [page](http://onlinestatbook.com/2/transformations/box-cox.html) for more details on Box Cox Transformation as well as the [scipy function's page](https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.special.boxcox1p.html)[Reference - Box Cox Transformation](https://datascienceschool.net/view-notebook/3f485c426a4b49fc9de95a02137ca6b4/)<jupyter_code>skewness = skewness[abs(skewness) > 0.75]
print("There are {} skewed numerical features to Box Cox transform".format(skewness.shape[0]))

from scipy.special import boxcox1p
skewed_features = skewness.index
lam = 0.15
for feat in skewed_features:
    #all_data[feat] += 1
    all_data[feat] = boxcox1p(all_data[feat], lam)
    
#all_data[skewed_features] = np.log1p(all_data[skewed_features])<jupyter_output>There are 58 skewed numerical features to Box Cox transform<jupyter_text>2-3-3-5 Scaling featuresif need scaling, we can also use minmax scaler[Reference - Scikit-Learn Preprocessing fucntions](https://datascienceschool.net/view-notebook/f43be7d6515b48c0beb909826993c856/)<jupyter_code>#from sklearn.preprocessing import minmax_scale

#df['col_name']= minmax_scale(df['col_name'], axis=0, copy=True)<jupyter_output><empty_output><jupyter_text>we don't need scaling because of the we did log<jupyter_code>all_data[numeric_feats].std()<jupyter_output><empty_output><jupyter_text>2-3-3-6 Add/Drop more features  Check Correlation<jupyter_code>corr = train.corr(method='pearson').drop(['SalePrice']).sort_values('SalePrice', ascending=False)['SalePrice']
corr<jupyter_output><empty_output><jupyter_text>we better delete correlation range from -0.1 ~ 0.1<jupyter_code>all_data.drop('PoolArea', axis=1, inplace=True)
all_data.drop('MoSold', axis=1, inplace=True)
all_data.drop('3SsnPorch', axis=1, inplace=True)
all_data.drop('BsmtFinSF2', axis=1, inplace=True)
all_data.drop('BsmtHalfBath', axis=1, inplace=True)
all_data.drop('MiscVal', axis=1, inplace=True)
all_data.drop('LowQualFinSF', axis=1, inplace=True)
all_data.drop('YrSold', axis=1, inplace=True)
all_data.drop('OverallCond', axis=1, inplace=True)
all_data.drop('MSSubClass', axis=1, inplace=True)
all_data.shape #check deleted 10 cols<jupyter_output><empty_output><jupyter_text>Adding one more important featureSince area related features are very important to determine house prices, we add one more feature which is the total area of basement, first and second floor areas of each house<jupyter_code># Adding total sqfootage feature 
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']
print(all_data.shape)<jupyter_output>(2917, 69)<jupyter_text>Getting dummy categorical features<jupyter_code>all_data = pd.get_dummies(all_data)
print(all_data.shape)<jupyter_output>(2917, 210)<jupyter_text>Getting the new train and test sets.<jupyter_code>#Split train/test set
train = all_data[:ntrain]
test = all_data[ntrain:]<jupyter_output><empty_output><jupyter_text>Lesson 2-4: Model Selection  2-4-1 : Import libratries<jupyter_code>from sklearn.linear_model import ElasticNet, Lasso,  BayesianRidge, LassoLarsIC
from sklearn.ensemble import RandomForestRegressor,  GradientBoostingRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler
from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_squared_error
import xgboost as xgb
import lightgbm as lgb<jupyter_output><empty_output><jupyter_text>2-4-2 : Define a cross validation strategyWe use the cross_val_score function of Sklearn. However this function has not a shuffle attribute, we add then one line of code, in order to shuffle the dataset prior to cross-validation**Cross Validation**- 모델을 평가하는 하나의 방법이다.- K-겹 교차검증(K-fold Cross Validation)을 활용한다.**K-겹 교차검증**- K-겹 교차 검증은 모든 데이터가 최소 한 번은 테스트셋으로 쓰이도록 한다.![](https://static.packt-cdn.com/products/9781789617740/graphics/b04c27c5-7e3f-428a-9aa6-bb3ebcd3584c.png)<jupyter_code>#Validation function
n_folds = 5

def rmsle_cv(model):
    kf = KFold(n_folds, shuffle=True, random_state=42).get_n_splits(train.values)
    rmse= np.sqrt(-cross_val_score(model, train.values, y_train, scoring="neg_mean_squared_error", cv = kf))
    return(rmse)<jupyter_output><empty_output><jupyter_text>2-4-3 : Making Base Regression Model  2-4-3-1 : LASSO Regression This model may be very sensitive to outliers. So we need to made it more robust on them. For that we use the sklearn's Robustscaler() method on pipeline[Reference - RobustScaler](https://mkjjo.github.io/python/2019/01/10/scaler.html)<jupyter_code>lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))<jupyter_output><empty_output><jupyter_text>2-4-3-2 : Elastic Net Regressionagain made robust to outliers**l1_ratio (default=0.5)** * l1_ratio = 0 (L2 Regularization만 사용).  * l1_ratio = 1 (L1 Regularization만 사용).  * 0 < l1_ratio < 1 (L1 Regularization 과 L2 Regularization의 혼합사용)<jupyter_code>ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=3))<jupyter_output><empty_output><jupyter_text>2-4-3-3 : Kernel Ridge Regression[**Polynomial Features**]((https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html?highlight=poly%20featuresklearn.preprocessing.PolynomialFeatures))다항식의 계수를 사용해 **새로운 feature를 생성**한다.예를들면, [a, b] 2개의 feature가 존재한다고 가정하고,degree=2로 설정한다면, polynomial features 는 [1, a, b, a^2, ab, b^2] 가 된다.**Polynomial Features** 를 입력으로 학습하면 일반적으로 더 성능이 좋은 모델이 만들어 진다.<jupyter_code>KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)<jupyter_output><empty_output><jupyter_text>2-4-3-4 : Gradient Boosting RegressionWith huber loss that makes it robust to outliers<jupyter_code>GBoost = GradientBoostingRegressor(n_estimators=3000, learning_rate=0.05,
                                   max_depth=4, max_features='sqrt',
                                   min_samples_leaf=15, min_samples_split=10, 
                                   loss='huber', random_state =5)<jupyter_output><empty_output><jupyter_text>2-4-3-5 : XGBoost<jupyter_code>model_xgb = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, 
                             learning_rate=0.05, max_depth=3, 
                             min_child_weight=1.7817, n_estimators=2200,
                             reg_alpha=0.4640, reg_lambda=0.8571,
                             subsample=0.5213, silent=1,
                             random_state =7, nthread = -1)<jupyter_output><empty_output><jupyter_text>2-4-3-6 : LightGBM<jupyter_code>model_lgb = lgb.LGBMRegressor(objective='regression',num_leaves=5,
                              learning_rate=0.05, n_estimators=720,
                              max_bin = 55, bagging_fraction = 0.8,
                              bagging_freq = 5, feature_fraction = 0.2319,
                              feature_fraction_seed=9, bagging_seed=9,
                              min_data_in_leaf =6, min_sum_hessian_in_leaf = 11)<jupyter_output><empty_output><jupyter_text>Lesson 2-5: Model Evaluation & Prediction  2-5-1 : Base model EvaluationLet's see how these base models perform on the data by evaluating the cross-validation **rmsle** error  2-5-1-1 : LASSO Regression<jupyter_code>score = rmsle_cv(lasso)
print("\nLasso score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))<jupyter_output>Lasso score: 0.1172 (0.0076)<jupyter_text>2-5-1-2 : Elastic Net Regression<jupyter_code>score = rmsle_cv(ENet)
print("ElasticNet score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))<jupyter_output>ElasticNet score: 0.1172 (0.0077)<jupyter_text>2-5-1-3 : Kernel Ridge Regression<jupyter_code>score = rmsle_cv(KRR)
print("Kernel Ridge score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))<jupyter_output>Kernel Ridge score: 0.1211 (0.0078)<jupyter_text>2-5-1-4 : Gradient Boosting Regression<jupyter_code>score = rmsle_cv(GBoost)
print("Gradient Boosting score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))<jupyter_output>Gradient Boosting score: 0.1224 (0.0076)<jupyter_text>2-5-1-5 : XGBoost<jupyter_code>score = rmsle_cv(model_xgb)
print("Xgboost score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))<jupyter_output>Xgboost score: 0.1207 (0.0061)<jupyter_text>2-5-1-6 : LightGBM<jupyter_code>score = rmsle_cv(model_lgb)
print("LGBM score: {:.4f} ({:.4f})\n" .format(score.mean(), score.std()))<jupyter_output>LGBM score: 0.1203 (0.0055)<jupyter_text>2-5-2 : Ensembling앙상블(Ensemble) 학습은 여러 개의 학습 알고리즘을 사용하고, 그 예측을 결합함으로써 보다 정확한 최종 예측을 도출하는 기법입니다. 하나의 강한 머신러닝 알고리즘보다 여러 개의 약한 머신러닝 알고리즘이 낫다 라는 아이디어를 가지고 이해하면 좋습니다.이미지, 영상, 음성 등의 비정형 데이터의 분류는 딥러닝이 뛰어난 성능을 보이고 있지만, 대부분의 정형 데이터 분류 시에는 앙상블이 뛰어난 성능을 나타내고 있습니다.문제와 데이터에 따라 단일 모델의 성능이 더 좋은 경우도 있습니다. 하지만 앙상블 기법을 사용하면 더 유연성있는 모델을 만들며 더 좋은 예측 결과를 기대할 수 있습니다.Regression에서 앙상블 학습의 유형은 가장 많이 알려진 Averaging, Stacking 등이 있으며, 그 외에도 다양한 앙상블 학습의 유형이 연구되고 있습니다.1) Averaging2) Stacking  2-5-2-1 : Averaging에버리징(Averaging)은 가장 쉬운 앙상블 기법입니다. 서로 다른 알고리즘을 가진 Regression 모델를 결합하는 방식입니다.1. Simple Averaging    회귀 문제에서 사용하는 방법으로, 각 예측값을 평균내어 사용합니다. 이 방법은 경우에 따라 과대적합을 줄여주고, 더 부드러운 회귀모델을 만들어줍니다.    2. Weighted Averaging    위에서 평균을 낼 때, 각 모델별 가중치를 두어 평균내는 방식입니다.  Averaged base models class<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Averaged base models scoreWe just average four models here ENet, GBoost, KRR and lasso. Of course we could easily add more models in the mix.<jupyter_code>averaged_models = AveragingModels(models = (ENet, GBoost, KRR, lasso))

score = rmsle_cv(averaged_models)
print(" Averaged base models score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))<jupyter_output>Averaged base models score: 0.1142 (0.0074)<jupyter_text>It seems even the averaging approach really improve the score.  2-5-2-2 : Stacking In this approach, we add a meta-model on averaged base models and use the out-of-folds predictions of these base models to train our meta-model.The procedure, for the training part, may be described as follows:Split the total training set into two disjoint sets (here train and .holdout )Train several base models on the first part (train)Test these base models on the second part (holdout)Use the predictions from 3) (called out-of-folds predictions) as the inputs, and the correct responses (target variable) as the outputs to train a higher level learner called meta-model.The first three steps are done iteratively . If we take for example a 5-fold stacking , we first split the training data into 5 folds. Then we will do 5 iterations. In each iteration, we train every base model on 4 folds and predict on the remaining fold (holdout fold).So, we will be sure, after 5 iterations , that the entire data is used to get out-of-folds predictions that we will then use as new feature to train our meta-model in the step 4.For the prediction part , We average the predictions of all base models on the test data and used them as meta-features on which, the final prediction is done with the meta-model. ![Stacking](http://i.imgur.com/QBuDOjs.jpg)from [Faron](https://www.kaggle.com/getting-started/18153post103381)other reference : [KazAnova's interview](http://blog.kaggle.com/2017/06/15/stacking-made-easy-an-introduction-to-stacknet-by-competitions-grandmaster-marios-michailidis-kazanova/)On this gif, the base models are algorithms 0, 1, 2 and the meta-model is algorithm 3. The entire training dataset is A+B (target variable y known) that we can split into train part (A) and holdout part (B). And the test dataset is C.B1 (which is the prediction from the holdout part) is the new feature used to train the meta-model 3 and C1 (which is the prediction from the test dataset) is the meta-feature on which the final prediction is done.  Stacking averaged Models Class<jupyter_code>class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):
    def __init__(self, base_models, meta_model, n_folds=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.n_folds = n_folds
   
    # We again fit the data on clones of the original models
    def fit(self, X, y):
        self.base_models_ = [list() for x in self.base_models]
        self.meta_model_ = clone(self.meta_model)
        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)
        
        # Train cloned base models then create out-of-fold predictions
        # that are needed to train the cloned meta-model
        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
        for i, model in enumerate(self.base_models):
            for train_index, holdout_index in kfold.split(X, y):
                instance = clone(model)
                self.base_models_[i].append(instance)
                instance.fit(X[train_index], y[train_index])
                y_pred = instance.predict(X[holdout_index])
                out_of_fold_predictions[holdout_index, i] = y_pred
                
        # Now train the cloned  meta-model using the out-of-fold predictions as new feature
        self.meta_model_.fit(out_of_fold_predictions, y)
        return self
   
    #Do the predictions of all base models on the test data and use the averaged predictions as 
    #meta-features for the final prediction which is done by the meta-model
    def predict(self, X):
        meta_features = np.column_stack([
            np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)
            for base_models in self.base_models_ ])
        return self.meta_model_.predict(meta_features)<jupyter_output><empty_output><jupyter_text>Stacking Averaged models ScoreTo make the two approaches comparable (by using the same number of models) , we just average Enet KRR and Gboost, then we add lasso as meta-model.<jupyter_code>stacked_averaged_models = StackingAveragedModels(base_models = (ENet, GBoost, KRR),
                                                 meta_model = lasso)

score = rmsle_cv(stacked_averaged_models)
print("Stacking Averaged models score: {:.4f} ({:.4f})".format(score.mean(), score.std()))<jupyter_output>Stacking Averaged models score: 0.1141 (0.0072)<jupyter_text>We get again a better score by stacking approach.  2-5-2-3 : Ensembling StackedRegressor, XGBoost and LightGBMWe add XGBoost and LightGBM to the StackedRegressor defined previously.  Define rmsle evalution fuction<jupyter_code>def rmse(y, y_pred):
    return np.sqrt(mean_squared_error(y, y_pred))<jupyter_output><empty_output><jupyter_text>StackedRegressor Training<jupyter_code>stacked_averaged_models.fit(train.values, y_train)
stacked_train_pred = stacked_averaged_models.predict(train.values)
stacked_pred = np.expm1(stacked_averaged_models.predict(test.values))
print(rmse(y_train, stacked_train_pred))<jupyter_output>0.08280455767256527<jupyter_text>XGBoost Training<jupyter_code>model_xgb.fit(train, y_train)
xgb_train_pred = model_xgb.predict(train)
xgb_pred = np.expm1(model_xgb.predict(test))
print(rmse(y_train, xgb_train_pred))<jupyter_output>0.08136603761705513<jupyter_text>LightGBM Training<jupyter_code>model_lgb.fit(train, y_train)
lgb_train_pred = model_lgb.predict(train)
lgb_pred = np.expm1(model_lgb.predict(test.values))
print(rmse(y_train, lgb_train_pred))<jupyter_output>0.07712779216682149<jupyter_text>RMSLE on the entire Train data when averaging<jupyter_code>print('RMSLE score on train data:')
print(rmse(y_train,stacked_train_pred*0.70 +
               xgb_train_pred*0.15 + lgb_train_pred*0.15 ))<jupyter_output>RMSLE score on train data:
0.07970002244907942<jupyter_text>2-5-3 : Prediction  Ensemble prediction**Weighted Blending** : 각 모델의 예측값에 대하여 weight를 곱하여 최종 output 계산- 모델에 대한 가중치를 조절하여, 최종 output을 산출- 주의!!! **가중치의 합은 1.0**<jupyter_code>ensemble = stacked_pred*0.70 + xgb_pred*0.15 + lgb_pred*0.15<jupyter_output><empty_output><jupyter_text>Making Prediction Result into CSV File<jupyter_code>result = pd.DataFrame()
result['Id'] = test_ID
result['SalePrice'] = ensemble
result.to_csv(dirname + '/result.csv',index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12247407249011012, 'fraction_numerical': 0.02467122848283973, 'mean_word_length': 5.326175177544809, 'pattern_counts': {'":': 0, '<': 267, '<?xml version=': 0, '>': 267, 'https://': 17, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '70', 'dir': 'jupyter-structured-clean-dedup', 'id': '371673', 'license': 'Apache-2.0', 'n_tokens_mistral': 15621, 'n_tokens_neox': 15770, 'n_words': 3292, 'path': 'step2 - Regression.ipynb', 'repo_name': 'ttc-edu-dev-team/da-structured-data', 'size': '37422'}