# Install required packages !pip install --upgrade pip !pip install -q -r requirements.txt

# sklearn modules import sklearn from sklearn.linear_model import LinearRegression from sklearn.preprocessing import PolynomialFeatures from sklearn.metrics import mean_squared_error from sklearn.pipeline import make_pipeline from sklearn.compose import ColumnTransformer, make_column_transformer from sklearn.preprocessing import StandardScaler from sklearn.model_selection import GridSearchCV, KFold from sklearn.linear_model import Ridge from sklearn.linear_model import RidgeCV from sklearn.linear_model import Lasso # pandas, numpy, and plotting import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import scipy.stats as st from scipy.stats import norm # schrutepy from schrutepy import schrutepy np.random.seed(44)

# Display plots inline %matplotlib inline # Plotting defaults plt.rcParams['figure.figsize'] = (8,5) plt.rcParams['figure.dpi'] = 80 df = schrutepy.load_schrute()

#Combining number of reviews and average review rating, deleting previous columns d = pd.read_csv("the_office.csv") tv = d["total_votes"].mean() d["Episode Scor"] = 0.7*d["imdb_rating"] + 10*(0.3)*(1-np.exp((d["total_votes"])/np.log(1/2)*tv)) d['Episode Score'] = d['Episode Scor'].apply(lambda x: round(x, 2)) del d['total_votes'] del d['imdb_rating'] del d['Episode Scor']

plt.hist(d['Episode Score'],bins=10,alpha=0.7,color = 'green', density = True) plt.xlabel('Episode Score') plt.ylabel('Frequency Density') plt.title('Histogram of Episode Scores') score_fit = norm.fit(d['Episode Score']) x_vals = np.linspace(7.4,10.4,1000) pdf = norm.pdf(x_vals,score_fit[0],score_fit[1]) plt.plot(x_vals,pdf)

nam = list(d["episode_name"]) cam_1 = ["Hot Girl","Office Olympics","The Fire","The Client","E-Mail Surveilance","Valentine's Day","The Convention","The Return","Launch Party (Parts 1&2)","Business Trip","New Boss","Two Weeks","Dream Team","The Michael Scott Paper Company","Stairmageddon","Promos","Livin' the Dream","Heavy Competition","Broke","Sabre","Mafia","The Delivery (Parts 1&2)","Nepotism","The Sting","The Seminar","Training Day","Michael's Last Dundies","Goodbye, Michael","The Inner Circle","Mrs. California","Turf War","Here Comes Treble","The Boat","Lice","Paper Airplane"] cam_2 = ["Booze Cruise","Garden Party"] cam_3 = ["Stress Relief (Parts 1&2)"] cam_4 = ["Finale"] cam_5 = ["Search Committee (Parts 1&2)"] cameos = [0] * 186 for i in cam_1: cameos[nam.index(i)] = 1 for i in cam_2: cameos[nam.index(i)] = 2 for i in cam_3: cameos[nam.index(i)] = 3 for i in cam_4: cameos[nam.index(i)] = 4 for i in cam_5: cameos[nam.index(i)] = 5 d["Number of Cameos"] = cameos

d["Episode Length (Minutes)"] = (23,23,22,23,23,23,21,22,22,22,30,22,22,22,23,22,22,22,22,22,22,22,22,21,22,22,22,22, 23,22,22,22,22,22,31,30,30,42,21,22,28,21,21,22,21,30,21,21,29,28,42,42,42,42,42,22,30,30,30,30,30,22,30,42,42,43,22, 22,30,22,22,30,30,22,30,30,22,60,60,60,22,30,30,30,30,21,22,30,22,30,30,30,30,30,30,49,30,30,30,23,30,30,30,30,30,30, 60,30,23,30,30,30,30,30,30,30,30,30,30,30,30,23,30,30,30,60,30,30,30,22,30,30,30,30,20,50,23,23,42,23,23,23,23,23,23, 23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,22,22,22,22,22,22,22,22,22,22,22,22,42,22,22,22,22,42, 43,51,) d["Episode Location (Binary)"] = (0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,1,0,0,0,1,0,0,0,0, 1,1,1,0,0,0,1,1,0,1,0,0,0,0,1,1,1,1,0,1,0,1,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1, 0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,1,1,1,1,1,1,1,1,0,0,1, 0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1)

#finding best director unique_directors = d['director'].unique().tolist() directors = list(d["director"]) director_count = [] for i in unique_directors: dir_count = 0 for j in d.director: if j == i: dir_count += 1 director_count.append(dir_count) count = 0 names = [] for i in range(0,len(unique_directors)): if director_count[i] >5: names.append(unique_directors[i]) count+=1 for i in directors: if i not in names: directors[directors.index(i)] = "others" directors1 = pd.get_dummies(directors) directors1.columns = ['D: Charles McDougall', 'D: David Rogers','D: Greg Daniels','D: Jeffrey Blitz','D: Ken Kwapis', 'D: Ken Whittingham' , 'D: Matt Sohn','D: Paul Feig','D: Paul Lieberstein','D: Randall Einhorn','D: Others']

writers = list(d.writer) for j in writers: if j == 'Lee Eisenberg;Gene Stupnitsky': writers[writers.index(j)] = "Gene Stupnitsky;Lee Eisenberg"

#finding best writers unique_writers = d['writer'].unique().tolist() writer_count = [] for i in unique_writers: wri_count = 0 for j in d.writer: if j == i: wri_count += 1 writer_count.append(wri_count) count = 0 names = [] for i in range(0,len(unique_writers)): if writer_count[i] >6: names.append(unique_writers[i]) count+=1 for i in writers: if i not in names: writers[writers.index(i)] = "others" writers1 = pd.get_dummies(writers) writers1.columns = ['W: B.J. Novak', 'W:Brent Forrester','W:Charlie Grandy','W:Gene Stupnitsky;Lee Eisenberg', 'W:Greg Daniels', 'W:Jennifer Celotta' , 'W:Justin Spitzer','W:Michael Schur','W:Mindy Kaling','W:Paul Lieberstein', 'W: Others']

d = pd.concat([d, writers1], axis=1) d = pd.concat([d, directors1], axis=1)

# Lineshare (number of lines spoken by characters) # Remove records that don't contain dialogue df = df.dropna() df episode_names = d.episode_name # Main 39 characters main_chars = ["Michael","Dwight","Jim","Pam","Andy","Angela","Kevin","Erin","Oscar","Ryan"] n_main_chars = len(main_chars) # Number of main characters n_episodes = len(episode_names) # Number of episodes array = np.zeros((n_episodes,n_main_chars)) # Create array of zeros with episode number as rows and character names as columns # Create array of number of lines spoken by each characters for i in episode_names: for j in main_chars: episode_names_list = episode_names.values.tolist() row = episode_names_list.index(i) col = main_chars.index(j) episode = df.loc[df['episode_name'] == i] character_lines = episode.loc[episode['character'] == j] number_lines = len(character_lines) array[row,col] = number_lines # Create array with proportion of lines spoken by each character

array_proportion = np.zeros((186,10)) for i in range(0,186): row_tot = sum(array[i,:]) for j in range(0,10): array_proportion[i,j] = array[i,j]/row_tot # To produce models using these variables without skewing the dependency on variables, we reduce the number of characters we consider ten_main_prop = array_proportion # To show that all columns are indepent, we show ten_main_prop has full rank print("The rank of ten_main_prop is:",np.linalg.matrix_rank(ten_main_prop)) # Add characters to dataframe for i in range(0,10): d[main_chars[i]] = ten_main_prop[:,i]

pd.DataFrame(new_d)

pd.DataFrame(new_d) new_d = d.drop(columns=['season','episode','director','writer','air_date','main_chars','episode_name'])

new_d_main_vars = new_d[['n_lines','n_directions','n_words','n_speak_char','Episode Score','Number of Cameos','Episode Length (Minutes)']] plt.figure(figsize=(12, 8)) mask = np.triu(np.ones_like(new_d_main_vars.corr(), dtype=np.bool)) corr = new_d_main_vars.corr() kot = corr[abs(corr)>=0.5] ax = plt.axes() sns.heatmap(kot, mask=mask, vmin=-1, vmax=1, annot=True,ax=ax) ax.set_title("Correlation Heatmap between Main Variables", fontsize = 16) plt.show()

new_d_characters = new_d[["Michael","Dwight","Jim","Pam","Andy","Angela","Kevin","Erin","Oscar","Ryan"]] plt.figure(figsize=(12,8)) mask = np.triu(np.ones_like(new_d_characters.corr(),dtype=np.bool)) corr1 = new_d_characters.corr() kot1 = corr1[abs(corr1)>=0.35] ax = plt.axes() sns.heatmap(kot1, mask=mask, vmin=-1, vmax=1, annot=True, ax = ax) ax.set_title("Heatmap of correlation between main characters", fontsize = 16) plt.show()

new_d = new_d.drop(columns=["n_lines","n_words","Dwight","Oscar","Erin","Andy","Kevin"])

# List of non binary and none dropped after heat map correlation test features. non_binary_column_names = ['n_directions','n_speak_char','Number of Cameos','Episode Length (Minutes)', 'Michael','Jim','Pam','Angela','Ryan']

fig, axes = plt.subplots(3,3, figsize=(12, 12), sharey=True) degree_range = 5 colours = ['red','green','black','magenta','cyan'] for i in range(3): for j in range(3): x = new_d[non_binary_column_names[3*i+j]] y = new_d['Episode Score'] sns.scatterplot(ax=axes[i,j],x = new_d[non_binary_column_names[3*i+j]],y = new_d['Episode Score']) for r in range(1,degree_range): degree = r X = pd.DataFrame(x) X_seq = np.linspace(X.min(),X.max(),186).reshape(-1,1) polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression()) polyreg.fit(X,y) axes[i,j].plot(X_seq,polyreg.predict(X_seq),color=colours[r-1], label = 'Degree =' + str(r)) if j == 0: if i == 0: axes[i,j].legend()

new_d = new_d.drop(columns=["Jim","Ryan","Pam","Angela"])

#helper functions def get_coefs(m): """Returns the model coefficients from a Scikit-learn model object as an array, includes the intercept if available. """ # If pipeline, use the last step as the model if (isinstance(m, sklearn.pipeline.Pipeline)): m = m.steps[-1][1] if m.intercept_ is None: return m.coef_ return np.concatenate([[m.intercept_], m.coef_]) def model_fit(m, X, y, title, plot = False): """Returns the root mean squared error of a fitted model based on provided X and y values. Args: m: sklearn model object X: model matrix to use for prediction y: outcome vector to use to calculating rmse and residuals plot: boolean value, should fit plots be shown """ y_hat = m.predict(X) rmse = mean_squared_error(y, y_hat, squared=False) res = pd.DataFrame( data = {'y': y, 'y_hat': y_hat, 'resid': y - y_hat} ) if plot: plt.figure(figsize=(12, 6)) plt.subplot(121) sns.lineplot(x='y', y='y_hat', color="grey", data = pd.DataFrame(data={'y': [min(y),max(y)], 'y_hat': [min(y),max(y)]})) sns.scatterplot(x='y', y='y_hat', data=res).set_title("Fit plot") plt.subplot(122) sns.scatterplot(x='y', y='resid', data=res).set_title("Residual plot") plt.subplots_adjust(left=0.0) plt.suptitle(title + str(round(rmse, 4)), fontsize=16) plt.show() return rmse

kf = KFold(n_splits = 5, shuffle = True,random_state=1) result = next(kf.split(new_d), None) train = new_d.iloc[result[0]] test = new_d.iloc[result[1]] x_train = train.drop(columns=['Episode Score']) y_train = train['Episode Score']

# Non scaled lm = LinearRegression(fit_intercept = False).fit(x_train, y_train) #From now on we will only fit to x_train with scaled categorical features scaled_features = x_train col_names = ['n_directions', 'n_speak_char', 'Michael', 'Number of Cameos','Episode Length (Minutes)'] features = scaled_features[col_names] scaler = StandardScaler().fit(features.values) features = scaler.transform(features.values) scaled_features[col_names] = features x_train_scaled = scaled_features lm_scaled = make_pipeline( LinearRegression(fit_intercept=False) ).fit(x_train_scaled, y_train)

# Tuning linear ridge model for best alpha alphas = np.linspace(0.0000001, 0.01, num=200) lm_scaled_ridge = GridSearchCV( make_pipeline( Ridge(fit_intercept = False) ), param_grid={'ridge__alpha': alphas}, cv=KFold(5, shuffle=True), scoring="neg_root_mean_squared_error" ).fit(x_train_scaled, y_train)

get_coefs(lm_scaled_ridge.best_estimator_)[1:30]

alphas = np.linspace(0.0000001, 1, num=100) lm_scaled_lasso = GridSearchCV( make_pipeline( Lasso(fit_intercept = False) ), param_grid={'lasso__alpha': alphas}, cv=KFold(5, shuffle=True), scoring="neg_root_mean_squared_error" ).fit(x_train_scaled, y_train)

get_coefs(lm_scaled_lasso.best_estimator_)[1:30]

poly_model = make_pipeline( make_column_transformer( (PolynomialFeatures(include_bias=False), ['n_directions']), (PolynomialFeatures(include_bias=False), ['n_speak_char']), (PolynomialFeatures(include_bias=False), ['Number of Cameos']), remainder='passthrough'), LinearRegression(fit_intercept=False)) parameters_poly_model = { 'columntransformer__polynomialfeatures-1__degree': np.arange(1,4,1), 'columntransformer__polynomialfeatures-2__degree': np.arange(1,3,1), 'columntransformer__polynomialfeatures-3__degree': np.arange(1,4,1), } kf = KFold(n_splits=5, shuffle=True, random_state=0) grid_search_poly_model = GridSearchCV(poly_model, parameters_poly_model, cv=kf, scoring="neg_root_mean_squared_error").fit(x_train_scaled, y_train)

get_coefs(grid_search_poly_model.best_estimator_)[1:30]

alphas_poly = np.linspace(0.014, 0.015, num=10) poly_model_lasso = make_pipeline( make_column_transformer( (PolynomialFeatures(include_bias=False), ['n_directions']), (PolynomialFeatures(include_bias=False), ['n_speak_char']), (PolynomialFeatures(include_bias=False), ['Number of Cameos']), remainder = 'passthrough'), Lasso(fit_intercept = False)) parameters_poly_model_lasso = { 'columntransformer__polynomialfeatures-1__degree': np.arange(1,4,1), 'columntransformer__polynomialfeatures-2__degree': np.arange(1,3,1), 'columntransformer__polynomialfeatures-3__degree': np.arange(1,4,1), 'lasso__alpha': alphas_poly, } grid_search_poly_model_lasso = GridSearchCV(poly_model_lasso, parameters_poly_model_lasso, cv=kf, scoring="neg_root_mean_squared_error", ).fit(x_train_scaled, y_train)

get_coefs(grid_search_poly_model_lasso.best_estimator_)[1:30]

alphas = np.linspace(2, 5, num=10) poly_model_ridge = make_pipeline( make_column_transformer( (PolynomialFeatures(include_bias=False), ['n_directions']), (PolynomialFeatures(include_bias=False), ['n_speak_char']), (PolynomialFeatures(include_bias=False), ['Number of Cameos']), remainder = 'passthrough'), Ridge(fit_intercept = False)) parameters_poly_model_ridge = { 'columntransformer__polynomialfeatures-1__degree': np.arange(1,4,1), 'columntransformer__polynomialfeatures-2__degree': np.arange(1,3,1), 'columntransformer__polynomialfeatures-3__degree': np.arange(1,4,1), 'ridge__alpha': alphas, } grid_search_poly_model_ridge = GridSearchCV(poly_model_ridge, parameters_poly_model_ridge, cv=kf, scoring="neg_root_mean_squared_error", ).fit(x_train_scaled, y_train)

get_coefs(grid_search_poly_model_ridge.best_estimator_)[1:30]

x_test = test.drop(columns=['Episode Score']) y_test =test['Episode Score']

rmse_nslm = model_fit(lm, x_test, y_test, 'Non Scaled Linear Model RMSE: ', plot=False)

rmse_slm = model_fit(lm_scaled, x_test, y_test, 'Scaled Linear Model RMSE: ', plot=False)

rmse_rr_slm = model_fit(lm_scaled_ridge.best_estimator_, x_test, y_test, 'Ridge Regression on Scaled Linear Model RMSE: ', plot=False)

rmse_l_slm = model_fit(lm_scaled_lasso.best_estimator_, x_test, y_test, 'Lasso on Scaled Linear Model RMSE: ', plot=False)

rmse_std_pm = model_fit(grid_search_poly_model.best_estimator_, x_test, y_test, 'Standard Polynomial Model RMSE: ',plot=False)

rmse_l_pm = model_fit(grid_search_poly_model_lasso.best_estimator_, x_test, y_test, 'Lasso Polynomial Model RMSE: ',plot=False)

rmse_rr_pm = model_fit(grid_search_poly_model_ridge.best_estimator_, x_test, y_test, 'Ridge Polynomial Model RMSE: ',plot=False)

rmse_results ={'Method:':['Non-Scaled Linear','Scaled Linear','Ridge Scaled Linear', 'Lasso Scaled Linear','Polynomial','Ridge Polynomial','Lasso Polynomial'], 'RMSE':np.around([rmse_nslm, rmse_slm, rmse_rr_slm, rmse_l_slm, rmse_std_pm, rmse_rr_pm, rmse_l_pm], decimals = 3) } rmse_results = pd.DataFrame(rmse_results) rmse_results.sort_values(by = ['RMSE'])

coeffs

array_prop_df = pd.DataFrame(array_proportion) # Find episode in which Michael has the proportion of lines index = np.argmax(array_prop_df[[0]]) best_line = {'character:': main_chars, 'proportion' : np.around(array_prop_df.iloc[index,:],2) } best_line_df = pd.DataFrame(best_line).transpose() best_line_df

# Writer writer_coeffs = lm_scaled_ridge.best_estimator_.named_steps['ridge'].coef_[5:16] best_writer = variables[np.argmax(writer_coeffs)+5] best_writer_coeffs = max(writer_coeffs) # Director director_coeffs = lm_scaled_ridge.best_estimator_.named_steps['ridge'].coef_[16:27] best_director = variables[np.argmax(director_coeffs)+16] best_director_coeffs = max(director_coeffs) print('The best writer is', best_writer, 'with coefficient', best_writer_coeffs) print('The best director is', best_director, 'with coefficient', best_director_coeffs)

fig = plt.figure() ax = plt.axes() x = np.linspace(0, 60, 1000) ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,1]*x , color = 'red',label='n_directions'); ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,2]*x , color = 'blue',label='n_speak_char'); ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,4]*x , color = 'green',label='Episode Length'); ax.plot(x,chosen_model.iloc[1,0]+chosen_model.iloc[1,3]*x , color = 'cyan',label='Number of Cameos'); plt.ylim(8.6, 10) plt.xlim(0, max(x)) plt.legend() plt.title("Continuous Variable in their own dimensional plane") plt.xlabel("Value") plt.ylabel("Model Predicted Episode Score") plt.show()

# Run the following to render to PDF !crew install pandoc !jupyter nbconvert --to pdf project1.ipynb