<jupyter_start><jupyter_text>**Artificial Intelligence - MSc** ET5003 - MACHINE LEARNING APPLICATIONS  Instructor: <NAME> ET5003_Etivity-2<jupyter_code>#@title Current Date
Today = '2021-08-22' #@param {type:"date"}

#@markdown ---
#@markdown ### Enter your details here:
Student_ID = "20157347" #@param {type:"string"}
Student_full_name = "<NAME>" #@param {type:"string"}
#@markdown ---
#@title Notebook information
Notebook_type = 'Example' #@param ["Example", "Lab", "Practice", "Etivity", "Assignment", "Exam"]
Version = 'Draft' #@param ["Draft", "Final"] {type:"raw"}
Submission = False #@param {type:"boolean"}<jupyter_output><empty_output><jupyter_text>INTRODUCTION **Piecewise regression**, extract from [Wikipedia](https://en.wikipedia.org/wiki/Segmented_regression):Segmented regression, also known as piecewise regression or broken-stick regression, is a method in regression analysis in which the independent variable is partitioned into intervals and a separate line segment is fit to each interval. * Segmented regression analysis can also be performed on multivariate data by partitioning the various independent variables. * Segmented regression is useful when the independent variables, clustered into different groups, exhibit different relationships between the variables in these regions. * The boundaries between the segments are breakpoints.* Segmented linear regression is segmented regression whereby the relations in the intervals are obtained by linear regression.  ***The goal is to use advanced Machine Learning methods to predict House price.***  Imports<jupyter_code># Suppressing Warnings:
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
import arviz as az
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import RobustScaler
from sklearn.preprocessing import MinMaxScaler
import seaborn as sb
from sklearn.model_selection import train_test_split
# to plot
import matplotlib.colors
from mpl_toolkits.mplot3d import Axes3D

# to generate classification, regression and clustering datasets
import sklearn.datasets as dt

# to create data frames
from pandas import DataFrame

# to generate data from an existing dataset
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
# Define the seed so that results can be reproduced
seed = 11
rand_state = 11

# Define the color maps for plots
color_map = plt.cm.get_cmap('RdYlBu')
color_map_discrete = matplotlib.colors.LinearSegmentedColormap.from_list("", ["red","cyan","magenta","blue"])<jupyter_output><empty_output><jupyter_text>DATASET Extract from this [paper](https://ieeexplore.ieee.org/document/9300074):* House prices are a significant impression of the economy, and its value ranges are of great concerns for the clients and property dealers. * Housing price escalate every year that eventually reinforced the need of strategy or technique that could predict house prices in future. * There are certain factors that influence house prices including physical conditions, locations, number of bedrooms and others. 1. [Download the dataset](https://github.com/UL-ET5003/ET5003_SEM1_2021-2/tree/main/Week-3). 2. Upload the dataset into your folder.<jupyter_code>from google.colab import drive
drive.mount('/content/drive')
path = '/content/sample_data/'

col_names =['ad_id', 'area', 'bathrooms', 'beds', 
            'ber_classification', 'county', 'description_block','environment', 
            'facility', 'features', 'latitude', 'longitude', 'no_of_units',
            'price', 'property_category', 'property_type', 'surface']

train_file = path+'house_train.csv'
test_file = path+'house_test.csv'
true_file = path+'true_price.csv'

train_data = pd.read_csv(train_file, usecols=[i for i in col_names if i not in ['features', 'description_block']])
test_data = pd.read_csv(test_file, usecols=[i for i in col_names if i not in ['features', 'description_block', 'price']])
true_data = pd.read_csv(true_file)<jupyter_output><empty_output><jupyter_text>CommentsWe removed 'features' and ' description_block' from the files as these were free-text fields with no useful information.<jupyter_code>print(train_data.shape)
print(test_data.shape)
print(true_data.shape)<jupyter_output>(2982, 15)
(500, 14)
(500, 2)<jupyter_text>Exploratory Data Analysis  Analysis of Training Data - dealing with numeric values<jupyter_code>train_data.head(5)
train_data.describe()<jupyter_output><empty_output><jupyter_text>Analysis of resultsThe above table provides a lot of information for data preparation:1. **count** is the count of non-NaN values and tells us if there are missing values. We can see a number of the above features have missing values and need further analysis.2. **min** and **max** show different ranges and as such need to be scaled to the same range. We will do this as part of Data Preparations.<jupyter_code># Check % of missing values in each feature

train_missing = train_data.isnull().sum()*100/len(train_data)
print(train_missing)
# Let's check for uniqueness in the columns as it seems some columns have only one unique value

print("Unique values in 'county':", train_data.county.unique())
print("Unique values in 'environment':", train_data.environment.unique())<jupyter_output>Unique values in 'county': ['Dublin']
Unique values in 'environment': ['prod']<jupyter_text>Both of these columns have only one value in them and as such we will also remove these.  Analysis of resultsThe above provides some insight into features that can be dropped. Features that are almost empty can be dropped and features with significant number of missing values can also be dropped. **List of features to drop**1. "ber_classification"2. "facilities"3. "no_of_units"4. 'surface'Plus the following clumns with only one unique value.5. 'county'6. 'environment'The 'ad_id' is an identifier and, as we will see below in the correlation map, it is not relevant for training to predict a price. It also will be removed from the training set.<jupyter_code># Drop columns where there is significant missing data

train_data.drop(['ad_id','ber_classification', 'facility', 'no_of_units', 'surface', 'county', 'environment'], axis=1, inplace=True)

# List columns again
train_missing = train_data.isnull().sum()*100/len(train_data)
print(train_missing)<jupyter_output>area                 0.000000
bathrooms            1.710262
beds                 1.710262
latitude             0.000000
longitude            0.000000
price                3.018109
property_category    0.000000
property_type        1.710262
dtype: float64<jupyter_text>The columns now left are in good shape except for a few that have a small percentage of missing vales. Because the % missing is small, we will just remove these rows before further processing.<jupyter_code># Remove full row instances where there are missing values

df_train = train_data.dropna()
# Let's look at features that may have nonsense values

# Remove instances where 'bedrooms' or 'bathrooms' is zero

df_train = df_train.drop(df_train[(df_train.bathrooms==0) | (df_train.beds==0)].index)
# Look at values of 'longitude' and 'latitude' that do not make sense, i.e. outliers

outliers = df_train[(df_train.longitude > -6) | (df_train.latitude < 53)]

outliers.head()<jupyter_output><empty_output><jupyter_text>**Comment**The locations with these coordinates are not in the areas indicated.The Clondalkin location has coordinates for Birmingham.The Glenageary location has coordinates for Bristol.We will remove these two ouliers.<jupyter_code>df_train = df_train.drop(df_train[(df_train.longitude > -6) | (df_train.latitude < 53)].index)
# Let's see how many training instances are left after clean up

print("Number of training instances remaining: ", df_train.shape)

train_missing = df_train.isnull().sum()*100/len(df_train)
print(train_missing)<jupyter_output>Number of training instances remaining:  (2826, 8)
area                 0.0
bathrooms            0.0
beds                 0.0
latitude             0.0
longitude            0.0
price                0.0
property_category    0.0
property_type        0.0
dtype: float64<jupyter_text>CommentsFrom above we have removed columns and rows to improve the quality of the data set. We have 1. Reduced the number of columns from 16 to 82. We have reduced the number of records/instances from 2982 to 2826.The reduction in number of data instances is small and is unlikely to affect the outcome. Reducing the number of columns is a good start.<jupyter_code># Let's look at any correlations between numeric values in the df_train data

sb.heatmap(df_train.corr(), annot=True)
# Let's take a closer look at som eof these correlations

sb.pairplot(data=df_train[['beds','bathrooms', 'latitude', 'longitude', 'price']])<jupyter_output><empty_output><jupyter_text>CommentsWe can see a positive correlation between the 'price' and 'bathrooms' and 'beds'. These are clearly important for price.It is interesting to see that the latitude seems to have a slight negative impact on price and the longitude has a slightly positive impact (from Heatmap). This can be attributed to the dynamics of Dublin prices impacted by location, i.e. there is a distinct cultural driver in North/West versus South Dublin prices.We will extact 'price' as our target (y) value for the model training.  Analysis of Training Data - dealing with categorical values<jupyter_code># Let's deal with Categorical data

# Feature map 'area' first
unique_areas = df_train['area'].unique() # create a list of unique areas
unique_areas_list = np.append(unique_areas, "Kilbarrack")
unique_areas_index = []

# Create a dictionary lookup with Area Name: index pairs for lookup
i = [i+1 for i in range(len(unique_areas_list))]
unique_areas_index = dict(zip(unique_areas_list, i))

# Replace area name with the index from the 'unique_areas_index' dictionary 
df_train_1 = df_train.replace({'area': unique_areas_index})

# Do the same with 'property_type' - there are 12 values here so one-hot-encoding may be a good option here - reduces need for scaling
df_train_1 = pd.get_dummies(df_train_1, columns=['property_type'])

# We are going to map 'property_category' to binary values: it only have two values, it could be one-hot-encoded but it doesn't matter if values are 0,1
mapper_prop_category = {"new_development_parent": 0, "sale": 1}
df_train_1['property_category'] = df_train_1['property_category'].replace(mapper_prop_category)

df_train = df_train_1.copy()
df_train.head()<jupyter_output><empty_output><jupyter_text>Feature distribution and scaling<jupyter_code># Let's look at feature distributions and scaling

# First we will look at the distribution of key predictor features

fig = plt.figure(figsize=(12, 16))

ax1 = fig.add_subplot(5, 2, 1)
ax1.set_title("Plot of 'area'")
df_train.boxplot(column=['area'])

ax2 = fig.add_subplot(5, 2, 2)
ax2.set_title("Histogram of 'area'")
df_train['area'].hist(bins=20)

ax3 = fig.add_subplot(5, 2, 3)
ax3.set_title("Plot of 'bathrooms'")
df_train.boxplot(column=['bathrooms'])

ax4 = fig.add_subplot(5, 2, 4)
ax4.set_title("Histogram of 'bathrooms'")
df_train['bathrooms'].hist(bins=20)

ax5 = fig.add_subplot(5, 2, 5)
ax5.set_title("Plot of 'beds'")
df_train.boxplot(column=['beds'])

ax6 = fig.add_subplot(5, 2, 6)
ax6.set_title("Histogram of 'beds'")
df_train['beds'].hist(bins=20)

ax7 = fig.add_subplot(5, 2, 7)
ax7.set_title("Plot of 'latitude'")
df_train.boxplot(column=['latitude'])

ax8 = fig.add_subplot(5, 2, 8)
ax8.set_title("Histogram of 'latitude'")
df_train['latitude'].hist(bins=20)

ax9 = fig.add_subplot(5, 2, 9)
ax9.set_title("Plot of 'longitude'")
df_train.boxplot(column=['longitude'])

ax10 = fig.add_subplot(5, 2, 10)
ax10.set_title("Histogram of 'logitude'")
df_train['longitude'].hist(bins=20)<jupyter_output><empty_output><jupyter_text>Analysis of distributions and outliersThe above plots indicate outliers (boxplots) and distributions (histograms) of features. The shapes determine the approach to scaling and outlier handling as follows:**area** is not Gaussian and can be scaled using a power transform such as Yeo-Johnson to make it more Gaussian-like**bathrooms** and **beds** have relatively significant outliers and should be treated with the Log Transform**latitude** and **longitude** have outliers but can be treated with the Standard Scaler.We will apply these treatments next.<jupyter_code># Apply scalers for distribution and outlier handling

# Scaling 'area' with Box-Cox to rescale from left skewed to more normal distribution (perfect for non-negative values)
from sklearn.preprocessing import PowerTransformer

bc_scaler = PowerTransformer(method='box-cox')
df_train[['area']] = bc_scaler.fit_transform(df_train[['area']])

# Scaling 'bathrooms' and 'beds'
df_train['bathrooms'] = np.log(df_train['bathrooms'])
df_train['beds'] = np.log(df_train['beds'])

# Scaling 'latitude' and 'longitude'

standard_scaler = StandardScaler()

df_train[['latitude', 'longitude']] = standard_scaler.fit_transform(df_train[['latitude', 'longitude']])

df_train.head()<jupyter_output><empty_output><jupyter_text>CommentsThe above data is now processed to ensure all predictor features are scaled them same.<jupyter_code># Apply all preprocessing to the test dataset

# Drop irrelevant columns from test dataset
test_data.drop(['ad_id','ber_classification', 'facility', 'no_of_units', 'surface', 'county', 'environment'], axis=1, inplace=True)

# Drop samples/instances with null values
df_test = test_data.dropna()

# Remove samples/instances with zero values in 'bathrooms' and 'beds'
df_test = df_test.drop(df_test[(df_test.bathrooms==0) | (df_test.beds==0)].index)

# Remove samples/instances where the coordinates are not correct
df_test = df_test.drop(df_test[(df_test.longitude > -6) | (df_test.latitude < 53)].index)

# Replace area name with the index from the 'unique_areas_index' dictionary y
df_test_1 = df_test.replace({'area': unique_areas_index})

# Do the same with 'property_type' - there are 12 values here so one-hot-encoding may be a good option here - reduces need for scaling
df_test_1 = pd.get_dummies(df_test_1, columns=['property_type'])

# We are going to map 'property_category' to binary values: it only have two values, it could be one-hot-encoded but it doesn't matter if values are 0,1
mapper_prop_category = {"new_development_parent": 0, "sale": 1}
df_test_1['property_category'] = df_test_1['property_category'].replace(mapper_prop_category)

df_test = df_test_1.copy()

# Scaling 'area' with Box-Cox to rescale from left skewed to more normal distribution (perfect for non-negative values)
bc_scaler = PowerTransformer(method='box-cox')
df_test[['area']] = bc_scaler.fit_transform(df_test[['area']])

# Scaling 'bathrooms' and 'beds'
df_test['bathrooms'] = np.log(df_test['bathrooms'])
df_test['beds'] = np.log(df_test['beds'])

# Scaling 'latitude' and 'longitude'
standard_scaler = StandardScaler()

df_test[['latitude', 'longitude']] = standard_scaler.fit_transform(df_test[['latitude', 'longitude']])

df_test.head()
# Let's look at any correlations between numeric values in the df_train data

sb.heatmap(df_train.corr())<jupyter_output><empty_output><jupyter_text>CommentsBased in the correlation matrix above, it is clear that the following may have an imapct on the price of a property:1. property_type_detached2. property_type_apartment3. longitude4. latitude5. beds6. bathrooms7. areaWe wil use these to test the baseline pediction of price.  PIECEWISE REGRESSION  Full Model<jupyter_code># select some features columns just for the baseline model
# assume not all of the features are informative or useful
# in this exercise you could try all of them if possible

featrain = ['area', 'bathrooms', 'beds', 'latitude', 'longitude',
            'property_type_apartment', 'property_type_detached', 'price']
featest = ['area','bathrooms', 'beds', 'latitude', 'longitude',
            'property_type_apartment', 'property_type_detached']

df_subset_train = df_train[featrain]
df_subset_test = df_test[featest]
df_cost =  true_data[true_data.index.isin(df_subset_test.index)]
# Let's look at the distribution of property prices from training data

fig = plt.figure(figsize=(10, 6))

ax0 = fig.add_subplot(1, 2, 1)
ax0.set_title("Plot of 'Price'")
df_subset_train.boxplot(column=['price'])

ax00 = fig.add_subplot(1, 2, 2)
ax00.set_title("Histogram of 'Price'")
df_subset_train['price'].hist(bins=20)<jupyter_output><empty_output><jupyter_text>CommentThese data will need to be log and standard scaled.<jupyter_code># Separate training data into target 'y' and predictor 'X' values

# train set, input columns
X_train = df_subset_train.iloc[:,0:-1].values 
# train set, output column, cost
y_train = df_subset_train.iloc[:,-1].values.reshape(-1,1)   

# test set, input columns
X_test = df_subset_test.iloc[:,0:].values 
# test set, output column, cost
y_test = df_cost.Expected.values
# Scale y_train 'price' values

ylog = np.log(y_train.astype('float'))
yscaler = StandardScaler().fit(ylog)
y_train = yscaler.transform(ylog)
# model
with pm.Model() as model:
    #prior over the parameters of linear regression
    alpha = pm.Normal('alpha', mu=0, sigma=30)
    #we have one beta for each column of Xn
    beta = pm.Normal('beta', mu=0, sigma=30, shape=X_train.shape[1])
    #prior over the variance of the noise
    sigma = pm.HalfCauchy('sigma_n', 5)
    #linear regression model in matrix form
    mu = alpha + pm.math.dot(beta, X_train.T)
    #likelihood, be sure that observed is a 1d vector
    like = pm.Normal('like', mu=mu, sigma=sigma, observed=y_train[:,0])
    
#number of iterations of the algorithms
iter = 50000

# run the model
with model:
    approximation = pm.fit(iter,method='advi')
    
# check the convergence
plt.plot(approximation.hist);
# samples from the posterior
posterior = approximation.sample(draws=2000)
# prediction
ll=np.mean(posterior['alpha']) + np.dot(np.mean(posterior['beta'],axis=0), X_test.T)
y_pred_BLR = np.exp(yscaler.inverse_transform(ll.reshape(-1,1)))[:,0]
print("MAE = ",(np.mean(abs(y_pred_BLR - y_test))))
print("MAPE = ",(np.mean(abs(y_pred_BLR - y_test) / y_test)))<jupyter_output>MAE =  203412.09453346493
MAPE =  0.2643648015412629<jupyter_text>Analysis of first modelThe results for MAE and MAPE are very similar for both the reduced feature set above and the full set of features. Adding in or reoving any other feature reduced the performance of the model, sometimes considerably. We will use this measure of performance a our baseline.  Clustering  Full Model<jupyter_code># training gaussian mixture model 
from sklearn.mixture import GaussianMixture

gmm = GaussianMixture(n_components=4)<jupyter_output><empty_output><jupyter_text>Clusters<jupyter_code># train clusters


# test clusters<jupyter_output><empty_output><jupyter_text>Piecewise Model<jupyter_code># model_0
with pm.Model() as model_0:
  # prior over the parameters of linear regression
  alpha = pm.Normal('alpha', mu=0, sigma=30)
  # we have a beta for each column of Xn0
  beta = pm.Normal('beta', mu=0, sigma=30, shape=Xn0.shape[1])
  # prior over the variance of the noise
  sigma = pm.HalfCauchy('sigma_n', 5)
  # linear regression relationship
  #linear regression model in matrix form
  mu = alpha + pm.math.dot(beta, Xn0.T)
  # likelihood, be sure that observed is a 1d vector
  like = pm.Normal('like', mu=mu, sigma=sigma, observed=yn0[:,0])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10365509816988773, 'fraction_numerical': 0.020197877684933613, 'mean_word_length': 5.619613165931455, 'pattern_counts': {'":': 4, '<': 77, '<?xml version=': 0, '>': 77, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '54529', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6383, 'n_tokens_neox': 6016, 'n_words': 1999, 'path': 'ET5003_Etivity2_RobertBarrett_20157347 (5).ipynb', 'repo_name': 'rb2022/ET5003_SEM1_2021-2', 'size': '19523'}