Due Friday, May 10 @ 11:59PM Wednesday, May 15 @ 11:59PM
People often think linear regression on a dataset with millions of observations is a big data problem. Now we learnt various methods for solving linear regression and should realize that, with right choice of algorithm, it is a problem that can be handled by any moderate computer.
Download the flight data from http://stat-computing.org/dataexpo/2009/the-data.html. For this exercise, we only need data from years 2003-2008. If you are using Mac or Linux, you can run the following Julia script to download and unzip files for all years.
;cat getflightdata.jl
Please do not put data files in Git. For grading purpose (reproducibility), we will assume that data files are in the same directory as the Jupyter notebook file.
Find out how many data points are in each year.
We are interested in how the time gain of a flight, defined as DepDelay - ArrDelay, depends on the distance traveled (Distance), departure delay (DepDelay), and carrier (UniqueCarrier).
We want to fit a linear regression Gain ~ 1 + Distance + DepDelay + UniqueCarrier using data from 2003-2008. Note UniqueCarrier is a factor with 23 levels: "9E", "AA", "AQ", "AS", "B6", "CO", "DH", "DL", "EV", "F9", "FL", "HA", "HP", "MQ", "NW", "OH", "OO", "TZ", "UA", "US", "WN", "XE", and "YV". We use the dummy coding with "9E" as base level.
Will the design matrix $\mathbf{X}$ (in double precision) fit into the memory of you computer?
Review the Summary of Linear Regression and choose one method in the table to solve the linear regression.
Report the estimated regression coefficients $\widehat \beta$, estimated variance $\widehat \sigma^2 = \sum_i (y_i - \widehat y_i)^2 / (n - p)$, and standard errors of $\widehat \beta$.
Hint: It took my laptop less than 3 minutes to import data and fit linear regression.
Go to your resume/cv and claim you have experience performing analytics on data with millions of observations.
Following code explores the data in 2003 and generates the design matrix and responses for that year. Feel free to use the code in your solution.
versioninfo()
Print first 10 lines of 2003 data.
;head 2003.csv
# how many data points
countlines("2003.csv")
# import data from csv
using JuliaDB
@time yrtable = loadtable("2003.csv",
datacols = ["DepDelay", "ArrDelay", "UniqueCarrier", "Distance"])
# drop rows with missing values
yrtable = dropmissing(yrtable)
# mapping from variable names to X columns
# carrier "9E" is used as base level
const var2col = Dict(
"Intercept" => 1,
"Distance" => 2,
"DepDelay" => 3,
"AA" => 4,
"AQ" => 5,
"AS" => 6,
"B6" => 7,
"CO" => 8,
"DH" => 9,
"DL" => 10,
"EV" => 11,
"F9" => 12,
"FL" => 13,
"HA" => 14,
"HP" => 15,
"MQ" => 16,
"NW" => 17,
"OH" => 18,
"OO" => 19,
"TZ" => 20,
"UA" => 21,
"US" => 22,
"WN" => 23,
"XE" => 24,
"YV" => 25,
"Gain" => 26)
# mapping from column to variable names
const col2var = Dict(value => key for (key, value) in var2col)
# a custom function to generate [X y] from data table
function generate_xy(tbl::IndexedTable)
# X matrix
XY = zeros(length(tbl), 26)
# intercept term
XY[:, 1] .= 1
# Distance term
XY[:, 2] = columns(tbl, :Distance)
# DepDelay term
XY[:, 3] = columns(tbl, :DepDelay)
# Dummy coding for airline
@inbounds for i in 1:length(tbl)
tbl[i][:UniqueCarrier] == "9E" && continue # base level
XY[i, var2col[tbl[i][:UniqueCarrier]]] = 1
end
# last column is response: gain = depdelay - arrdelay
XY[:, 26] = select(tbl,
(:DepDelay, :ArrDelay) => p -> Float64(p.DepDelay - p.ArrDelay))
# return
XY
end
@time xy = generate_xy(yrtable)
We are going to try different numerical methods learnt in class on the Google PageRank problem.
Let $\mathbf{A} \in \{0,1\}^{n \times n}$ be the connectivity matrix of $n$ web pages with entries $$ \begin{eqnarray*} a_{ij}= \begin{cases} 1 & \text{if page $i$ links to page $j$} \\ 0 & \text{otherwise} \end{cases}. \end{eqnarray*} $$ $r_i = \sum_j a_{ij}$ is the out-degree of page $i$. That is $r_i$ is the number of links on page $i$. Imagine a random surfer exploring the space of $n$ pages according to the following rules.
The process defines a Markov chain on the space of $n$ pages. Write the transition matrix $\mathbf{P}$ of the Markov chain as a sparse matrix plus rank 1 matrix.
According to standard Markov chain theory, the (random) position of the surfer converges to the stationary distribution $\mathbf{x} = (x_1,\ldots,x_n)^T$ of the Markov chain. $x_i$ has the natural interpretation of the proportion of times the surfer visits page $i$ in the long run. Therefore $\mathbf{x}$ serves as page ranks: a higher $x_i$ means page $i$ is more visited. It is well-known that $\mathbf{x}$ is the left eigenvector corresponding to the top eigenvalue 1 of the transition matrix $\mathbf{P}$. That is $\mathbf{P}^T \mathbf{x} = \mathbf{x}$. Therefore $\mathbf{x}$ can be solved as an eigen-problem. Show that it can also be cast as solving a linear system. Since the row sums of $\mathbf{P}$ are 1, $\mathbf{P}$ is rank deficient. We can replace the first equation by the $\sum_{i=1}^n x_i = 1$.
Obtain the connectivity matrix A from the SNAP/web-Google data in the MatrixDepot package.
using MatrixDepot
md = mdopen("SNAP/web-Google")
md.A
# I found I need to run following line at end of notebook to avoid a bug
Base.Filesystem.rm(dirname(pathof(MatrixDepot)) * "/../data/db.data", force=true)
Compute summary statistics:
Hint: For plots, you can use the UnicodePlots.jl package.
Set the teleportation parameter at $p = 0.85$. Try the following methods to solve the PageRank problem using the SNAP/web-Google data.
For iterative methods, you can use the IterativeSolvers.jl and Arpack.j packages. Make sure to utilize the special structure of $\mathbf{P}$ (sparse + rank 1) to speed up the matrix-vector multiplication. (Hint: The LinearMaps.jl packages exists exactly for this purpose.)
List the top 20 nodes you found.
As of Monday Apr 29, 2019, there are at least 5.49 billion indexed webpages on internet according to http://www.worldwidewebsize.com/. Explain whether each of these methods works for the PageRank problem at this scale.
Go to your resume/cv and claim you have experience performing analysis on a network of one million nodes.