<jupyter_start><jupyter_text>- Homework source https://github.com/alexeygrigorev/mlbookcamp-code/blob/master/course-zoomcamp/06-trees/homework.md- Lecture https://github.com/alexeygrigorev/mlbookcamp-code/blob/master/chapter-06-trees/06-trees.ipynb 6.10 HomeworkThe goal of this homework is to create a tree-based regression model for prediction apartment prices (column `'price'`).In this homework we'll again use the New York City Airbnb Open Data dataset - the same one we used in homework 2 and 3.You can take it from [Kaggle](https://www.kaggle.com/dgomonov/new-york-city-airbnb-open-data?select=AB_NYC_2019.csv)or download from [here](https://raw.githubusercontent.com/alexeygrigorev/datasets/master/AB_NYC_2019.csv)if you don't want to sign up to Kaggle.For this homework, we prepared a [starter notebook](homework-6-starter.ipynb).  Loading the data* Use only the following columns:    * `'neighbourhood_group',`    * `'room_type',`    * `'latitude',`    * `'longitude',`    * `'minimum_nights',`    * `'number_of_reviews','reviews_per_month',`    * `'calculated_host_listings_count',`    * `'availability_365',`    * `'price'`* Fill NAs with 0* Apply the log tranform to `price`* Do train/validation/test split with 60%/20%/20% distribution. * Use the `train_test_split` function and set the `random_state` parameter to 1* Use `DictVectorizer` to turn the dataframe into matrices<jupyter_code>col = ['neighbourhood_group',
    'room_type',
    'latitude',
    'longitude',
    'minimum_nights',
    'number_of_reviews','reviews_per_month',
    'calculated_host_listings_count',
    'availability_365',
    'price']
df = (
 pd.read_csv('../input/new-york-city-airbnb-open-data/AB_NYC_2019.csv')
[col]
.fillna(0)
) 
df['price'] = np.log1p(df['price'])
df.head()
y='price'
test=0.2
val=0.2
seed=1

df_train_full, df_test = train_test_split(df, test_size=test, random_state=seed)
df_train, df_val = train_test_split(df_train_full, test_size=val/(1-test), random_state=seed)

y_test = df_test[y].copy().values
y_val = df_val[y].copy().values
y_train = df_train[y].copy().values
del df_test[y]
del df_val[y]
del df_train[y]
# hot encoding
dict_train = df_train.to_dict(orient='records')
dict_val = df_val.to_dict(orient='records')
dv = DictVectorizer(sparse=False)
X_train = dv.fit_transform(dict_train)
X_val = dv.transform(dict_val)<jupyter_output><empty_output><jupyter_text>Question 1Let's train a decision tree regressor to predict the price variable. * Train a model with `max_depth=1`Which feature is used for splitting the data?* `room_type`* `neighbourhood_group`* `number_of_reviews`* `reviews_per_month`<jupyter_code># https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html
from sklearn.tree import DecisionTreeRegressor
dt = DecisionTreeRegressor(max_depth=1)
dt.fit(X_train, y_train)
from sklearn.tree import export_text
print(export_text(dt, feature_names=dv.get_feature_names()))
from sklearn.tree import plot_tree
plot_tree(dt, feature_names=dv.get_feature_names())
plt.show()
# first node
feature_id = dt.tree_.feature[0] # [12, -2, -2]
feature_name = dv.get_feature_names()[feature_id] # 'room_type=Entire home/apt'
md(f'### Which feature is used for splitting the data?: **{feature_name.split("=")[0]}**')<jupyter_output><empty_output><jupyter_text>Question 2Train a random forest model with these parameters:* `n_estimators=10`* `random_state=1`* `n_jobs=-1`  (optional - to make training faster)What's the RMSE of this model on validation?* 0.059* 0.259* 0.459* 0.659<jupyter_code># https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html?highlight=randomforest#sklearn.ensemble.RandomForestRegressor
from sklearn.ensemble import RandomForestRegressor
def get_rmse(y_pred, y_true):
    mse = ((y_pred - y_true) ** 2).mean()
    return np.sqrt(mse)
#enddef
rf = RandomForestRegressor(n_estimators=10, random_state=1) # n_jobs=-1
rf.fit(X_train, y_train)
md(f"### What's the RMSE of this model on validation? : **{get_rmse(rf.predict(X_val), y_val):.4f}**")
# 0.4599<jupyter_output><empty_output><jupyter_text>Question 3Now let's experiment with the `n_estimators` parameter* Try different values of this parameter from 10 to 200 with step 10* Set `random_state` to `1`* Evaluate the model on the validation datasetAfter which value of `n_estimators` does RMSE stop improving?- 10- 50- 70- 120<jupyter_code>with shelve.open(savefile, 'c') as save: 
    k = 'rmse01'
    if k not in save:
        rmse_list = {}
        for n in np.linspace(10,200,10).astype(int):
            rf = RandomForestRegressor(n_estimators=n, random_state=1,n_jobs=-1) 
            rf.fit(X_train, y_train)
            rmse_list[n] = get_rmse(rf.predict(X_val), y_val)
            print(n,rmse_list[n])
        #endfor
        save[k] = rmse_list
    else:
        rmse_list = save[k]
    #endif
#endwith        
pd.Series(rmse_list).plot()
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>After which value of n_estimators does RMSE stop improving? **120**  Question 4Let's select the best `max_depth`:* Try different values of `max_depth`: `[10, 15, 20, 25]`* For each of these values, try different values of `n_estimators` from 10 till 200 (with step 10)* Fix the random seed: `random_state=1`What's the best `max_depth`:* 10* 15* 20* 25<jupyter_code>rmse_list02 = None
with shelve.open(savefile, 'c') as save: 
    k = 'rmse02'
    if k not in save:
        rmse_list02 = {}
        for d in [10, 15, 20, 25]:
            rmse_list02[d] = rmse_list02.get(d,{}) # create empty Dictionary if key doesn't exist yet
            for n in np.linspace(10,200,10).astype(int):
                if n not in rmse_list02[d]:
                    rf = RandomForestRegressor(n_estimators=n
                                               ,max_depth=d
                                               ,random_state=1
                                              ,n_jobs=-1) # 
                    rf.fit(X_train, y_train)
                    rmse_list02[d][n] = get_rmse(rf.predict(X_val), y_val)
                #endif
                print(d,n,rmse_list02[d][n])
            #endfor
        #endfor
        save[k]= rmse_list02
    else:
        rmse_list02 = save[k]
    #endif
#endwith
plt.figure(figsize=(6, 4))
for d in [10, 15, 20, 25]:
    x = rmse_list02[d].keys()
    y = [rmse_list02[d][n] for n in x]
    plt.plot(x, y, label=f'depth={d}')
#endfor
plt.xticks(range(0, 201, 10))
plt.grid()
plt.legend()
plt.xlabel('n_estimators')
plt.ylabel('rmse')
plt.show()
res = { min([rmse_list02[d][n] for n in rmse_list02[d]]) : d for d in rmse_list02 }
md(f"### What's the best `max_depth`? : **{res[sorted(res)[0]]}**") # 15<jupyter_output><empty_output><jupyter_text>**Bonus question (not graded):**Will the answer be different if we change the seed for the model?**Answer**: it should *not*, since n_estimators is sufficently high (>100).  Question 5We can extract feature importance information from tree-based models. At each step of the decision tree learning algorith, it finds the best split. When doint it, we can calculate "gain" - the reduction in impurity before and after the split. This gain is quite useful in understanding what are the imporatant features for tree-based models.In Scikit-Learn, tree-based models contain this information in the[`feature_importances_`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.htmlsklearn.ensemble.RandomForestRegressor.feature_importances_)field. For this homework question, we'll find the most important feature:* Train the model with these parametes:    * `n_estimators=10`,    * `max_depth=20`,    * `random_state=1`,    * `n_jobs=-1` (optional)* Get the feature importance information from this modelWhat's the most important feature? * `neighbourhood_group=Manhattan`* `room_type=Entire home/apt`	* `longitude`* `latitude`<jupyter_code>rf = RandomForestRegressor(n_estimators=10
                           ,max_depth=20
                           ,random_state=1
                          ,n_jobs=-1) # 
rf.fit(X_train, y_train)
importances = list(zip(dv.feature_names_, rf.feature_importances_))
df_importance = (
    pd.DataFrame(importances, columns=['feature', 'gain'])
#       [lambda x : x['gain'] > 0]
     .sort_values(by='gain', ascending=False)
)
md(f"### What's the most important feature? : **{df_importance['feature'].iloc[0]}**")
# room_type=Entire home/apt<jupyter_output><empty_output><jupyter_text>Question 6Now let's train an XGBoost model! For this question, we'll tune the `eta` parameter* Install XGBoost* Create DMatrix for train and validation* Create a watchlist* Train a model with these parameters for 100 rounds:```xgb_params = {    'eta': 0.3,     'max_depth': 6,    'min_child_weight': 1,        'objective': 'reg:squarederror',    'nthread': 8,        'seed': 1,    'verbosity': 1,}```Now change `eta` first to `0.1` and then to `0.01`What's the best eta?* 0.3* 0.1* 0.01<jupyter_code>import xgboost as xgb # Install XGBoost
def parse_xgb_output(output):
    tree = []
    p_train = []
    p_val = []

    for line in output.stdout.strip().split('\n'):
        it_line, train_line, val_line = line.split('\t')

        it = int(it_line.strip('[]'))
        train = float(train_line.split(':')[1])
        val = float(val_line.split(':')[1])

        tree.append(it)
        p_train.append(train)
        p_val.append(val)

    return tree, p_train, p_val
#enddef
# Create DMatrix for train and validation
dtrain = xgb.DMatrix(X_train, label=y_train, feature_names=dv.feature_names_)
dval = xgb.DMatrix(X_val, label=y_val, feature_names=dv.feature_names_)
# Create a watchlist
watchlist = [(dtrain, 'train'), (dval, 'val')]
# Train a model with these parameters for 100 rounds:
xgb_params = {
    'eta': 0.3, 
    'max_depth': 6,
    'min_child_weight': 1,

    'objective': 'reg:squarederror',
    'nthread': 8,

    'seed': 1,
    'verbosity': 1,
}
%%capture output
# capture instruction that saves the result to output 
xgb_params['eta'] = 0.3
model = xgb.train(xgb_params, dtrain,
                  num_boost_round=100,
                  evals=watchlist, verbose_eval=10)
tree, p_train, p_val = parse_xgb_output(output)
print(f'Eta={xgb_params["eta"]} : Best performance (squarederror, number of trees) ', min(zip(p_val, tree)))

plt.figure(figsize=(6, 4))
plt.plot(tree, p_train, color='black', linestyle='dashed', label='Train Loss')
plt.plot(tree, p_val, color='black', linestyle='solid', label='Validation Loss')
# plt.xticks(range(0, 101, 25))
plt.legend()
plt.title('XGBoost: number of trees vs "squarederror"')
plt.xlabel('Number of trees')
plt.ylabel('squarederror')
plt.yscale('log')
plt.show()
%%capture output_010
# capture instruction that saves the result to output 
xgb_params['eta'] = 0.1
model = xgb.train(xgb_params, dtrain,
                  num_boost_round=100,
                  evals=watchlist, verbose_eval=10)
tree, _, p_val = parse_xgb_output(output_010)
print(f'Eta={xgb_params["eta"]} : Best performance (squarederror, number of trees) ', min(zip(p_val, tree)))
%%capture output_001
# capture instruction that saves the result to output 
xgb_params['eta'] = 0.01
model = xgb.train(xgb_params, dtrain,
                  num_boost_round=100,
                  evals=watchlist, verbose_eval=10)
tree, _, p_val = parse_xgb_output(output_001)
print(f'Eta={xgb_params["eta"]} : Best performance (squarederror, number of trees) ', min(zip(p_val, tree)))
plt.figure(figsize=(6, 4))
for eta, out in zip([0.3,0.1,0.01],[output,output_010,output_001]):
    tree, _, p_val = parse_xgb_output(out)
    #plt.plot(tree, p_train, color='black', linestyle='dashed', label='eta=eta, Train Loss')
    plt.plot(tree, p_val, linestyle='solid', label=f'eta={eta}')
    print(f'Eta={eta} : Best performance (squarederror, number of trees) ', min(zip(p_val, tree)))
    
# plt.xticks(range(0, 101, 25))
plt.legend()
plt.title('XGBoost: number of trees vs Validation "squarederror"')
plt.xlabel('Number of trees')
plt.ylabel('Validation "squarederror"')
# plt.yscale('log')
plt.ylim(0.43,0.46)
plt.grid()
plt.show()
md("### What's the best eta? **0.1**")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13675493973928016, 'fraction_numerical': 0.027383782897433796, 'mean_word_length': 4.195059625212947, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 31, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '9228', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4600, 'n_tokens_neox': 4234, 'n_words': 1018, 'path': '06-ml-zoomcamp-homework-trees.ipynb', 'repo_name': 'alexkolo/202109-ML-Zoomcamp-by-Alexey-Grigorev', 'size': '12197'}