<jupyter_start><jupyter_text>Hello professor,

PROJECT:

I did some experiments on the data. Before I tell you the result, I would like to say what I have done.
1- I put a filter to consider the data between 2019-01-01 and 2020-03-16 to avoid the covid effects on the model.
2- Consequently, I changed the STC's, except for our usual constraints, I added a filter that cancels out the STC days which fall in the covid time (after 2020-03-16).
3- I changed the days, months, DOW, and LU categories to ONE-HOT vectors. I found out that the days, months, and DOW also should be considered categorical data. In our case, increasing or decreasing their number would not have any numerical order. Assume we have a function that takes day number as one of the inputs and then give the possible salary. Here, the number of days or months have meaning, but in our case, it is not like this due to my understanding. So, I treat them as categorical data. FOR EXAMPLE, 001000000000 mean month 3.
4 - The total number of valid data points is 4946. I divided to test and train with the ratio of 20/80 and then split to train and validation set with the ratio of 15/85. For the test set, the samples are in order in their date for 15 specific sites to keep track of the AADP for further tests. The validation set is chosen randomly.
5- I have two sets of features. One of them is based on the LU categories, and the other is based on LU characteristics. LU category has three types (residential, employment and commercial area, and others). LU characteristics are the number of shops, banks, etc.
6 - For the first attempt to select the best set of features, I used ridge regression.
7- I used the variance inflation factor (VIF) to drop the features that have colinearity with others. You are right that we can not have LU categories and LU characteristics together. Although I set the VIF threshold to be 5, in one case (Considering the samples before covid), I changed the VIF threshold to 6.5 from 5 to not drop the LU categories. Other VIF's are all less than 5, and due to my search, the best threshold is 5, but features less than 10 could also work without problem in the data set. In Think 6.3, for LU categories in the pre-covid data set would not be problematic.
8- DOW is 0 for weekdays and 1 for weekends. I included the days in my models.
9- When I changed the STC's days and canceled out after covid days, the STC's make more sense now. I added the new STC's columns in the attached excel file(arizona.xlsx)

REGRESSION RESULTS:
Using all of the data: R2 = 0.32, MSE = 2000
Using data before pandemic(feature set one= using LU
  NN<jupyter_code>import keras
from keras.datasets import cifar10
from keras.preprocessing.image import ImageDataGenerator
from keras.models import Sequential, load_model
from keras.layers import Dense, Dropout, Activation, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras.callbacks import ModelCheckpoint
import math
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler

# parameters for this script
batch_size = 12
epochs = 120

from sklearn.model_selection import train_test_split

# FESTURES SET ONE
X1 = data_features.drop(data_features.columns[0:3],1)

# FESTURES SET TWO
X2 = data_features2.drop(data_features2.columns[0:3],1)
X2 = X2.drop(X2.columns[3],1)

# Normalizing the data X1 and X2
# X1
scaler = MinMaxScaler()
X1 = scaler.fit_transform(X1)
X2 = scaler.fit_transform(X2)

X1 = pd.DataFrame(X1)
X2 = pd.DataFrame(X2)





Y = pd.DataFrame(data['PEDESTRIAN COUNTS'])

X1_fit = X1.drop(X1.index[test_stations_ind])
Y_fit = Y.drop(Y.index[test_stations_ind])

X1_test = pd.DataFrame(X1, index=test_stations_ind )
Y_test = pd.DataFrame(Y,index = test_stations_ind)


X2_fit = X2.drop(X2.index[test_stations_ind])
X2_test = pd.DataFrame(X2, index=test_stations_ind )






X1
model = Sequential()

model.add(Dense(12, input_dim=X1.shape[1], kernel_initializer='he_uniform'))
model.add(Activation('sigmoid'))
model.add(Dropout(0))








model.add(Dense(1))
model.add(Activation('linear'))
print(model.summary())
# initiate Adam optimizer
opt = keras.optimizers.Adam()
opt2 = keras.optimizers.SGD(lr=0.01, momentum=0.9)
# Compile the model before using it
model.compile(loss='mean_squared_error',
              optimizer=opt,
              metrics=['accuracy','mse'] ,)
print(model.summary())

# create a callback that will save the best model while training
save_best_model = ModelCheckpoint('best_model3.h5', monitor='val_mse', mode='min', save_best_only=True, verbose=1, )

history = model.fit(X1_fit.values, Y_fit.values,
                  batch_size=batch_size,
                  epochs=epochs,
                  validation_data=(X1_test.values, Y_test.values),
                  
                  shuffle=True,
                  callbacks=[save_best_model])

# history = model.fit(X1.values, Y.values,
#                   batch_size=batch_size,
#                   epochs=epochs,
#                   validation_split= 0.2 ,
#                   shuffle=True,
#                   callbacks=[save_best_model], )


# Plot training accuracy
plt.plot(history.history['mse'])
plt.title('training accuracy')
plt.ylabel('training accuracy')
plt.xlabel('epoch')
plt.legend(['Conv'], loc='upper left')
plt.show()

# Plot validation accuracy
plt.plot(history.history['val_mse'])
plt.title('validation accuracy')
plt.ylabel('validation accuracy')
plt.xlabel('epoch')
plt.legend(['Conv'], loc='upper left')
plt.show()

# Evaluate the best model saved (i.e., model with best validation accuracy) on the test set
saved_model = load_model('best_model3.h5')
scores = saved_model.evaluate(X1_test.values, Y_test.values, verbose=1)
print('Test accuracy:', scores)
from sklearn.linear_model import Ridge
from sklearn.metrics import r2_score
def ridge_regression(xtrain,ytrain,xval,yval, alpha):
    #Fit the model
    ridgreg = Ridge(alpha=alpha,normalize=False, max_iter=1e5)
    ridgreg.fit(xtrain,ytrain)
    y_pred = ridgreg.predict(xval)

    #Return the result in pre-defined format
    mss = sum((y_pred-yval)**2) / y_pred.shape[0]
    r2 = r2_score(yval, y_pred)

    return mss,r2,yval,y_pred
mss,r2,yval,ypred = ridge_regression(X1_fit.values,Y_fit.values,X1_test.values,Y_test.values, 20)
print(mss),print(r2)
saved_model = load_model('best_model.h5')
pred = saved_model.predict(X1_test.values, verbose=1)
r2 = r2_score(Y_test,pred)
r2

X2_fit<jupyter_output><empty_output><jupyter_text>Gaussian Process Regression<jupyter_code>def predict_gaussian_process(inputs, posterior):

  mean , variance = posterior(inputs)
  return mean, variance
def eval_gaussian_process(inputs, posterior, targets):
  mean , variance = posterior(inputs)
  diff = mean-targets
  mean_squared_error = np.mean(diff**2)

  return mean_squared_error
def train_gaussian_process(train_inputs, train_targets, measurement_variance, kernel, kernel_param):

  inverse_regularized_gram_matrix = np.linalg.inv(kernel(train_inputs,train_inputs,kernel_param)+measurement_variance*np.eye(train_inputs.shape[0]))
  
  # define a function that computes a posterior distribution over the outputs given some inputs
  # inputs: 
  #     inputs: matrix of input points for which we would like to predict a Gaussian over the outputs (numpy array of N' data points x M features)
  #     train_inputs: matrix of train inputs that serve as evidence for the posterior (numpy array of N data points x M features)
  #     train_targets: vector of train targets that serve as evidence for the posterior (numpy array of N targets)
  #     inverse_regularized_gram_matrix: (K + sigma^2 I)^{-1}
  # outputs:
  #     mean: vector of Gaussian means for each prediction (numpy array of N' means)
  #     variance: vector of Gaussian variances for each prediction (numpy array of N' variances)
  def posterior(inputs, train_inputs=train_inputs, train_targets=train_targets, inverse_regularized_gram_matrix=inverse_regularized_gram_matrix, kernel=kernel):
    mean = kernel(inputs, train_inputs, kernel_param) @ inverse_regularized_gram_matrix @ train_targets
    variance = kernel(inputs, inputs, kernel_param) + (measurement_variance ** 2) * np.eye(inputs.shape[0]) - kernel(inputs, train_inputs, kernel_param) @ inverse_regularized_gram_matrix @ kernel(train_inputs, inputs, kernel_param)

    return mean, variance

  return posterior
def identity_kernel(inputs1,inputs2,dummy_param=None):
  gram_matrix = inputs1 @ inputs2.T
  return gram_matrix
def polynomial_kernel(inputs1,inputs2,degree):
  gram_matrix = (inputs1 @ inputs2.T + 1) ** degree
  return gram_matrix
def gaussian_kernel(inputs1,inputs2,width):

    gram_matrix = np.zeros([inputs1.shape[0], inputs2.shape[0]], dtype=float)
    nd = inputs2.shape[0]
    for i in range(nd):
        diff2 = np.sum((inputs1 - inputs2[i, :]) ** 2, axis=-1)
        gram_matrix[:, i] = np.exp(-diff2 / (2 * width ** 2))
  
    return gram_matrix
def cross_validation_gaussian_process(k_folds, kernel, hyperparameters, inputs, targets, measurement_variance):

# making validation and test matrices
  fold_num = k_folds
  validation_inputs = np.zeros( ( inputs.shape[0]//k_folds , inputs.shape[1] , k_folds ) )
  train_inputs_ = np.zeros((  inputs.shape[0] - inputs.shape[0]//k_folds , inputs.shape[1] , k_folds  ))
  validation_targets = np.zeros(( inputs.shape[0]//k_folds , k_folds))
  train_targets_ = np.zeros(( inputs.shape[0] - inputs.shape[0]//k_folds , k_folds))
  # filling validation and test matrices

  for i in range(0 , fold_num):
   validation_inputs [ : , : , i] = inputs [  i * inputs.shape[0] // k_folds : (inputs.shape[0]//k_folds * (i + 1)) , 0 : inputs.shape[1]]
   validation_targets [: , i] = targets [( i * inputs.shape[0] // k_folds) : (train_inputs.shape[0]//k_folds) * (i + 1)]

   train_inputs_[ : , : , i] = np.delete ( inputs , np.s_[( i * inputs.shape[0] // k_folds) : (inputs.shape[0]//k_folds) * (i + 1)] , axis=0  )
   train_targets_[: , i] = np.delete ( targets , np.s_[ ( i * inputs.shape[0] // k_folds) : (inputs.shape[0]//k_folds) * (i + 1) ] , axis = 0 )

  mse_matrix = np.zeros(( k_folds , len(hyperparameters) ))    # filling the accuracies matrix. x axes is # of hyperparameters (neighbours) and y axes is each fold of our train_inputs
  
  for ii in range(k_folds):
    for j in range(len(hyperparameters)):
     posterior = train_gaussian_process(train_inputs_[:,:,ii], train_targets_[:,ii], measurement_variance, kernel, hyperparameters[j] )
     mse_matrix[ii,j] = eval_gaussian_process(validation_inputs[:,:,ii], posterior, validation_targets[:,ii])

  mean_squared_errors = mse_matrix.mean(0)
  best_hyperparam = hyperparameters[np.argmin(mean_squared_errors)]
  best_mean_squared_error = np.min(mean_squared_errors)


  return best_hyperparam, best_mean_squared_error, mean_squared_errors
def plot_gaussian_process_mean_squared_errors(mean_squared_errors,hyperparams,xlabel):
  plt.plot(hyperparams,mean_squared_errors)
  plt.ylabel('mean squared error')
  plt.xlabel(xlabel)
  plt.show()
# load data
train_inputs = X1_fit.values
train_targets = Y_fit.values.flatten()
test_inputs = X1_test.values
test_targets = Y_test.values.flatten()

# likelihood parameter (see definition of train_gaussian_process for more details)
measurement_variance = 1

# number of folds in k-fold cross validation
k_folds = 9

# identity kernel
print('\nIdentity Kernel')
posterior = train_gaussian_process(train_inputs, train_targets, measurement_variance, identity_kernel, None)
mean_squared_error = eval_gaussian_process(test_inputs, posterior, test_targets)
print('test mean squared error with identity kernel: ' + str(mean_squared_error))

# # polynomial kernel
# print('\nPolynomial Kernel')
# # degrees to be evaluated by cross validation
# hyperparams = range(1,21)
# best_degree, best_mean_squared_error, mean_squared_errors = \
#   cross_validation_gaussian_process(k_folds,polynomial_kernel,hyperparams,train_inputs,train_targets,measurement_variance)
# plot_gaussian_process_mean_squared_errors(mean_squared_errors,hyperparams,'degree')
# print('best degree: ' + str (best_degree))
# print('best cross validation mean squared error: ' + str(best_mean_squared_error))
# posterior = train_gaussian_process(train_inputs, train_targets, measurement_variance, polynomial_kernel, best_degree)
# mean_squared_error = eval_gaussian_process(test_inputs, posterior, test_targets)
# print('test mean squared error with polynomial kernel: ' + str(mean_squared_error))

# gaussian kernel
print('\nGaussian Kernel')
# gaussian widths to be evaluated by cross validation
hyperparams = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
best_width, best_mean_squared_error, mean_squared_errors = \
  cross_validation_gaussian_process(k_folds,gaussian_kernel,hyperparams,train_inputs,train_targets,measurement_variance)
plot_gaussian_process_mean_squared_errors(mean_squared_errors,hyperparams,'width')
print('best width: ' + str (best_width))
print('best cross validation mean squared error: ' + str(best_mean_squared_error))
posterior = train_gaussian_process(train_inputs, train_targets, measurement_variance, gaussian_kernel, best_width)
mean_squared_error = eval_gaussian_process(test_inputs, posterior, test_targets)
print('test mean squared error with Gaussian kernel: ' + str(mean_squared_error))<jupyter_output>Identity Kernel
test mean squared error with identity kernel: 1860.3122794471678

Gaussian Kernel
{'fraction_non_alphanumeric': 0.09518248175182482, 'fraction_numerical': 0.019854014598540148, 'mean_word_length': 4.452049343414246, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '456640', 'license': 'Unlicense', 'n_tokens_mistral': 4648, 'n_tokens_neox': 4338, 'n_words': 1266, 'path': 'GPR_LR_NN.ipynb', 'repo_name': 'Amir-mjafari/AADP-ESTIMATION', 'size': '13700'}