<jupyter_start><jupyter_text>Demand forecasting with BigQuery and TensorFlowIn this notebook, we will develop a machine learning model to predict the demand for taxi cabs in New York.To develop the model, we will need to get historical data of taxicab usage. This data exists in BigQuery. Let's start by looking at the schema.<jupyter_code>import gcp.bigquery as bq
import pandas as pd
import numpy as np
%%bigquery schema --table "nyc-tlc:green.trips_2015"<jupyter_output><empty_output><jupyter_text>Analyzing taxicab demand Let's pull the number of trips for each day in the 2015 dataset.  We can use the BigQuery built-in Date-and-time function DAYOFYEAR (See  BigQuery query reference ) for a full list of such functions.<jupyter_code>%%sql
SELECT DAYOFYEAR(pickup_datetime) AS daynumber FROM [nyc-tlc:green.trips_2015] LIMIT 10<jupyter_output><empty_output><jupyter_text>Modular queries and Pandas dataframe Let's use the total number of trips as our proxy for taxicab demand (other reasonable alternatives are total trip_distance or total fare_amount).  It is possible to predict multiple variables using Tensorflow, but for simplicity, we will stick to just predicting the number of trips.We will give our query a name 'taxiquery' and have it use an input variable '$YEAR'. We can then invoke the 'taxiquery' by giving it a YEAR.  The to_dataframe() converts the BigQuery result into a Pandas dataframe.<jupyter_code>%%sql --module taxiquery
SELECT daynumber, COUNT(*) AS numtrips FROM
    (SELECT DAYOFYEAR(pickup_datetime) AS daynumber FROM [nyc-tlc:green.trips_$YEAR])
GROUP BY daynumber ORDER BY daynumber
trips = bq.Query(taxiquery, YEAR=2015).to_dataframe()
trips[:5]<jupyter_output><empty_output><jupyter_text>Benchmark Often, a reasonable estimate of something is its historical average. We can therefore benchmark our machine learning model against the historical average.<jupyter_code>avg = np.mean(trips['numtrips'])
print 'Just using average={0} has RMSE of {1}'.format(avg, np.sqrt(np.mean((trips['numtrips'] - avg)**2)))<jupyter_output>Just using average=54674.0994475 has RMSE of 10163.4654442<jupyter_text>The mean here is about 55,000 and the root-mean-square-error (RMSE) in this case is about 10,000. In other words, if we were to estimate that there are 55,000 taxi trips on any given day, that estimate is will be off on average by about 10,000 in either direction.  Let's see if we can do better than this -- our goal is to make predictions of taxicab demand whose RMSE is lower than 10,000.What kinds of things affect people's use of taxicabs?  Weather data We suspect that weather influences how often people use a taxi. Perhaps someone who'd normally walk to work would take a taxi if it is very cold or rainy.<NAME> has made weather observations from the US National Oceanic and Atmospheric Administration publicly available in BigQuery. Let's use that dataset and find the station number corresponding to New York's La Guardia airport.<jupyter_code>%%sql
SELECT * FROM [fh-bigquery:weather_gsod.stations]
WHERE state = 'NY' AND wban != '99999' AND name contains 'LA GUARDIA'<jupyter_output><empty_output><jupyter_text>Variables Let's pull out the minimum and maximum daily temperature (in Fahrenheit) as well as the amount of rain (in inches) for La Guardia airport.<jupyter_code>%%sql --module wxquery
SELECT DAYOFYEAR(TIMESTAMP('$YEAR'+mo+da)) daynumber,
       FIRST(DAYOFWEEK(TIMESTAMP('$YEAR'+mo+da))) dayofweek,
       MIN(min) mintemp, MAX(max) maxtemp, MAX(IF(prcp=99.99,0,prcp)) rain
FROM [fh-bigquery:weather_gsod.gsod$YEAR]
WHERE stn='725030' GROUP BY 1 ORDER BY daynumber DESC
weather = bq.Query(wxquery, YEAR=2015).to_dataframe()
weather[:5]<jupyter_output><empty_output><jupyter_text>Merge datasets Let's use Pandas to merge (combine) the taxi cab and weather datasets day-by-day.<jupyter_code>data = pd.merge(weather, trips, on='daynumber')
data[:5]<jupyter_output><empty_output><jupyter_text>Exploratory analysis Is there a relationship between maximum temperature and the number of trips?<jupyter_code>j = data.plot(kind='scatter', x='maxtemp', y='numtrips')<jupyter_output><empty_output><jupyter_text>The scatterplot above doesn't look very promising. Is there a relationship between the day of the week and the number of trips?<jupyter_code>j = data.plot(kind='scatter', x='dayofweek', y='numtrips')<jupyter_output><empty_output><jupyter_text>Hurrah, we seem to have found a predictor. It appears that people use taxis more later in the week. Perhaps New Yorkers make weekly resolutions to walk more and then lose their determination later in the week, or maybe it reflects tourism dynamics in New York City.Perhaps if we took out the confounding effect of the day of the week, maximum temperature will start to have an effect. Let's see if that's the case:<jupyter_code>j = data[data['dayofweek'] == 7].plot(kind='scatter', x='maxtemp', y='numtrips')<jupyter_output><empty_output><jupyter_text>Removing the confounding factor does seem to reflect an underlying trend around temperature. But ... the data are a little sparse, don't you think?  This is something that you have to keep in mind -- the more predictors you start to consider (here we are using two: day of week and maximum temperature), the more rows you will need so as to avoid  overfitting  the model.  Adding 2014 data Let's add in 2014 data to the Pandas dataframe.  Note how useful it was for us to modularize our queries around the YEAR.<jupyter_code>trips = bq.Query(taxiquery, YEAR=2014).to_dataframe()
weather = bq.Query(wxquery, YEAR=2014).to_dataframe()
data2014 = pd.merge(weather, trips, on='daynumber')
data2014[:5]
data2 = pd.concat([data, data2014])
data2.describe()
j = data2[data2['dayofweek'] == 7].plot(kind='scatter', x='maxtemp', y='numtrips')<jupyter_output><empty_output><jupyter_text>The data do seem a bit more robust.  If we had more data, it would be better of course. But in this case, we only have 2014 and 2015 data, so that's what we will go with.  Machine Learning with Tensorflow We'll use 80% of our dataset for training and 20% of the data for testing the model we have trained. Let's shuffle the rows of the Pandas dataframe so that this division is random.  The predictor (or input) columns will be every column in the database other than the number-of-trips (which is our target, or what we want to predict).The machine learning models that we will use -- linear regression and neural networks -- both require that the input variables are numeric in nature.The day of the week, however, is a categorical variable (i.e. Tuesday is not really greater than Monday). So, we should create separate columns for whether it is a Monday (with values 0 or 1), Tuesday, etc.Against that, we do have limited data (remember: the more columns you use as input features, the more rows you need to have in your training dataset), and it appears that there is a clear linear trend by day of the week. So, we will opt for simplicity here and use the data as-is.  Try uncommenting the code that creates separate columns for the days of the week and re-run the notebook if you are curious about the impact of this simplification.<jupyter_code>import tensorflow as tf
shuffled = data2.sample(frac=1)
# It would be a good idea, if we had more data, to treat the days as categorical variables
# with the small amount of data, we have though, the model tends to overfit
#predictors = shuffled.iloc[:,2:5]
#for day in xrange(1,8):
#  matching = shuffled['dayofweek'] == day
#  predictors.loc[matching, 'day_' + str(day)] = 1
#  predictors.loc[~matching, 'day_' + str(day)] = 0
predictors = shuffled.iloc[:,1:5]
predictors[:5]
shuffled[:5]
targets = shuffled.iloc[:,5]
targets[:5]<jupyter_output><empty_output><jupyter_text>Let's update our benchmark based on the 80-20 split and the larger dataset.<jupyter_code>trainsize = int(len(shuffled['numtrips']) * 0.8)
avg = np.mean(shuffled['numtrips'][:trainsize])
rmse = np.sqrt(np.mean((shuffled['numtrips'][trainsize:] - avg)**2))
print 'Just using average={0} has RMSE of {1}'.format(avg, rmse)<jupyter_output>Just using average=47147.5206422 has RMSE of 12732.1343487<jupyter_text>Linear regression with Tensorflow The key statements here are the model that we are creating:model = (tf.matmul(feature_data, weights) + biases) * SCALE_NUM_TRIPSwhich is a linear model and the error measure that we are minimizing:cost = tf.nn.l2_loss(model - target_data)which is the square-error (so what we are doing is least-squares linear regression).The rest of this is pretty boiler plate.  The scaling is so that the optimization goes a lot better (when the weights are small numbers, things are more likely to converge, i.e., less likely to blow up).  We save the weights into /tmp/trained_model and display the root mean square error on the test dataset.<jupyter_code>SCALE_NUM_TRIPS = 100000
trainsize = int(len(shuffled['numtrips']) * 0.8)
testsize = len(shuffled['numtrips']) - trainsize
npredictors = len(predictors.columns)
noutputs = 1
numiter = 10000
modelfile = '/tmp/trained_model'
with tf.Session() as sess:
  feature_data = tf.placeholder("float", [None, npredictors])
  target_data = tf.placeholder("float", [None, noutputs])
  weights = tf.Variable(tf.truncated_normal([npredictors, noutputs], stddev=0.01))
  biases = tf.Variable(tf.ones([noutputs]))
  model = (tf.matmul(feature_data, weights) + biases) * SCALE_NUM_TRIPS
  cost = tf.nn.l2_loss(model - target_data)
  training_step = tf.train.AdamOptimizer(learning_rate=0.0001).minimize(cost)
  init = tf.initialize_all_variables()
  sess.run(init)

  saver = tf.train.Saver({'weights' : weights, 'biases' : biases})
  for iter in xrange(0, numiter):
    sess.run(training_step, feed_dict = {
        feature_data : predictors[:trainsize].values,
        target_data : targets[:trainsize].values.reshape(trainsize, noutputs)
      })
    
    if iter%1000 == 0:
      print '{0} error={1}'.format(iter, np.sqrt(cost.eval(feed_dict = {
          feature_data : predictors[:trainsize].values,
          target_data : targets[:trainsize].values.reshape(trainsize, noutputs)
      }) / trainsize))
    
  filename = saver.save(sess, modelfile, global_step=numiter)
  print 'Model written to {0}'.format(filename)

  print 'testerror={0}'.format(np.sqrt(cost.eval(feed_dict = {
          feature_data : predictors[trainsize:].values,
          target_data : targets[trainsize:].values.reshape(testsize, noutputs)
      }) / testsize))<jupyter_output>0 error=37226.769573
1000 error=11397.6734661
2000 error=10759.6089285
3000 error=10037.6474454
4000 error=9306.57428869
5000 error=8646.24931995
6000 error=8103.11729799
7000 error=7699.99690218
8000 error=7440.12137614
9000 error=7305.52026543
Model written to /tmp/trained_model-10000
testerror=8386.6236006<jupyter_text>The test error of about 8400 indicates that we are doing better with the machine learning model than we would be if we were to just use the historical average (our benchmark).  Neural network with Tensorflow To change from linear regression to neural network regression, change the model to have more layers, and insert functional transformations of the linear combinations (here we do relu, but we could also do tanh or sigmoid):model = (tf.matmul(tf.nn.relu(tf.matmul(feature_data, weights1) + biases1), weights2) + biases2) * SCALE_NUM_TRIPSThere are more weights, which need to be saved and retrieved, but the boilerplate code remains mostly the same as the linear regression example.<jupyter_code>SCALE_NUM_TRIPS = 100000
trainsize = int(len(shuffled['numtrips']) * 0.8)
testsize = len(shuffled['numtrips']) - trainsize
npredictors = len(predictors.columns)
noutputs = 1
nhidden = 5
numiter = 10000
modelfile = '/tmp/trained_model'
with tf.Session() as sess:
  feature_data = tf.placeholder("float", [None, npredictors])
  target_data = tf.placeholder("float", [None, noutputs])
  weights1 = tf.Variable(tf.truncated_normal([npredictors, nhidden], stddev=0.01))
  weights2 = tf.Variable(tf.truncated_normal([nhidden, noutputs], stddev=0.01))
  biases1 = tf.Variable(tf.ones([nhidden]))
  biases2 = tf.Variable(tf.ones([noutputs]))
  model = (tf.matmul(tf.nn.relu(tf.matmul(feature_data, weights1) + biases1), weights2) + biases2) * SCALE_NUM_TRIPS
  cost = tf.nn.l2_loss(model - target_data)
  training_step = tf.train.AdamOptimizer(learning_rate=0.0001).minimize(cost)
  init = tf.initialize_all_variables()
  sess.run(init)

  saver = tf.train.Saver({'weights1' : weights1, 'biases1' : biases1, 'weights2' : weights2, 'biases2' : biases2})
  for iter in xrange(0, numiter):
    sess.run(training_step, feed_dict = {
        feature_data : predictors[:trainsize].values,
        target_data : targets[:trainsize].values.reshape(trainsize, noutputs)
      })
    
    if iter%1000 == 0:
      print '{0} error={1}'.format(iter, np.sqrt(cost.eval(feed_dict = {
          feature_data : predictors[:trainsize].values,
          target_data : targets[:trainsize].values.reshape(trainsize, noutputs)
      }) / trainsize))
    
  filename = saver.save(sess, modelfile, global_step=numiter)
  print 'Model written to {0}'.format(filename)

  print 'testerror={0}'.format(np.sqrt(cost.eval(feed_dict = {
          feature_data : predictors[trainsize:].values,
          target_data : targets[trainsize:].values.reshape(testsize, noutputs)
      }) / testsize))<jupyter_output>0 error=38225.6868073
1000 error=10581.3908453
2000 error=8598.12454908
3000 error=7775.54995489
4000 error=7289.19077071
5000 error=7166.95208592
6000 error=7115.87436393
7000 error=7073.102492
8000 error=7059.94672122
9000 error=7048.97842816
Model written to /tmp/trained_model-10000
testerror=8181.71200578<jupyter_text>Using a neural network does seem to get us an additional boost in performance (8400 to 8200 when I ran it -- because of random seeds, your actual answers might be different).  Running a trained model So, we have trained a model, and saved it to a file. Let's use this model to predict taxicab demand given the expected weather for three days.Here we make a Dataframe out of those inputs, load up the saved model (note that we have to know the model equation -- it's not saved in the model file) and use it to predict the taxicab demand.<jupyter_code>input = pd.DataFrame.from_dict(data = 
                               {'dayofweek' : [4, 5, 6],
                                'mintemp' : [50, 36, 54],
                                'maxtemp' : [80, 40, 75],
                                'rain' : [0, 0.8, 0.3]})
with tf.Session() as sess:
    filename = modelfile + '-' + str(numiter)
    saver = tf.train.Saver({'weights1' : weights1, 'biases1' : biases1, 'weights2' : weights2, 'biases2' : biases2})
    saver.restore(sess, filename)
    feature_data = tf.placeholder("float", [None, npredictors])
    predict_operation = (tf.matmul(tf.nn.relu(tf.matmul(feature_data, weights1) + biases1), weights2) + biases2) * SCALE_NUM_TRIPS
    predicted = sess.run(predict_operation, feed_dict = {
        feature_data : input.values
      })
    print predicted<jupyter_output>[[ 32465.81054688]
 [ 49777.13671875]
 [ 41207.85546875]]
{'fraction_non_alphanumeric': 0.09770760388452138, 'fraction_numerical': 0.04677280835039968, 'mean_word_length': 4.865168539325842, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '93093', 'license': 'Apache-2.0', 'n_tokens_mistral': 5200, 'n_tokens_neox': 4657, 'n_words': 1639, 'path': 'CPB100/lab4a/demandforecast.ipynb', 'repo_name': 'rituc/training-data-analyst', 'size': '15151'}