Anirudh Jayaraman

Category: Technical

  • Deterministic Selection Algorithm Python Code

    Deterministic Selection Algorithm Python Code

    Through this post, I’m sharing Python code implementing the median of medians algorithm, an algorithm that resembles quickselect, differing only in the way in which the pivot is chosen, i.e, deterministically, instead of at random.

    Its best case complexity is O(n) and worst case complexity O(nlog2n)

    I don’t have a formal education in CS, and came across this algorithm while going through Tim Roughgarden’s Coursera MOOC on the design and analysis of algorithms. Check out my implementation in Python.

    def merge_tuple(a,b):
    """ Function to merge two arrays of tuples """
    c = []
    while len(a) != 0 and len(b) != 0:
    if a[0][0] < b[0][0]:
    c.append(a[0])
    a.remove(a[0])
    else:
    c.append(b[0])
    b.remove(b[0])
    if len(a) == 0:
    c += b
    else:
    c += a
    return c
    def mergesort_tuple(x):
    """ Function to sort an array using merge sort algorithm """
    if len(x) == 0 or len(x) == 1:
    return x
    else:
    middle = len(x)/2
    a = mergesort_tuple(x[:middle])
    b = mergesort_tuple(x[middle:])
    return merge_tuple(a,b)
    def lol(x,k):
    """ Function to divide a list into a list of lists of size k each. """
    return [x[i:i+k] for i in range(0,len(x),k)]
    def preprocess(x):
    """ Function to assign an index to each element of a list of integers, outputting a list of tuples"""
    return zip(x,range(len(x)))
    def partition(x, pivot_index = 0):
    """ Function to partition an unsorted array around a pivot"""
    i = 0
    if pivot_index !=0: x[0],x[pivot_index] = x[pivot_index],x[0]
    for j in range(len(x)-1):
    if x[j+1] < x[0]:
    x[j+1],x[i+1] = x[i+1],x[j+1]
    i += 1
    x[0],x[i] = x[i],x[0]
    return x,i
    def ChoosePivot(x):
    """ Function to choose pivot element of an unsorted array using 'Median of Medians' method. """
    if len(x) <= 5:
    return mergesort_tuple(x)[middle_index(x)]
    else:
    lst = lol(x,5)
    lst = [mergesort_tuple(el) for el in lst]
    C = [el[middle_index(el)] for el in lst]
    return ChoosePivot(C)
    def DSelect(x,k):
    """ Function to """
    if len(x) == 1:
    return x[0]
    else:
    xpart = partition(x,ChoosePivot(preprocess(x))[1])
    x = xpart[0] # partitioned array
    j = xpart[1] # pivot index
    if j == k:
    return x[j]
    elif j > k:
    return DSelect(x[:j],k)
    else:
    k = k – j – 1
    return DSelect(x[(j+1):], k)
    arr = range(100,0,-1)
    print DSelect(arr,50)
    %timeit DSelect(arr,50)
    view raw DSelect.py hosted with ❤ by GitHub

    I get the following output:

    51
100 loops, best of 3: 2.38 ms per loop

    Note that on the same input, quickselect is faster, giving us:

    1000 loops, best of 3: 254 µs per loop
  • scikit-learn Linear Regression Example

    scikit-learn Linear Regression Example

    Here’s a quick example case for implementing one of the simplest of learning algorithms in any machine learning toolbox – Linear Regression. You can download the IPython / Jupyter notebook here so as to play around with the code and try things out yourself.

    I’m doing a series of posts on scikit-learn. Its documentation is vast, so unless you’re willing to search for a needle in a haystack, you’re better off NOT jumping into the documentation right away. Instead, knowing chunks of code that do the job might help.

    Loading
    Sorry, something went wrong. Reload?
    Sorry, we cannot display this file.
    Sorry, this file is invalid so it cannot be displayed.
    view raw linear_regression.ipynb hosted with ❤ by GitHub
  • Sharing IPython / Jupyter Notebooks via WordPress

    Sharing IPython / Jupyter Notebooks via WordPress

    In order to share (a static version of) your IPython / Jupyter notebook on your WordPress site, follow three straightforward steps.

    Step 1: Let’s say your Jupyter Notebook looks like this:

    blog_item_20160718_01

    Open this notebook in a text editor and copy the content which may look like so:

    blog_item_20160718_02

    Step 2: Ctrl + A and Ctrl + C this content. Then Ctrl + V this to a GitHub Gist that you should create, like so:

    blog_item_20160718_03

    Step 3: Now simply Create public gist and embed the gist like you always embed gists on WordPress, viz., go to the HTML editor and add like so:

    blog_item_20160718_04

    I followed the exact steps that I’ve mentioned above to get the following result:

    Loading
    Sorry, something went wrong. Reload?
    Sorry, we cannot display this file.
    Sorry, this file is invalid so it cannot be displayed.
    view raw knn.ipynb hosted with ❤ by GitHub

     

  • Randomized Selection Algorithm (Quickselect) – Python Code

    Randomized Selection Algorithm (Quickselect) – Python Code

    Find the kth smallest element in an array without sorting.

    That’s basically what this algorithm does. It piggybacks on the partition subroutine from the Quick Sort. If you don’t know what that is, you can check out more about the Quick Sort algorithm here and here, and understand the usefulness of partitioning an unsorted array around a pivot.

    Selecting_quickselect_frames
    Animated visualization of the randomized selection algorithm selecting the 22nd
     smallest value

    Python Implementation

    from random import randrange
    def partition(x, pivot_index = 0):
    i = 0
    if pivot_index !=0: x[0],x[pivot_index] = x[pivot_index],x[0]
    for j in range(len(x)-1):
    if x[j+1] < x[0]:
    x[j+1],x[i+1] = x[i+1],x[j+1]
    i += 1
    x[0],x[i] = x[i],x[0]
    return x,i
    def RSelect(x,k):
    if len(x) == 1:
    return x[0]
    else:
    xpart = partition(x,randrange(len(x)))
    x = xpart[0] # partitioned array
    j = xpart[1] # pivot index
    if j == k:
    return x[j]
    elif j > k:
    return RSelect(x[:j],k)
    else:
    k = k – j – 1
    return RSelect(x[(j+1):], k)
    x = [3,1,8,4,7,9]
    for i in range(len(x)):
    print RSelect(x,i),
    view raw RSelect.py hosted with ❤ by GitHub

     

    Related Posts
    Quick Sort Python Code
    Computing Work Done (Total Pivot Comparisons) by Quick Sort

  • Computing Work Done (Total Pivot Comparisons) by Quick Sort

    A key aspect of the Quick Sort algorithm is how the pivot element is chosen. In my earlier post on the Python code for Quick Sort, my implementation takes the first element of the unsorted array as the pivot element.

    However with some mathematical analysis it can be seen that such an implementation is O(n2) in complexity while if a pivot is randomly chosen, the Quick Sort algorithm is O(nlog2n).

    To witness this in action, one can measure the work done by the algorithm comparing two cases, one with a randomized pivot choice – and one with a fixed pivot choice, say the first element of the array (or the last element of the array).

    Implementation

    A decent proxy for the amount of work done by the algorithm would be the number of pivot comparisons. These comparisons needn’t be computed one-by-one, rather when there is a recursive call on a subarray of length m, you should simply add m−1 to your running total of comparisons.

    3 Cases

    To put things in perspective, let’s look at 3 cases. (This is basically straight out of a homework assignment from Tim Roughgarden’s course on the Design and Analysis of Algorithms).
    Case I with the pivot being the first element.
    Case II with the pivot being the last element.
    Case III using the “median-of-three” pivot rule. The primary motivation behind this rule is to do a little bit of extra work to get much better performance on input arrays that are nearly sorted or reverse sorted.

    Median-of-Three Pivot Rule

    Consider the first, middle, and final elements of the given array. (If the array has odd length it should be clear what the “middle” element is; for an array with even length 2k, use the kth element as the “middle” element. So for the array 4 5 6 7, the “middle” element is the second one —- 5 and not 6! Identify which of these three elements is the median (i.e., the one whose value is in between the other two), and use this as your pivot.

    Python Code

    This file contains all of the integers between 1 and 10,000 (inclusive, with no repeats) in unsorted order. The integer in the ith row of the file gives you the ith entry of an input array. I downloaded this file and named it QuickSort_List.txt

    You can run the code below and see for yourself that the number of comparisons for Case III are 138,382 compared to 162,085 and 164,123 for Case I and Case II respectively. You can play around with the code in an IPython / Jupyter notebook here.

    #!/usr/bin/env
    # Case I
    # First element of the unsorted array is chosen as pivot element for sorting using Quick Sort
    def countComparisonsWithFirst(x):
    """ Counts number of comparisons while using Quick Sort with first element of unsorted array as pivot """
    global count_pivot_first
    if len(x) == 1 or len(x) == 0:
    return x
    else:
    count_pivot_first += len(x)-1
    i = 0
    for j in range(len(x)-1):
    if x[j+1] < x[0]:
    x[j+1],x[i+1] = x[i+1], x[j+1]
    i += 1
    x[0],x[i] = x[i],x[0]
    first_part = countComparisonsWithFirst(x[:i])
    second_part = countComparisonsWithFirst(x[i+1:])
    first_part.append(x[i])
    return first_part + second_part
    # Case II
    # Last element of the unsorted array is chosen as pivot element for sorting using Quick Sort
    def countComparisonsWithLast(x):
    """ Counts number of comparisons while using Quick Sort with last element of unsorted array as pivot """
    global count_pivot_last
    if len(x) == 1 or len(x) == 0:
    return x
    else:
    count_pivot_last += len(x)-1
    x[0],x[-1] = x[-1],x[0]
    i = 0
    for j in range(len(x)-1):
    if x[j+1] < x[0]:
    x[j+1],x[i+1] = x[i+1], x[j+1]
    i += 1
    x[0],x[i] = x[i],x[0]
    first_part = countComparisonsWithLast(x[:i])
    second_part = countComparisonsWithLast(x[i+1:])
    first_part.append(x[i])
    return first_part + second_part
    # Case III
    # Median-of-three method used to choose pivot element for sorting using Quick Sort
    def middle_index(x):
    """ Returns the index of the middle element of an array """
    if len(x) % 2 == 0:
    middle_index = len(x)/2 – 1
    else:
    middle_index = len(x)/2
    return middle_index
    def median_index(x,i,j,k):
    """ Returns the median index of three when passed an array and indices of any 3 elements of that array """
    if (x[i]-x[j])*(x[i]-x[k]) < 0:
    return i
    elif (x[j]-x[i])*(x[j]-x[k]) < 0:
    return j
    else:
    return k
    def countComparisonsMedianOfThree(x):
    """ Counts number of comparisons while using Quick Sort with median-of-three element is chosen as pivot """
    global count_pivot_median
    if len(x) == 1 or len(x) == 0:
    return x
    else:
    count_pivot_median += len(x)-1
    k = median_index(x, 0, middle_index(x), -1)
    if k != 0: x[0], x[k] = x[k], x[0]
    i = 0
    for j in range(len(x)-1):
    if x[j+1] < x[0]:
    x[j+1],x[i+1] = x[i+1], x[j+1]
    i += 1
    x[0],x[i] = x[i],x[0]
    first_part = countComparisonsMedianOfThree(x[:i])
    second_part = countComparisonsMedianOfThree(x[i+1:])
    first_part.append(x[i])
    return first_part + second_part
    #####################################################################
    # initializing counts
    count_pivot_first = 0; count_pivot_last = 0; count_pivot_median = 0
    #####################################################################
    # Cast I
    # Read the contents of the file into a Python list
    NUMLIST_FILENAME = "QuickSort_List.txt"
    inFile = open(NUMLIST_FILENAME, 'r')
    with inFile as f: numList = [int(integers.strip()) for integers in f.readlines()]
    # call functions to count comparisons
    countComparisonsWithFirst(numList)
    #####################################################################
    # Read the contents of the file into a Python list
    NUMLIST_FILENAME = "QuickSort_List.txt"
    inFile = open(NUMLIST_FILENAME, 'r')
    with inFile as f: numList = [int(integers.strip()) for integers in f.readlines()]
    # call functions to count comparisons
    countComparisonsWithLast(numList)
    #####################################################################
    # Read the contents of the file into a Python list
    NUMLIST_FILENAME = "QuickSort_List.txt"
    inFile = open(NUMLIST_FILENAME, 'r')
    with inFile as f: numList = [int(integers.strip()) for integers in f.readlines()]
    # call functions to count comparisons
    countComparisonsMedianOfThree(numList)
    #####################################################################
    print count_pivot_first, count_pivot_last, count_pivot_median
    view raw countComparisons.py hosted with ❤ by GitHub
  • Quick Sort Python Code

    Quick Sort Python Code

    Sorting_quicksort_anim

    Yet another post for the crawlers to better index my site for algorithms and as a repository for Python code. The quick sort algorithm is well explained in the topmost Google search result for ‘Quick Sort Python Code’, but the code is unnecessarily convoluted. Instead, go with the code below.

    In it, I assume the pivot to be the first element. You can easily add a function to  randomize selection of the pivot. Choosing a random pivot minimizes the chance that you will encounter worst-case O(n2) performance. Always choosing first or last would cause worst-case performance for nearly-sorted or nearly-reverse-sorted data.

    def quicksort(x):
    if len(x) == 1 or len(x) == 0:
    return x
    else:
    pivot = x[0]
    i = 0
    for j in range(len(x)-1):
    if x[j+1] < pivot:
    x[j+1],x[i+1] = x[i+1], x[j+1]
    i += 1
    x[0],x[i] = x[i],x[0]
    first_part = quicksort(x[:i])
    second_part = quicksort(x[i+1:])
    first_part.append(x[i])
    return first_part + second_part
    alist = [54,26,93,17,77,31,44,55,20]
    quicksort(alist)
    print(alist)
    view raw quicksort.py hosted with ❤ by GitHub

    Also read:
    Computing Work Done (Total Pivot Comparisons) by Quick Sort
    Karatsuba Multiplication Algorithm – Python Code
    Merge Sort

  • Detecting Structural Breaks in China’s FX Regime

    Detecting Structural Breaks in China’s FX Regime

    Edit: This post is in its infancy. Work is still ongoing as far as deriving insight from the data is concerned. More content and economic insight is expected to be added to this post as and when progress is made in that direction.

    This is an attempt to detect structural breaks in China’s FX regime using Frenkel Wei regression methodology (this was later improved by Perron and Bai). I came up with the motivation to check for these structural breaks while attending a guest lecture on FX regimes by Dr. Ajay Shah delivered at IGIDR. This is work that I and two other classmates are working on as a term paper project under the supervision of Dr. Rajeswari Sengupta.

    The code below can be replicated and run as is, to get same results.

    ## if fxregime or strucchange package is absent from installed packages, download it and load it
    if(!require('fxregime')){
    install.packages("fxregime")
    }
    if(!require('strucchange')){
    install.packages("strucchange")
    }
    ## load packages
    library("fxregime")
    library('strucchange')
    # load the necessary data related to exchange rates – 'FXRatesCHF'
    # this dataset treats CHF as unit currency
    data("FXRatesCHF", package = "fxregime")
    ## compute returns for CNY (and explanatory currencies)
    ## since China abolished fixed USD regime
    cny <- fxreturns("CNY", frequency = "daily",
    start = as.Date("2005-07-25"), end = as.Date("2010-02-12"),
    other = c("USD", "JPY", "EUR", "GBP"))
    ## compute all segmented regression with minimal segment size of
    ## h = 100 and maximal number of breaks = 10
    regx <- fxregimes(CNY ~ USD + JPY + EUR + GBP,
    data = cny, h = 100, breaks = 10, ic = "BIC")
    ## Print summary of regression results
    summary(regx)
    ## minimum BIC is attained for 2-segment (1-break) model
    plot(regx)
    round(coef(regx), digits = 3)
    sqrt(coef(regx)[, "(Variance)"])
    ## inspect associated confidence intervals
    cit <- confint(regx, level = 0.9)
    cit
    breakdates(cit)
    ## plot LM statistics along with confidence interval
    flm <- fxlm(CNY ~ USD + JPY + EUR + GBP, data = cny)
    scus <- gefp(flm, fit = NULL)
    plot(scus, functional = supLM(0.1))
    ## add lines related to breaks to your plot
    lines(cit)
    view raw FrankelWei_2005_2010.r hosted with ❤ by GitHub

    As can be seen in the figure below, the structural breaks correspond to the vertical bars. We are still working on understanding the motivations of China’s central bank in varying the degree of the managed float exchange rate.

    strucchange_china_2006_2010

    EDIT (May 16, 2016):

    The code above uses data provided by the package itself. If you wished to replicate this analysis on data after 2010, you will have to use your own data. We used Quandl, which lets you get 10 premium datasets for free. An API key (for only 10 calls on premium datasets) is provided if you register there. Foreign exchange rate data (2000 onward till date) apparently, is premium data. You can find these here.

    Here are the (partial) results and code to work the same methodology on the data from 2010 to 2016:

    20102016

    ## if fxregime is absent from installed packages, download it and load it
    if(!require('fxregime')){
    install.packages("fxregime")
    }
    ## load package
    library("fxregime")
    # load the necessary data related to exchange rates – 'FXRatesCHF'
    # this dataset treats CHF as unit currency
    # install / load Quandl
    if(!require('Quandl')){
    install.packages("Quandl")
    }
    library(Quandl)
    # Extract and load currency data series with respect to CHF from Quandl
    # Extract data series from Quandl. Each Quandl user will have unique api_key
    # upon signing up. The freemium version allows access up to 10 fx rate data sets
    # USDCHF <- Quandl("CUR/CHF", api_key="p2GsFxccPGFSw7n1-NF9")
    # write.csv(USDCHF, file = "USDCHF.csv")
    # USDCNY <- Quandl("CUR/CNY", api_key="p2GsFxccPGFSw7n1-NF9")
    # write.csv(USDCNY, file = "USDCNY.csv")
    # USDJPY <- Quandl("CUR/JPY", api_key="p2GsFxccPGFSw7n1-NF9")
    # write.csv(USDJPY, file = "USDJPY.csv")
    # USDEUR <- Quandl("CUR/EUR", api_key="p2GsFxccPGFSw7n1-NF9")
    # write.csv(USDEUR, file = "USDEUR.csv")
    # USDGBP <- Quandl("CUR/GBP", api_key="p2GsFxccPGFSw7n1-NF9")
    # write.csv(USDGBP, file = "USDGBP.csv")
    # load the data sets into R
    USDCHF <- read.csv("G:/China's Economic Woes/USDCHF.csv")
    USDCHF <- USDCHF[,2:3]
    USDCNY <- read.csv("G:/China's Economic Woes/USDCNY.csv")
    USDCNY <- USDCNY[,2:3]
    USDEUR <- read.csv("G:/China's Economic Woes/USDEUR.csv")
    USDEUR <- USDEUR[,2:3]
    USDGBP <- read.csv("G:/China's Economic Woes/USDGBP.csv")
    USDGBP <- USDGBP[,2:3]
    USDJPY <- read.csv("G:/China's Economic Woes/USDJPY.csv")
    USDJPY <- USDJPY[,2:3]
    start = 1 # corresponds to 2016-05-12
    end = 2272 # corresponds to 2010-02-12
    dates <- as.Date(USDCHF[start:end,1])
    USD <- 1/USDCHF[start:end,2]
    CNY <- USDCNY[start:end,2]/USD
    JPY <- USDJPY[start:end,2]/USD
    EUR <- USDEUR[start:end,2]/USD
    GBP <- USDGBP[start:end,2]/USD
    # reverse the order of the vectors to reflect dates from 2005 – 2010 instead of
    # the other way around
    USD <- USD[length(USD):1]
    CNY <- CNY[length(CNY):1]
    JPY <- JPY[length(JPY):1]
    EUR <- EUR[length(EUR):1]
    GBP <- GBP[length(GBP):1]
    dates <- dates[length(dates):1]
    df <- data.frame(CNY, USD, JPY, EUR, GBP)
    df$weekday <- weekdays(dates)
    row.names(df) <- dates
    df <- subset(df, weekday != 'Sunday')
    df <- subset(df, weekday != 'Saturday')
    df <- df[,1:5]
    zoo_df <- as.zoo(df)
    # Code to replicate analysis
    cny_rep <- fxreturns("CNY", data = zoo_df, frequency = "daily",
    other = c("USD", "JPY", "EUR", "GBP"))
    time(cny_rep) <- as.Date(row.names(df)[2:1627])
    regx_rep <- fxregimes(CNY ~ USD + JPY + EUR + GBP,
    data = cny_rep, h = 100, breaks = 10, ic = "BIC")
    summary(regx_rep)
    ## minimum BIC is attained for 2-segment (5-break) model
    plot(regx_rep)
    round(coef(regx_rep), digits = 3)
    sqrt(coef(regx_rep)[, "(Variance)"])
    ## inspect associated confidence intervals
    cit_rep <- confint(regx_rep, level = 0.9)
    breakdates(cit_rep)
    ## plot LM statistics along with confidence interval
    flm_rep <- fxlm(CNY ~ USD + JPY + EUR + GBP, data = cny_rep)
    scus_rep <- gefp(flm_rep, fit = NULL)
    plot(scus_rep, functional = supLM(0.1))
    ## add lines related to breaks to your plot
    lines(cit_rep)
    apply(cny_rep,1,function(x) sum(is.na(x)))
    view raw FrankelWei_2010_2016.r hosted with ❤ by GitHub

    We got breaks in 2010 and in 2015 (when China’s stock markets crashed). We would have hoped for more breaks (we can still get them), but that would depend on the parameters chosen for our regression.

     

  • Google’s New Deep Learning MOOC Using TensorFlow

    Google’s New Deep Learning MOOC Using TensorFlow

    Deep learning became a hot topic in machine learning in the last 3-4 years (see inset below) and recently, Google released TensorFlow (a Python based deep learning toolkit) as an open source project to bring deep learning to everyone.

    deep_learning_google_trends
    Interest in the Google search term Deep Learning over time

    If you have wanted to get your hands dirty with TensorFlow or needed more direction with that, here’s some good news – Google is offering an open MOOC on deep learning methods using TensorFlow here. This course has been developed with Vincent Vanhoucke, Principal Scientist at Google, and technical lead in the Google Brain team. However, this is an intermediate to advanced level course and assumes you have taken a first course in machine learning, or that you are at least familiar with supervised learning methods.

    Google’s overall goal in designing this course is to provide the machine learning enthusiast a rapid and direct path to solving real and interesting problems with deep learning techniques.

    What is Deep Learning?

    Course Overview

  • Data Manipulation in R with dplyr – Part 3

    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

     

     

  • Data Manipulation in R with dplyr – Part 2

    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.