Due Friday, May 25 @ 11:59PM
Those who took my 203B: Introduction to Data Science last quarter had a (painful) experience of wrangling an Apache Spark cluster to do linear regression on a dataset with more than 100 million observations. 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 Bash script, which downloads and unzips files for all years.
# Download flight data by year
for i in {1987..2008}
do
echo "$(date) $i Download"
fnam=$i.csv.bz2
wget -O ./$fnam http://stat-computing.org/dataexpo/2009/$fnam
echo "$(date) $i unzip"
bzip2 -d ./$fnam
done
# Download airline carrier data
wget -O ./airlines.csv http://www.transtats.bts.gov/Download_Lookup.asp?Lookup=L_UNIQUE_CARRIERS
# Download airports data
wget -O ./airports.csv https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat
Find out how many data points 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 (in double precision) fit into the memory of you computer?
Review the Summary of Linear Regression and devise a strategy 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 - 1)$, and coefficient standard errors.
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 hundred 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.
;head 2003.csv
# how many data points
countlines("2003.csv")
# import data from csv
using JuliaDB
# only need columns: DepDelay, ArrDelay, UniqueCarrier, Distance
@time yrtable = loadtable(
"2003.csv",
datacols = ["DepDelay", "ArrDelay", "UniqueCarrier", "Distance"])
# drop rows with missing values
yrtable = dropna(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 = map(reverse, var2col)
# a custom function to generate [X y] from data table
function generate_xy(tbl::NextTable)
# X matrix
XY = zeros(length(tbl), 26)
# intercept term
@views fill!(XY[:, 1], 1)
# Distance term
@views copy!(XY[:, 2], columns(tbl, :Distance))
# DepDelay term
@views copy!(XY[:, 3], columns(tbl, :DepDelay))
# Dummy coding for airline
@inbounds for i in 1:length(tbl)
yrtable[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$.
Download the ucla.zip package from course webpage. Unzip the package, which contains two files U.txt and A.txt. U.txt lists the 500 URL names. A.txt is the $500 \times 500$ connectivity matrix. Read data into Julia. Compute summary statistics:
Set the teleportation parameter at $p = 0.85$. Try the following methods to solve for $\mathbf{x}$ using the ucla.zip data.
For iterative methods, you can use the IterativeSolvers.jl package. Make sure to utilize the special structure of $\mathbf{P}$ (sparse + rank 1) to speed up the matrix-vector multiplication.
List the top 20 ranked URLs you found.
As of Monday May 11 2018, there are at least 1.83 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.