Technical – Page 3 – Anirudh Jayaraman

Data Manipulation in R with dplyr – Part 3

This happens to be my 50th blog post – and my blog is 8 months old.

🙂

This post is the third and last post in in a series of posts (Part 1 – Part 2) on data manipulation with dlpyr. Note that the objects in the code may have been defined in earlier posts and the code in this post is in continuation with code from the earlier posts.

Although datasets can be manipulated in sophisticated ways by linking the 5 verbs of dplyr in conjunction, linking verbs together can be a bit verbose.

Creating multiple objects, especially when working on a large dataset can slow you down in your analysis. Chaining functions directly together into one line of code is difficult to read. This is sometimes called the Dagwood sandwich problem: you have too much filling (too many long arguments) between your slices of bread (parentheses). Functions and arguments get further and further apart.

The %>% operator allows you to extract the first argument of a function from the arguments list and put it in front of it, thus solving the Dagwood sandwich problem.

	# %>% OPERATOR ———————————————————————-

	# with %>% operator
	hflights %>%
	mutate(diff = TaxiOut – TaxiIn) %>%
	filter(!is.na(diff)) %>%
	summarise(avg = mean(diff))

	# without %>% operator
	# arguments get further and further apart
	summarize(filter(mutate(hflights, diff = TaxiOut – TaxiIn),!is.na(diff)),
	avg = mean(diff))


	# with %>% operator
	d <- hflights %>%
	select(Dest, UniqueCarrier, Distance, ActualElapsedTime) %>%
	mutate(RealTime = ActualElapsedTime + 100, mph = Distance/RealTime*60)

	# without %>% operator
	d <- mutate(select(hflights, Dest, UniqueCarrier, Distance, ActualElapsedTime),
	RealTime = ActualElapsedTime + 100, mph = Distance/RealTime*60)

	# Filter and summarise d
	d %>%
	filter(!is.na(mph), mph < 70) %>%
	summarise(n_less = n(), n_dest = n_distinct(Dest),
	min_dist = min(Distance), max_dist = max(Distance))

	# Let's define preferable flights as flights that are 150% faster than driving,
	# i.e. that travel 105 mph or greater in real time. Also, assume that cancelled or
	# diverted flights are less preferable than driving.


	# ADVANCED PIPING EXERCISES
	# Use one single piped call to print a summary with the following variables:

	# n_non – the number of non-preferable flights in hflights,
	# p_non – the percentage of non-preferable flights in hflights,
	# n_dest – the number of destinations that non-preferable flights traveled to,
	# min_dist – the minimum distance that non-preferable flights traveled,
	# max_dist – the maximum distance that non-preferable flights traveled

	hflights %>%
	mutate(RealTime = ActualElapsedTime + 100, mph = Distance/RealTime*60) %>%
	filter(mph < 105 \| Cancelled == 1 \| Diverted == 1) %>%
	summarise(n_non = n(), p_non = 100*n_non/nrow(hflights), n_dest = n_distinct(Dest),
	min_dist = min(Distance), max_dist = max(Distance))

	# Use summarise() to create a summary of hflights with a single variable, n,
	# that counts the number of overnight flights. These flights have an arrival
	# time that is earlier than their departure time. Only include flights that have
	# no NA values for both DepTime and ArrTime in your count.

	hflights %>%
	mutate(overnight = (ArrTime < DepTime)) %>%
	filter(overnight == TRUE) %>%
	summarise(n = n())

view raw pipingOperator.r hosted with ❤ by GitHub

group_by()

group_by() defines groups within a data set. Its influence becomes clear when calling summarise() on a grouped dataset. Summarizing statistics are calculated for the different groups separately.

	# group_by() ————————————————————————-

	# Generate a per-carrier summary of hflights with the following variables: n_flights,
	# the number of flights flown by the carrier; n_canc, the number of cancelled flights;
	# p_canc, the percentage of cancelled flights; avg_delay, the average arrival delay of
	# flights whose delay does not equal NA. Next, order the carriers in the summary from
	# low to high by their average arrival delay. Use percentage of flights cancelled to
	# break any ties. Which airline scores best based on these statistics?

	hflights %>%
	group_by(UniqueCarrier) %>%
	summarise(n_flights = n(), n_canc = sum(Cancelled), p_canc = 100*n_canc/n_flights,
	avg_delay = mean(ArrDelay, na.rm = TRUE)) %>% arrange(avg_delay)

	# Generate a per-day-of-week summary of hflights with the variable avg_taxi,
	# the average total taxiing time. Pipe this summary into an arrange() call such
	# that the day with the highest avg_taxi comes first.

	hflights %>%
	group_by(DayOfWeek) %>%
	summarize(avg_taxi = mean(TaxiIn + TaxiOut, na.rm = TRUE)) %>%
	arrange(desc(avg_taxi))

view raw group_by.R hosted with ❤ by GitHub

Combine group_by with mutate

group_by() can also be combined with mutate(). When you mutate grouped data, mutate() will calculate the new variables independently for each group. This is particularly useful when mutate() uses the rank() function, that calculates within group rankings. rank() takes a group of values and calculates the rank of each value within the group, e.g.

rank(c(21, 22, 24, 23))

has output

[1] 1 2 4 3

As with arrange(), rank() ranks values from the largest to the smallest and this behaviour can be reversed with the desc() function.

	# Combine group_by with mutate—–

	# First, discard flights whose arrival delay equals NA. Next, create a by-carrier
	# summary with a single variable: p_delay, the proportion of flights which are
	# delayed at arrival. Next, create a new variable rank in the summary which is a
	# rank according to p_delay. Finally, arrange the observations by this new rank
	hflights %>%
	filter(!is.na(ArrDelay)) %>%
	group_by(UniqueCarrier) %>%
	summarise(p_delay = sum(ArrDelay >0)/n()) %>%
	mutate(rank = rank(p_delay)) %>%
	arrange(rank)

	# n a similar fashion, keep flights that are delayed (ArrDelay > 0 and not NA).
	# Next, create a by-carrier summary with a single variable: avg, the average delay
	# of the delayed flights. Again add a new variable rank to the summary according to
	# avg. Finally, arrange by this rank variable.
	hflights %>%
	filter(!is.na(ArrDelay), ArrDelay > 0) %>%
	group_by(UniqueCarrier) %>%
	summarise(avg = mean(ArrDelay)) %>%
	mutate(rank = rank(avg)) %>%
	arrange(rank)

	# Advanced group_by exercises——————————————————-

	# Which plane (by tail number) flew out of Houston the most times? How many times?
	# Name the column with this frequency n. Assign the result to adv1. To answer this
	# question precisely, you will have to filter() as a final step to end up with only
	# a single observation in adv1.
	# Which plane (by tail number) flew out of Houston the most times? How many times? adv1
	adv1 <- hflights %>%
	group_by(TailNum) %>%
	summarise(n = n()) %>%
	filter(n == max(n))

	# How many airplanes only flew to one destination from Houston? adv2
	# How many airplanes only flew to one destination from Houston?
	# Save the resulting dataset in adv2, that contains only a single column,
	# named nplanes and a single row.
	adv2 <- hflights %>%
	group_by(TailNum) %>%
	summarise(n_dest = n_distinct(Dest)) %>%
	filter(n_dest == 1) %>%
	summarise(nplanes = n())

	# Find the most visited destination for each carrier and save your solution to adv3.
	# Your solution should contain four columns:
	# UniqueCarrier and Dest,
	# n, how often a carrier visited a particular destination,
	# rank, how each destination ranks per carrier. rank should be 1 for every row,
	# as you want to find the most visited destination for each carrier.

	adv3 <- hflights %>%
	group_by(UniqueCarrier, Dest) %>%
	summarise(n = n()) %>%
	mutate(rank = rank(desc(n))) %>%
	filter(rank == 1)

	# Find the carrier that travels to each destination the most: adv4
	# For each destination, find the carrier that travels to that destination the most.
	# Store the result in adv4. Again, your solution should contain 4 columns:
	# Dest, UniqueCarrier, n and rank.

	adv4 <- hflights %>%
	group_by(Dest, UniqueCarrier) %>%
	summarise(n = n()) %>%
	mutate(rank = rank(desc(n))) %>%
	filter(rank == 1)

view raw group_by_miscellaneous.R hosted with ❤ by GitHub

December 22, 2015

Data Manipulation in R with dplyr – Part 2

Note that this post is in continuation with Part 1 of this series of posts on data manipulation with dplyr in R. The code in this post carries forward from the variables / objects defined in Part 1.

In the previous post, I talked about how dplyr provides a grammar of sorts to manipulate data, and consists of 5 verbs to do so:

The 5 verbs of dplyr
select – removes columns from a dataset
filter – removes rows from a dataset
arrange – reorders rows in a dataset
mutate – uses the data to build new columns and values
summarize – calculates summary statistics

I went on to discuss examples using select() and mutate(). Let’s now talk about filter(). R comes with a set of logical operators that you can use inside filter(). These operators are:
x < y, TRUE if x is less than y
x <= y, TRUE if x is less than or equal to y
x == y, TRUE if x equals y
x != y, TRUE if x does not equal y
x >= y, TRUE if x is greater than or equal to y
x > y, TRUE if x is greater than y
x %in% c(a, b, c), TRUE if x is in the vector c(a, b, c)

The following call, for example, filters df such that only the observations where the variable a is greater than the variable b:
filter(df, a > b)

	# Print out all flights in hflights that traveled 3000 or more miles
	filter(hflights, Distance > 3000)

	# All flights flown by one of JetBlue, Southwest, or Delta
	filter(hflights, UniqueCarrier %in% c('JetBlue', 'Southwest', 'Delta'))

	# All flights where taxiing took longer than flying
	filter(hflights, TaxiIn + TaxiOut > AirTime)

view raw verbs05.r hosted with ❤ by GitHub

Combining tests using boolean operators
R also comes with a set of boolean operators that you can use to combine multiple logical tests into a single test. These include & (and), | (or), and ! (not). Instead of using the & operator, you can also pass several logical tests to filter(), separated by commas. The following calls equivalent:

filter(df, a > b & c > d)
filter(df, a > b, c > d)

The is.na() will also come in handy very often. This expression, for example, keeps the observations in df for which the variable x is not NA:

filter(df, !is.na(x))

	# Combining tests using boolean operators

	# All flights that departed before 5am or arrived after 10pm
	filter(hflights, DepTime < 500 \| ArrTime > 2200 )

	# All flights that departed late but arrived ahead of schedule
	filter(hflights, DepDelay > 0 & ArrDelay < 0)

	# All cancelled weekend flights
	filter(hflights, DayOfWeek %in% c(6,7) & Cancelled == 1)

	# All flights that were cancelled after being delayed
	filter(hflights, Cancelled == 1, DepDelay > 0)

view raw verbs06.r hosted with ❤ by GitHub

A recap on select(), mutate() and filter():

	# Summarizing Exercise
	# Select the flights that had JFK as their destination: c1
	c1 <- filter(hflights, Dest == 'JFK')

	# Combine the Year, Month and DayofMonth variables to create a Date column: c2
	c2 <- mutate(c1, Date = paste(Year, Month, DayofMonth, sep = "-"))

	# Print out a selection of columns of c2
	select(c2, Date, DepTime, ArrTime, TailNum)

	# How many weekend flights flew a distance of more than 1000 miles
	# but had a total taxiing time below 15 minutes?
	nrow(filter(hflights, DayOfWeek %in% c(6,7), Distance > 1000, TaxiIn + TaxiOut < 15))

view raw verbs07.r hosted with ❤ by GitHub

Arranging Data
arrange() can be used to rearrange rows according to any type of data. If you pass arrange() a character variable, R will rearrange the rows in alphabetical order according to values of the variable. If you pass a factor variable, R will rearrange the rows according to the order of the levels in your factor (running levels() on the variable reveals this order).

By default, arrange() arranges the rows from smallest to largest. Rows with the smallest value of the variable will appear at the top of the data set. You can reverse this behaviour with the desc() function. arrange() will reorder the rows from largest to smallest values of a variable if you wrap the variable name in desc() before passing it to arrange()

	# Definition of dtc
	dtc <- filter(hflights, Cancelled == 1, !is.na(DepDelay))

	# Arrange dtc by departure delays
	arrange(dtc, DepDelay)

	# Arrange dtc so that cancellation reasons are grouped
	arrange(dtc, CancellationCode)

	# Arrange dtc according to carrier and departure delays
	arrange(dtc, UniqueCarrier, DepDelay)

	# Arrange according to carrier and decreasing departure delays
	arrange(hflights, UniqueCarrier, desc(DepDelay))

	# Arrange flights by total delay (normal order).
	arrange(hflights, DepDelay + ArrDelay)

	# Keep flights leaving to DFW before 8am and arrange according to decreasing AirTime
	arrange(filter(hflights, Dest == 'DFW', DepTime < 800), desc(AirTime))

view raw verbs08.r hosted with ❤ by GitHub

Summarizing Data

summarise(), the last of the 5 verbs, follows the same syntax as mutate(), but the resulting dataset consists of a single row instead of an entire new column in the case of mutate().

In contrast to the four other data manipulation functions, summarise() does not return an altered copy of the dataset it is summarizing; instead, it builds a new dataset that contains only the summarizing statistics.

Note: summarise() and summarize() both work the same!

You can use any function you like in summarise(), so long as the function can take a vector of data and return a single number. R contains many aggregating functions. Here are some of the most useful:

min(x) – minimum value of vector x.
max(x) – maximum value of vector x.
mean(x) – mean value of vector x.
median(x) – median value of vector x.
quantile(x, p) – pth quantile of vector x.
sd(x) – standard deviation of vector x.
var(x) – variance of vector x.
IQR(x) – Inter Quartile Range (IQR) of vector x.
diff(range(x)) – total range of vector x.

	# Print out a summary with variables min_dist and max_dist
	summarize(hflights, min_dist = min(Distance), max_dist = max(Distance))

	# Print out a summary with variable max_div
	summarize(filter(hflights, Diverted == 1), max_div = max(Distance))

	# Remove rows that have NA ArrDelay: temp1
	temp1 <- filter(hflights, !is.na(ArrDelay))

	# Generate summary about ArrDelay column of temp1
	summarise(temp1, earliest = min(ArrDelay), average = mean(ArrDelay),
	latest = max(ArrDelay), sd = sd(ArrDelay))

	# Keep rows that have no NA TaxiIn and no NA TaxiOut: temp2
	temp2 <- filter(hflights, !is.na(TaxiIn), !is.na(TaxiOut))

	# Print the maximum taxiing difference of temp2 with summarise()
	summarise(temp2, max_taxi_diff = max(abs(TaxiIn – TaxiOut)))

view raw verbs09.r hosted with ❤ by GitHub

dplyr provides several helpful aggregate functions of its own, in addition to the ones that are already defined in R. These include:

first(x) – The first element of vector x.
last(x) – The last element of vector x.
nth(x, n) – The nth element of vector x.
n() – The number of rows in the data.frame or group of observations that summarise() describes.
n_distinct(x) – The number of unique values in vector x

	# Generate summarizing statistics for hflights
	summarise(hflights, n_obs = n(), n_carrier = n_distinct(UniqueCarrier),
	n_dest = n_distinct(Dest), dest100 = nth(Dest, 100))

	# Filter hflights to keep all American Airline flights: aa
	aa <- filter(hflights, UniqueCarrier == "American")

	# Generate summarizing statistics for aa
	summarise(aa, n_flights = n(), n_canc = sum(Cancelled),
	p_canc = 100*(n_canc/n_flights), avg_delay = mean(ArrDelay, na.rm = TRUE))

view raw verbs10.r hosted with ❤ by GitHub

This would be it for Part-2 of this series of posts on data manipulation with dplyr. Part 3 would focus on the pipe operator, Group_by and working with databases.

December 18, 2015

Data Manipulation in R with dplyr – Part 1

dplyr is one of the packages in R that makes R so loved by data scientists. It has three main goals:

Identify the most important data manipulation tools needed for data analysis and make them easy to use in R.
Provide blazing fast performance for in-memory data by writing key pieces of code in C++.
Use the same code interface to work with data no matter where it’s stored, whether in a data frame, a data table or database.

Introduction to the dplyr package and the tbl class
This post is mostly about code. If you’re interested in learning dplyr I recommend you type in the commands line by line on the R console to see first hand what’s happening.

	# INTRODUCTION TO dplyr AND tbls
	# Load the dplyr package
	library(dplyr)

	# Load the hflights package
	library(hflights)

	# Call both head() and summary() on hflights
	head(hflights)
	summary(hflights)


	# Convert the hflights data.frame into a hflights tbl
	hflights <- tbl_df(hflights)

	# Display the hflights tbl
	hflights

	# Create the object carriers, containing only the UniqueCarrier variable of hflights
	carriers <- hflights$UniqueCarrier


	# Use lut to translate the UniqueCarrier column of hflights and before doing so
	# glimpse hflights to see the UniqueCarrier variablle
	glimpse(hflights)

	lut <- c("AA" = "American", "AS" = "Alaska", "B6" = "JetBlue", "CO" = "Continental",
	"DL" = "Delta", "OO" = "SkyWest", "UA" = "United", "US" = "US_Airways",
	"WN" = "Southwest", "EV" = "Atlantic_Southeast", "F9" = "Frontier",
	"FL" = "AirTran", "MQ" = "American_Eagle", "XE" = "ExpressJet", "YV" = "Mesa")
	hflights$UniqueCarrier <- lut[hflights$UniqueCarrier]
	# Now glimpse hflights to see the change in the UniqueCarrier variable
	glimpse(hflights)

	# Fill up empty entries of CancellationCode with 'E'
	# To do so, first index the empty entries in CancellationCode
	cancellationEmpty <- hflights$CancellationCode == ""
	# Assign 'E' to the empty entries
	hflights$CancellationCode[cancellationEmpty] <- 'E'

	# Use a new lookup table to create a vector of code labels. Assign the vector to the CancellationCode column of hflights
	lut = c('A' = 'carrier', 'B' = 'weather', 'C' = 'FFA', 'D' = 'security', 'E' = 'not cancelled')
	hflights$CancellationCode <- lut[hflights$CancellationCode]

	# Inspect the resulting raw values of your variables
	glimpse(hflights)

view raw introduction.R hosted with ❤ by GitHub

Select and mutate
dplyr provides grammar for data manipulation apart from providing data structure. The grammar is built around 5 functions (also referred to as verbs) that do the basic tasks of data manipulation.

The 5 verbs of dplyr
select – removes columns from a dataset
filter – removes rows from a dataset
arrange – reorders rows in a dataset
mutate – uses the data to build new columns and values
summarize – calculates summary statistics

dplyr functions do not change the dataset. They return a new copy of the dataset to use.

To answer the simple question whether flight delays tend to shrink or grow during a flight, we can safely discard a lot of the variables of each flight. To select only the ones that matter, we can use select()

	hflights[c('ActualElapsedTime','ArrDelay','DepDelay')]
	# Equivalently, using dplyr:
	select(hflights, ActualElapsedTime, ArrDelay, DepDelay)

	# Print out a tbl with the four columns of hflights related to delay
	select(hflights, ActualElapsedTime, AirTime, ArrDelay, DepDelay)

	# Print out hflights, nothing has changed!
	hflights

	# Print out the columns Origin up to Cancelled of hflights
	select(hflights, Origin:Cancelled)

	# Find the most concise way to select: columns Year up to and
	# including DayOfWeek, columns ArrDelay up to and including Diverted
	# Answer to last question: be concise!
	# You may want to examine the order of hflight's column names before you
	# begin with names()
	names(hflights)
	select(hflights, -(DepTime:AirTime))

view raw verbs01.R hosted with ❤ by GitHub

dplyr comes with a set of helper functions that can help you select variables. These functions find groups of variables to select, based on their names. Each of these works only when used inside of select()

starts_with(“X”): every name that starts with “X”
ends_with(“X”): every name that ends with “X”
contains(“X”): every name that contains “X”
matches(“X”): every name that matches “X”, where “X” can be a regular expression
num_range(“x”, 1:5): the variables named x01, x02, x03, x04 and x05
one_of(x): every name that appears in x, which should be a character vector

	# Helper functions used with dplyr

	# Print out a tbl containing just ArrDelay and DepDelay
	select(hflights, ArrDelay, DepDelay)
	# Use a combination of helper functions and variable names to print out
	# only the UniqueCarrier, FlightNum, TailNum, Cancelled, and CancellationCode
	# columns of hflights
	select(hflights, UniqueCarrier, FlightNum, contains("Tail"), contains("Cancel"))

	# Find the most concise way to return the following columns with select and its
	# helper functions: DepTime, ArrTime, ActualElapsedTime, AirTime, ArrDelay,
	# DepDelay. Use only helper functions
	select(hflights, ends_with("Time"), ends_with("Delay"))

view raw verbs02.R hosted with ❤ by GitHub

In order to appreciate the usefulness of dplyr, here are some comparisons between base R and dplyr

	# Some comparisons to basic R
	# both hflights and dplyr are available

	ex1r <- hflights[c("TaxiIn","TaxiOut","Distance")]
	ex1d <- select(hflights, TaxiIn, TaxiOut, Distance)

	ex2r <- hflights[c("Year","Month","DayOfWeek","DepTime","ArrTime")]
	ex2d <- select(hflights, Year:ArrTime, -DayofMonth)

	ex3r <- hflights[c("TailNum","TaxiIn","TaxiOut")]
	ex3d <- select(hflights, TailNum, contains("Taxi"))

view raw comparisons01.R hosted with ❤ by GitHub

mutate() is the second of the five data manipulation functions. mutate() creates new columns which are added to a copy of the dataset.

	# Add the new variable ActualGroundTime to a copy of hflights and save the result as g1.
	g1 <- mutate(hflights, ActualGroundTime = ActualElapsedTime – AirTime)

	# Add the new variable GroundTime to a g1. Save the result as g2.
	g2 <- mutate(g1, GroundTime = TaxiIn + TaxiOut)

	# Add the new variable AverageSpeed to g2. Save the result as g3.
	g3 <- mutate(g2, AverageSpeed = Distance / AirTime * 60)

	# Print out g3
	g3

view raw verbs03.r hosted with ❤ by GitHub

So far we have added variables to hflights one at a time, but we can also use mutate() to add multiple variables at once.

	# Add a second variable loss_percent to the dataset: m1
	m1 <- mutate(hflights, loss = ArrDelay – DepDelay, loss_percent = ((ArrDelay – DepDelay)/DepDelay)*100)

	# mutate() allows you to use a new variable while creating a next variable in the same call
	# Copy and adapt the previous command to reduce redendancy: m2
	m2 <- mutate(hflights, loss = ArrDelay – DepDelay, loss_percent = (loss/DepDelay) * 100 )

	# Add the three variables as described in the third instruction: m3
	m3 <- mutate(hflights, TotalTaxi = TaxiIn + TaxiOut, ActualGroundTime = ActualElapsedTime – AirTime, Diff = TotalTaxi – ActualGroundTime)

view raw verbs04.r hosted with ❤ by GitHub

December 17, 2015

Generating Permutation Matrices in Octave / Matlab

I have been doing Gilbert Strang’s linear algebra assignments, some of which require you to write short scripts in MatLab, though I use GNU Octave (which is kind of like a free MatLab). I was trying out this problem:

To solve this quickly, it would have been nice to have a function that would give a list of permutation matrices for every n-sized square matrix, but there was none in Octave, so I wrote a function permMatrices which creates a list of permutation matrices for a square matrix of size n.

	% function to generate permutation matrices given the size of the desired permutation matrices
	function x = permMatrices(n)
	x = zeros(n,n,factorial(n));
	permutations = perms(1:n);
	for i = 1:size(x,3)
	x(:,:,i) = eye(n)(permutations(i,:),:);
	end
	endfunction

view raw permMatrices.m hosted with ❤ by GitHub

For example:

The MatLab / Octave code to solve this problem is shown below:

	% Solution for part (a)
	p = permMatrices(3);
	n = size(p,3); % number of permutation matrices
	v = zeros(n,1); % vector of zeros with dimension equalling number of permutation matrices
	% check for permutation matrices other than identity matrix with 3rd power equalling identity matrix
	for i = 1:n
	if p(:,:,i)^3 == eye(3)
	v(i,1) = 1;
	end
	end
	v(1,1) = 0; % exclude identity matrix
	ans1 = p(:,:,v == 1)

	% Solution for part (b)
	P = permMatrices(4);
	m = size(P,3); % number of permutation matrices
	t = zeros(m,1); % vector of zeros with dimension equalling number of permutation matrices
	% check for permutation matrices with 4th power equalling identity matrix
	for i = 1:m
	if P(:,:,i)^4 == eye(4)
	t(i,1) = 1;
	end
	end
	% print the permutation matrices
	ans2 = P(:,:,t == 0)

view raw Section2_7_13.m hosted with ❤ by GitHub

Output:

December 14, 2015

Sherlock and the Beast – HackerRank

I found myself stuck on this problem recently. I must confess, I lost a couple of hours trying to get to figure the logic for this one. Here’s the problem:

I’ve written 2 functions to solve this problem. The first one I used for smaller N, say N < 30 and the second one for N > 30. The second function is elegant, and it relies on the mathematical property that if a number N is not divisible by 3, it could either leave a remainder 1 or 2.

If it leaves a remainder 2, then subtracting 5 once would make the number divisible by 3. If it leaves a remainder 1, then subtracting 5 twice would make the number divisible by 3.

We subtract 5 from N iteratively and attempt to divide N into 2 parts, one divisible by 3 and the other divisible by 5. We want the part that is divisible by 3 to be the larger part, so that the associated Decent Number is the largest possible. This explanation might seem obtuse, but if you get pen down on paper, you’ll understand what I mean.

Solution

December 12, 2015

Supplementary Material to Andrew Ng’s Machine Learning MOOC

Although the lecture videos and lecture notes from Andrew Ng‘s Coursera MOOC are sufficient for the online version of the course, if you’re interested in more mathematical stuff or want to be challenged further, you can go through the following notes and problem sets from CS 229, a 10-week course that he teaches at Stanford (which also happens to be the most enrolled course on campus). It’s not hard to end up with a 100% score on his MOOC which is obviously a (much) watered down version of the course he teaches at Stanford, at least in terms of difficulty. If you don’t believe me, just have a go at the problem sets from the links below.

Lecture Notes

Lecture notes 1 (ps) (pdf) Supervised Learning, Discriminative Algorithms
Lecture notes 2 (ps) (pdf) Generative Algorithms
Lecture notes 3 (ps) (pdf) Support Vector Machines
Lecture notes 4 (ps) (pdf) Learning Theory
Lecture notes 5 (ps) (pdf) Regularization and Model Selection
Lecture notes 6 (ps) (pdf) Online Learning and the Perceptron Algorithm. (optional reading)
Lecture notes 7a (ps) (pdf) Unsupervised Learning, k-means clustering.
Lecture notes 7b (ps) (pdf) Mixture of Gaussians
Lecture notes 8 (ps) (pdf) The EM Algorithm
Lecture notes 9 (ps) (pdf) Factor Analysis
Lecture notes 10 (ps) (pdf) Principal Components Analysis
Lecture notes 11 (ps) (pdf) Independent Components Analysis
Lecture notes 12 (ps) (pdf) Reinforcement Learning and Control

Section Notes

Section notes 1 (pdf) Linear Algebra Review and Reference
Section notes 2 (pdf) Probability Theory Review
Files for the Matlab tutorial: sigmoid.m, logistic_grad_ascent.m, matlab_session.m
Section notes 4 (ps) (pdf) Convex Optimization Overview, Part I
Section notes 5 (ps) (pdf) Convex Optimization Overview, Part II
Section notes 6 (ps) (pdf) Hidden Markov Models
Section notes 7 (pdf) The Multivariate Gaussian Distribution
Section notes 8 (pdf) More on Gaussian Distribution
Section notes 9 (pdf) Gaussian Processes

Handouts and Problem Sets

Handout #1: Course Information (HTML) (pdf)
Handout #2: Course Schedule (HTML) (pdf)
Handout #3: Cover Sheet
Handout #4: Practice Midterm 1 Solution: Solution
Handout #5: Practice Midterm 2 Solution: Solution
Problem Set 1 (pdf) Data: q1x.dat, q1y.dat, q2x.dat, q2y.dat Solution: Solution (pdf)
Problem Set 2 (pdf) Data: ps2.zip Solution: Solution (pdf)
Problem Set 3 (pdf) Solution: Solution (pdf)
Problem Set 4 (pdf)

November 25, 2015

Solutions to Machine Learning Programming Assignments

This post contains links to a bunch of code that I have written to complete Andrew Ng’s famous machine learning course which includes several interesting machine learning problems that needed to be solved using the Octave / Matlab programming language. I’m not sure I’d ever be programming in Octave after this course, but learning Octave just so that I could complete this course seemed worth the time and effort. I would usually work on the programming assignments on Sundays and spend several hours coding in Octave, telling myself that I would later replicate the exercises in Python.

If you’ve taken this course and found some of the assignments hard to complete, I think it might not hurt to go check online on how a particular function was implemented. If you end up copying the entire code, it’s probably your loss in the long run. But then John Maynard Keynes once said, ‘In the long run we are all dead‘. Yeah, and we wonder why people call Economics the dismal science!

Most people disregard Coursera’s feeble attempt at reigning in plagiarism by creating an Honor Code, precisely because this so-called code-of-conduct can be easily circumvented. I don’t mind posting solutions to a course’s programming assignments because GitHub is full to the brim with such content. Plus, it’s always good to read others’ code even if you implemented a function correctly. It helps understand the different ways of tackling a given programming problem.

ex1
ex2
ex3
ex4
ex5
ex6
ex7
ex8

Enjoy!

November 24, 2015

Troubleshooting ‘Rattle’ (R library) Installation on Ubuntu

This post pertains to Ubuntu / Debian users only.

rattle is a free graphical interface for data mining with R. I wanted to visualize decision trees and had to install this library.
> install.packages('rattle')
got me the following error message:

configure: error: GTK version 2.8.0 required ERROR: configuration failed for package ‘RGtk2’

This error occurs when attempting to install the RGtk2 package. The install is looking for the header files for GTK. Possibly they are not yet. Luckily the problem can be solved quite easily. Open Terminal (Ctrl + Alt + T) and type in the following commands:

sudo apt-get update wajig install libgtk2.0-dev

Go back and try installing rattle now with the same command as earlier. It should work. It did for me! As you can see below, decision trees are visualized lot better with rattle than if you used just rpart.

November 23, 2015

Spot the Difference — It’s NumPy!

My first brush with NumPy happened over writing a block of code to make a plot using pylab. ⇣

pylab is part of matplotlib (in matplotlib.pylab) and tries to give you a MatLab like environment. matplotlib has a number of dependencies, among them numpy which it imports under the common alias np. scipy is not a dependency of matplotlib.

I had a tuple (of lows and highs of temperature) of lengh 2 with 31 entries in each (the number of days in the month of July), parsed from this text file:

	Boston July Temperatures
	————————-

	Day High Low
	————

	1 91 70
	2 84 69
	3 86 68
	4 84 68
	5 83 70
	6 80 68
	7 86 73
	8 89 71
	9 84 67
	10 83 65
	11 80 66
	12 86 63
	13 90 69
	14 91 72
	15 91 72
	16 88 72
	17 97 76
	18 89 70
	19 74 66
	20 71 64
	21 74 61
	22 84 61
	23 86 66
	24 91 68
	25 83 65
	26 84 66
	27 79 64
	28 72 63
	29 73 64
	30 81 63
	31 73 63

view raw julyTemps.txt hosted with ❤ by GitHub

Given below, are 2 sets of code that do the same thing; one without NumPy and the other with NumPy. They output the following graph using PyLab:

Code without NumPy

	import pylab

	def loadfile():
	inFile = open('julyTemps.txt', 'r')
	high =[]; low = []
	for line in inFile:
	fields = line.split()
	if len(fields) < 3 or not fields[0].isdigit():
	pass
	else:
	high.append(int(fields[1]))
	low.append(int(fields[2]))
	return low, high

	def producePlot(lowTemps, highTemps):
	diffTemps = [highTemps[i] – lowTemps[i] for i in range(len(lowTemps))]
	pylab.title('Day by Day Ranges in Temperature in Boston in July 2012')
	pylab.xlabel('Days')
	pylab.ylabel('Temperature Ranges')
	return pylab.plot(range(1,32),diffTemps)

	producePlot(loadfile()[1], loadfile()[0])

view raw withoutNumPy.py hosted with ❤ by GitHub

Code with NumPy

	import pylab
	import numpy as np

	def loadFile():
	inFile = open('julyTemps.txt')
	high = [];vlow = []
	for line in inFile:
	fields = line.split()
	if len(fields) != 3 or 'Boston' == fields[0] or 'Day' == fields[0]:
	continue
	else:
	high.append(int(fields[1]))
	low.append(int(fields[2]))
	return (low, high)

	def producePlot(lowTemps, highTemps):
	diffTemps = list(np.array(highTemps) – np.array(lowTemps))
	pylab.plot(range(1,32), diffTemps)
	pylab.title('Day by Day Ranges in Temperature in Boston in July 2012')
	pylab.xlabel('Days')
	pylab.ylabel('Temperature Ranges')
	pylab.show()


	(low, high) = loadFile()
	producePlot(low, high)

view raw withNumPy.py hosted with ❤ by GitHub

The difference in code lies in how the variable diffTemps is calculated.

diffTemps = list(np.array(highTemps) - np.array(lowTemps))

seems more readable than

diffTemps = [highTemps[i] - lowTemps[i] for i in range(len(lowTemps))]

Notice how straight forward it is with NumPy. At the core of the NumPy package, is the ndarray object. This encapsulates n-dimensional arrays of homogeneous data types, with many operations being performed in compiled code for performance. element-by-element operations are the “default mode” when an ndarray is involved, but the element-by-element operation is speedily executed by pre-compiled C code.

October 22, 2015

Karatsuba Multiplication Algorithm – Python Code

Motivation for this blog post

I’ve enrolled in Stanford Professor Tim Roughgarden’s Coursera MOOC on the design and analysis of algorithms, and while he covers the theory and intuition behind the algorithms in a surprising amount of detail, we’re left to implement them in a programming language of our choice.

~~And I’m ging to post Python code for all the algorithms covered during the course!~~

The Karatsuba Multiplication Algorithm

Karatsuba’s algorithm reduces the multiplication of two n-digit numbers to at most $n^{\log_23}\approx n^{1.585}$ single-digit multiplications in general (and exactly $n^{\log_23}$ when n is a power of 2). Although the familiar grade school algorithm for multiplying numbers is how we work through multiplication in our day-to-day lives, it’s slower ( $\Theta(n^2)\,\!$ ) in comparison, but only on a computer, of course!

Here’s how the grade school algorithm looks:
(The following slides have been taken from Tim Roughgarden’s notes. They serve as a good illustration. I hope he doesn’t mind my sharing them.)

…and this is how Karatsuba Multiplication works on the same problem:

A More General Treatment

Let $x$ and $y$ be represented as $n$ -digit strings in some base $B$ . For any positive integer $m$ less than $n$ , one can write the two given numbers as

$x = x_1B^m + x_0$
$y = y_1B^m + y_0$ ,

where $x_0$ and $y_0$ are less than $B^m$ . The product is then

$xy = (x_1B^m + x_0)(y_1B^m + y_0)$
$xy = z_2B^{2m} + z_1B^m + z_0$

where

$z_2 = x_1y_1$
$z_1 = x_1y_0 + x_0y_1$
$z_0 = x_0y_0$

These formulae require four multiplications, and were known to Charles Babbage. Karatsuba observed that $xy$ can be computed in only three multiplications, at the cost of a few extra additions. With $z_0$ and $z_2$ as before we can calculate

$z_1 = (x_1 + x_0)(y_1 + y_0) - z_2 - z_0$

which holds since

$z_1 = x_1y_0 + x_0y_1$
$z_1 = (x_1 + x_0)(y_1 + y_0) - x_1y_1 - x_0y_0$

A more efficient implementation of Karatsuba multiplication can be set as $xy = (b^2 + b)x_1y_1 - b(x_1 - x_0)(y_1 - y_0) + (b + 1)x_0y_0$ , where $b = B^m$ .

Example

To compute the product of 12345 and 6789, choose B = 10 and m = 3. Then we decompose the input operands using the resulting base (B^m = 1000), as:

12345 = 12 · 1000 + 345

6789 = 6 · 1000 + 789

Only three multiplications, which operate on smaller integers, are used to compute three partial results:

z₂ = 12 × 6 = 72

z₀ = 345 × 789 = 272205

z₁ = (12 + 345) × (6 + 789) − z₂ − z₀ = 357 × 795 − 72 − 272205 = 283815 − 72 − 272205 = 11538

We get the result by just adding these three partial results, shifted accordingly (and then taking carries into account by decomposing these three inputs in base 1000 like for the input operands):

result = z₂ · B^2m + z₁ · B^m + z₀, i.e.

result = 72 · 1000² + 11538 · 1000 + 272205 = 83810205.

Pseudocode and Python code

	procedure karatsuba(num1, num2)
	if (num1 < 10) or (num2 < 10)
	return num1*num2
	/* calculates the size of the numbers */
	m = max(size_base10(num1), size_base10(num2))
	m2 = m/2
	/* split the digit sequences about the middle */
	high1, low1 = split_at(num1, m2)
	high2, low2 = split_at(num2, m2)
	/* 3 calls made to numbers approximately half the size */
	z0 = karatsuba(low1,low2)
	z1 = karatsuba((low1+high1),(low2+high2))
	z2 = karatsuba(high1,high2)
	return (z210^(2m2))+((z1-z2-z0)*10^(m2))+(z0)

view raw karatsuba_pseudocode.txt hosted with ❤ by GitHub

	def karatsuba(x,y):
	"""Function to multiply 2 numbers in a more efficient manner than the grade school algorithm"""
	if len(str(x)) == 1 or len(str(y)) == 1:
	return x*y
	else:
	n = max(len(str(x)),len(str(y)))
	nby2 = n / 2

	a = x / 10**(nby2)
	b = x % 10**(nby2)
	c = y / 10**(nby2)
	d = y % 10**(nby2)

	ac = karatsuba(a,c)
	bd = karatsuba(b,d)
	ad_plus_bc = karatsuba(a+b,c+d) – ac – bd

	# this little trick, writing n as 2*nby2 takes care of both even and odd n
	prod = ac * 10*(2nby2) + (ad_plus_bc * 10**nby2) + bd

	return prod

view raw karatsuba.py hosted with ❤ by GitHub

October 13, 2015

Category: Technical

Share this:

Share this:

Share this:

Share this:

Share this:

Share this:

Share this:

Share this:

Share this:

Example

Share this: