additional.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "R (system-wide)",
      "language": "r",
      "metadata": {
        "cocalc": {
          "description": "R statistical programming language",
          "priority": 10,
          "url": "https://www.r-project.org/"
        }
      },
      "name": "ir",
      "resource_dir": "/ext/jupyter/kernels/ir"
    },
    "language_info": {
      "codemirror_mode": "r",
      "file_extension": ".r",
      "mimetype": "text/x-r-source",
      "name": "R",
      "pygments_lexer": "r",
      "version": "4.0.5"
    },
    "colab": {
      "provenance": [],
      "collapsed_sections": [],
      "private_outputs": true,
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/alicewu1/introR/blob/main/File2_Tutorial_1_Intro_R.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "g-1EX1EjUnB_"
      },
      "source": [
        "# Tutorial 1: Introduction to R\n",
        "### Summer School in Statistics for Astronomers XVII, June 2022  \n",
        "### Instructor: Eric Feigelson (Center for Astrostatistics, Penn State University)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "_FCJFBrqUnCG"
      },
      "source": [
        "R is a powerful software environment for data analysis, graphics, and especially statistical analysis. It is available free to the public at www.r-project.org with easily installed binaries for Linux, MacOS and Windows.  This notebook provides an introduction to R designed for students and researchers in astronomy.  No previous experience is needed, although some familiarity with scripting languages like Matlab or IDL is helpful.  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "f4ZwVqeEUnCH"
      },
      "source": [
        "### Running the Jupyter notebooks \n",
        "\n",
        "During the Summer School, we are running these Jupyter notebooks in Google's Colaboratory environment [here](https://colab.research.google.com/drive/1iz6ILnVGt8Qc6UR1l7oTPou4l6WSrw9S?usp=sharing).  To run them here, click the **Connect** button in the upper-right section of the page.  When you encounter a Jupyter cell with an R script, click the black circle with arrow button to run the code. You are encouraged to investigate options and learn further by editing the script; again click the black circle to run the cell. \n",
        "\n",
        "You can also run them in Python on your local computer; Python's *anaconda* distribution and *conda* package management system automatically include Jupyter Notebooks as an application accessed through the *Anaconda Navigator*.  However, it is possible that your Anaconda release did not automatically include the R kernel for Jupyter.  You can tell whether it is present if the R logo (blue R in front of gray ellipse) is shown at the top-right of the Jupyter page.  If it is missing, then the code cells below will fail upon execution. To install the R kernel within Python, type: **conda install -c r r-irkernel**.  Alternatively, you can make R available to Jupyter from an R console by typing: **install.packages('IRkernel')  ;  IRkernel::installspec()**.  For further discussion of the relationship between Python, R, Jupyter and other notebook environments, see [here](https://www.datacamp.com/community/blog/jupyter-notebook-r). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "W9c2hitGUnCH"
      },
      "source": [
        "Some basic operational information:\n",
        "- R is driven by the command line in the R console.  Double-click the R icon or (in Linux) type 'R' in the terminal to open the console. \n",
        "- R has >100,000 'functions' that perform various tasks. Each function uses () brackets to list arguments and parameters.  \n",
        "- Every function has a 'help(fn)' page telling how it is used, what operation is run, and what output it produces with references and examples that can be cut-and-paste into an R console. It is essential for both novice and experienced R programmers to read help files.\n",
        "- A hash mark ( # ) denotes comments in an R script.  A semi-colon ( ; ) has the same action as a carriage return\n",
        "- The commands in the Summer School tutorials can be run interactively in Jupyter, or cut-and-pasted into a separate R console."
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**In Section I** we examine and set up our software environment. *sessionInfo()* tells the operating system, and basic pre-installed packages. Try this both in Google Colaboratory and on your local R console.  *getwd()* and *setwd()* are commonly used functions at the beginning of a session -- *wd* stands for 'working directory'.  Google Colaboratory does not allow changing directories, but you will want to create a directory for our Summer School work on your local computer.  *library()* pops up a new window listing the CRAN packages currently available to you.  Remember there are ~19000 more!  Finally, we run *citation()* that gives the citation to R for publications. Every CRAN package has a corresponding *citation(package)* that should be used for publications."
      ],
      "metadata": {
        "id": "96opTdEahmMu"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rrFK4kPWUnCI",
        "collapsed": true
      },
      "source": [
        "# I.  Set up your session\n",
        "\n",
        "sessionInfo()  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "getwd()                     # find working directory.  \n",
        "#setwd('/Users/ericfeigelson/Desktop/Rdir')  # set working directory if you run \n",
        "                                             # locally on your computer\n",
        "#getwd()                     # see the working directory has changed"
      ],
      "metadata": {
        "id": "NJJBawVQjKaP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RV6FEm_xUnCJ"
      },
      "source": [
        "library()                   # see packages installed on your computer\n",
        "                            # ~30 are installed automatically with R\n",
        "citation()                  # quote this citation in any publication using R"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**Section II** creates a simple R object -- a vector of four integers -- and functions that tell its characteristics.  *ls()* is useful to see what objects are in your current session.  The *str* function is very useful for seeing the structure of complex R objects produced by advanced statistical functions. Note that the content of an object is displayed on the console by just stating its name; no *print* or *show* command is needed. The *cat* function gives a more elegant way to display on the console. \n",
        "\n",
        "*write* and *save* bring R objects out of your session onto your disk in ASCII and binary formats respectively. CRAN package *FITSio* provides output in astronomers' FITS format.   \n",
        "\n",
        "R vectors and arrays start at index 1, unlike Python where the first element has index 0.  This incompatibility can complicate programs that entangle the two languages.  It is very useful to [link text](https://)"
      ],
      "metadata": {
        "id": "7WdmwwGSjWE_"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cqdb23DtUnCJ"
      },
      "source": [
        "# II. Create and characterize a vector\n",
        "\n",
        "a <- c(33, 44, 92, 58)      # combine numbers into a vector\n",
        "length(a)\n",
        "sum(a)\n",
        "median(a)  ;  mad(a)\n",
        "ls()                        # list names of objects in your environment\n",
        "class(a)                    # state the `class' of an R object (described in III below)\n",
        "str(a)                      # state the structure of an R object\n",
        "summary(a)                  # many R objects have a built-in 'summary' function"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I3ZS9Gf-UnCK"
      },
      "source": [
        "a                           # display the contents of an R object on the console\n",
        "# Annotated write to console. The \\n symbol is a carriage return.\n",
        "cat('Sum of ', length(a), ' elements in the vector a = ', sum(a), '\\n')\n",
        "\n",
        "write(file='output', a)     # write an ASCII file into the working directory\n",
        "save(file='output_bin', a)  # write a binary file\n",
        "save(file='output_bin', a)  # error because 'Save' is not a known function.\n",
        "                            # R syntax is case sensitive."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1rmy1In3UnCL"
      },
      "source": [
        "# Manipulation of vector indices\n",
        "\n",
        "a[1:4]          # show all four elements of the vector\n",
        "a[3]            # note that R vectors start with index 1 not 0, unlike Python\n",
        "a > 40          # logical operation\n",
        "sum(a[a>40])    # note the self-referential use of vector/array indices here\n",
        "which.max(a)    # R has many built-in functions for vectors and arrays\n",
        "match(44, a)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "4WDGHRRNUnCL"
      },
      "source": [
        "****\n",
        "\n",
        "**Section III** highlights classes and packages.  Base-R places objects into several: numeric, character, logical, vector, matrix, factor, data.frame, list, and dozens of others designed by advanced R functions and CRAN packages. *plot, print, summary* functions are adapted to class objects; see e.g. *methods(summary)*.  \n",
        "\n",
        "Two particularly important classes are the 'data frame' used for tabular data and the 'list' used as a bucket with heterogeneous content.  The data frame is a 2-dimensional array with associated column names. The list class allows a hierarchical structure of R objects such as scalars, vectors, arrays, and attributes.  Here we make a hierarchical list, use 'str' (structure) to show its contents, and access an element of the list using the $ delimiter. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KiC-8vHEUnCM"
      },
      "source": [
        "# III.  R classes and packages\n",
        "\n",
        "# Make and write a data.frame, a 2D array with column names\n",
        "\n",
        "d <- data.frame(cbind(seq(1:4), a, a^3))  # Bind columns into data frame\n",
        "class(d)\n",
        "names(d) <- c('ID', 'a', 'a_cubed') # Column names for data frame                                         \n",
        "d2 <- d[-4,-1]                            # Remove 4th row and 1st column\n",
        "d ; d2\n",
        "write.table(d, file='d.txt', quote=FALSE, row.names=FALSE)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AEItDAn_UnCM"
      },
      "source": [
        "# Make and show a list.\n",
        "\n",
        "b_list <- list(star=c('Sirius', 'Procyon'), SpTy=c('O','B','A'), Hubble_km.s=68)\n",
        "str(b_list)\n",
        "b_list[['SpTy']] = list(subtype=seq(0.1:0.9, by=0.1))\n",
        "str(b_list)\n",
        "b_list$SpTy$subtype[1:3]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "tOey_6CGUnCM"
      },
      "source": [
        "The 19,000+ CRAN packages are obtained on-the-fly as you need them. The *install.packages()* function, needed only once in the lifetime of your computer, downloads the chosen package from *https://cloud.r-project.org* or a mirror site.  But the *library()* command is frequently needed to bring the package into the current R session.  \n",
        "\n",
        "There is no useful index of the 19,000+ CRAN packages and their many functions. The collection expanded exponentially during the 2000s and is continues to grow by several packages every day.  Expert volunteers in ~40 statistical areas update lists of CRAN packages in their area; these are accesses on the Web at [CRAN Task Views](https://cran.r-project.org/web/views/). Task Views of particular interest to astronomers include _Bayesian_, _Cluster_, _HighPerformanceComputing_, _MachineLearning_, _Multivariate_, _Optimization_, _Robust_, _Spatial_, _Survival_, and _TimeSeries_. \n",
        "\n",
        "Astronomy-specific packages (e.g. stellar evolutionary tracks, celestial mechanics) are listed in the _ChemPhys_ CRAN Task View. The package _FITSio_ reads astronomical FITS formatted files and headers, converting them into an R list. The package _astrolibR_ is a collection of long-established functionalities useful in astronomical research, translated from the _IDL Astronomy Library_.\n",
        "\n",
        "Unfortunately, the Google Colaboratory does not allow installation of new packages by the user.  So portions of our notebooks are best run by cut-and-paste into an R console on your local computer."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "Uz6Wyl8CUnCN"
      },
      "source": [
        "This is a good moment to examine some R help files.  Try *help(sum)*, *help(cbind)*, and *help(mad)* for functions we have already used.  All R users, both novice and expert, are constantly looking at help files to learn what functions do, the syntax of their use, and related functions. \n",
        "\n",
        "R help files give essential information on all functions in a standard format:  \n",
        "1. The top lines give the package where the function resides and a brief description.  \n",
        "2. _Usage_ gives the list of inputs for the function.  Input parameters with an _=_ have a default and do not need to be specified by the user.  \n",
        "3. _Arguments_ describes these input parameters.  \n",
        "4. _Details_ summarizes the functionality including formulae and algorithms. \n",
        "5. _Value_ gives the output of the function.  Typically the program specifies _outfn <- fn(x,y,z, option='special')_ so the full list of output values are available for use, such as _plot(outfn\\$x, outfn\\$y)_. \n",
        "6. _References_ to published literature where the function is described.\n",
        "7. _See also_ links to R functions with related purpose to the current function.\n",
        "8. _Examples_ show usage of the function, often with a built-in dataset.  Examples in R help files can always be cut-and-pasted into any R console.  "
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**Section IV** shows that R knows elementary math.\n"
      ],
      "metadata": {
        "id": "VU8gLgEnoPOS"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ut-dA1IMUnCO"
      },
      "source": [
        "# IV. Arithmetic, algebra, trigonometry, and formatting numerics\n",
        "\n",
        "5 + 10003\n",
        "5-3 ; 5*3 ; 5/3 ; 5^3\n",
        "x <- 5 ; y <- 3\n",
        "x+y "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wYpS8tY7UnCO"
      },
      "source": [
        "sin(0)  ; sin(pi/2)         # note angles are in radians\n",
        "ang <- seq(0, pi/2, length=30) \n",
        "sin(ang)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mGbJ2uXzUnCO"
      },
      "source": [
        "trunc(12345.6789) ; round(12345.6789)  \n",
        "format(12345.6789, digits=2, scientific=TRUE)   \n",
        "\n",
        "log(20)  ;  log10(20)   # log() in R is base-e. Use log10() for base-10 logarithms."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "1IuCkgR1UnCP"
      },
      "source": [
        "> **Exercise:** Practice R syntax.  Practice with some elementary R functions: arithmetic and algebra, 2-dimensional array manipulation, producing multi-element lists.  Write a brief program exercising program flow control (if, ifelse, when, repeat, break, stop, etc.). Use help files!"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Exercise: Practicing R Syntax\n",
        "# ELEMENTARY FUNCTIONS\n",
        "832 + 106\n",
        "23-9; 34*7; 100/5; 2^5\n",
        "a <- 15\n",
        "b <- 10\n",
        "a*b\n",
        "c <- c(23, 9, 45, 94, 106)\n",
        "\n",
        "# 2-D ARRAY MANIPULATION\n",
        "e <- data.frame(cbind(seq(1:5), a, a^3, b, c))  # Bind columns into data frame\n",
        "class(e)\n",
        "names(e) <- c('ID', 'a', 'a_cubed', 'b', 'c') # Column names for data frame                                         \n",
        "e2 <- e[-5,-1]                            # Remove 5th row and 1st column\n",
        "e ; e2\n",
        "write.table(d, file='d.txt', quote=FALSE, row.names=FALSE)\n",
        "\n",
        "\n",
        "# MULTI-ELEMENT LIST\n",
        "\n",
        "f_list <- list(music=c('Pop', 'Rock'), subgenre=c('alt rock','math rock','rock n roll'))\n",
        "str(b_list)\n",
        "f_list[['subgenre']] = list(subtype=seq(0.1:0.9, by=0.1))\n",
        "str(f_list)\n",
        "f_list$subgenre$subtype[1:4]"
      ],
      "metadata": {
        "id": "YDcot8JL5jIs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**Section V** shows that R can be used for calculations based on analytical equations from mathematics, physics or astrophysics.  Here we calculate galaxy luminosity distances for a Euclidean and a $\\Lambda$CDM relativistic universe.  This is our first exercise of the *function()* function, allowing the user to create a new procedure for later use.  The content of a new function lies within curly brackets, {}, and can be simple or complicated.   \n",
        "\n",
        "We also have our first experience with R's basic 2-dimensional plotting function, *plot()*.  It is wonderfully useful and flexible giving publication-quality graphics.   For its many option, see help(par) ... the longest help file I know!    "
      ],
      "metadata": {
        "id": "BcMZayxTotlr"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qhrHHou1UnCP"
      },
      "source": [
        "# V. Astrophysical calculations of galaxy distances \n",
        "\n",
        "# The `function' function: Note how one function uses another\n",
        "# This how R builds new capabilities based on old capabilities in a compact syntax.\n",
        "\n",
        "# First, make a simple calculation without functions: galaxy distances in Mpc for six redshifts  \n",
        "\n",
        "z <- seq(0.0, 0.5, 0.1)\n",
        "z\n",
        "H_0 <- 68  \t \t\t               \t# km/s/Mpc,  Planck value\n",
        "speed.light <- 3.0E5          \t# km/s\n",
        "dist <- speed.light*z / H_0     # in Mpc\n",
        "dist\n",
        "class(dist)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "t-fR30Z6UnCP"
      },
      "source": [
        "# Now, make a more complicated relativistic calculation using new functions.\n",
        "\n",
        "Omega_m <- (0.022068 + 0.12029) / (H_0/100)^2    \n",
        "Omega_Lambda <- 0.6825                    # Planck satellite values\n",
        "\n",
        "E.H0 <- function(redshift) {sqrt(Omega_m*(1+redshift)^3 + Omega_Lambda)}\n",
        "\n",
        "lum.dist <- function(redshift) {  \n",
        "\tluminosity.distance = (speed.light/H_0) * integrate(E.H0, 0, redshift)$value \n",
        "\treturn(luminosity.distance)  }\n",
        "\n",
        "distGR <- Vectorize(lum.dist)(z)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TQOCUjTGUnCP"
      },
      "source": [
        "# Plot the results\n",
        "# The 'plot' function has extensive options to change format; see 'help(par)'\n",
        "\n",
        "options(jupyter.plot_scale=1)\n",
        "options(repr.plot.width = 7, repr.plot.height = 5)\n",
        "\n",
        "plot(z, distGR, type='l', lty=2, lwd=2, ylab='Distance (Mpc)')\n",
        "lines(z, dist, lty=1, lwd=2)   \n",
        "legend(0.0, 2500, lty=c(1,2), lwd=c(2,2), title='Galaxy distances', \n",
        "     \tlegend=c('Euclidean', expression(Lambda*CDM)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "kT9aRCNHUnCQ"
      },
      "source": [
        "> **Exercise** Integrate an astrophysical function.  Estimate the age of the Universe as a function of redshift for a standard $\\Lambda$CDM universe model: \n",
        "$$ t(z) = H_0^{-1} \\int_{z}^{\\infty}\\frac{dz'}{(1+z')h(z')} $$ \n",
        "where $ h(z) = \\sqrt{(1-\\Omega_{total})(1+z)^2 + \\Omega_m(1+z)^3 + \\Omega_{\\Lambda} } $, $\\Omega_m$ is the matter density parameter, and $\\Omega_{\\Lambda}$ is the dark energy density parameter. Plot the age of the Universe vs. redshift ($z=0$ to 10) for three hypothetical universes: matter-dominated ($\\Omega_m=1.0$ and $\\Omega_{\\Lambda}=0.0$), dark-energy-dominated ($\\Omega_m=0.01$ and $\\Omega_{\\Lambda}=0.99$), and a realistic universe ($\\Omega_m=0.31$ and $\\Omega_{\\Lambda}=0.69$).  This problem and solution is courtesy of graduate student Phoebe Sandhaus, Penn State."
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**Section VI** gives our first astrostatistical calculation based on a (fictitious) dataset. Here we seek the distribution of galaxy X-ray luminosities from a sample of 12 galaxies where four are undetected. In statistical parlance, the astronomers' 'upper limits' are called 'left-censored' data points and are treated in a field of applied statistical called 'survival analysis'.  One of the most important results of 20th century statistics was the theorem by [Kaplan & Meier (1958)](https://www.jstor.org/stable/2281868) deriving the unique nonparametric maximum likelihood estimator for a randomly censored univariate dataset.  The KM paper has 62K citations and ~6M Google hits; the KM estimator has been used in ~800 astronomical studies.\n",
        "\n",
        "Although base-R has an implementation of the KM estimator, we will use the one provided by CRAN pacakage *NADA*.  The methodology and cookbook are detailed in the volume 'Statistics for Censored Environmental Data Using Minitab and R' (D.Helsel, 2nd ed, 2012) and is introduced in Chpt. 10 of our textbook 'Modern Statistical Methods for Astronomy with R Applications' (Feigelson & Babu 2012).  \n",
        "\n",
        "We start by installing *NADA*, one of many survival CRAN packages.  **This installation is not permitted by Google Colaboratory, so this section should be cut-and-pasted into your local R console.fco** We create a fictitious dataset of X-ray luminosities and apply NADA's *cenfit* function.  We use *str* to see the output from *cenfit* and various simple displays given in the *help(cenfit)* file.   "
      ],
      "metadata": {
        "id": "sJFAGDG6tOwq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# VI.  Calculate the Kaplan-Meier estimator of the galaxy X-ray Luminosity Function from a small dataset.\n",
        "\n",
        "install.packages('NADA')\n",
        "library(NADA)\n",
        "\n",
        "help(cenfit)\n",
        "galLx <- c(39.7, 42.3, 41.7, 41.6, 43.0, 40.2, 40.6, 41.1, 39.9, 40.7, 41.2, 42.9)\n",
        "cens <- c(F,F,F,F,F,F,F,F, T,T,T,T)  # FALSE=detected TRUE=left-censored. \n",
        "galKM <- cenfit(galLx, cens)\n",
        "str(galKM)   # The help(cenfit) should describe this output but does not!  \n",
        "\n",
        "summary(galKM)  # display KM estimator on console\n",
        "plot(galKM, xlab=expression(\"log L\"[x]~~\"(erg/s)\"))     # plot KM estimator with 95% confidence band\n",
        "quantile(galKM) # median and other quantiles\n",
        "mean(galKM)  ;  sd(galKM)  # mean, uncertainty of mean, standard deviation about the mean"
      ],
      "metadata": {
        "id": "PHB2RSJ_wY-C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**In Section VII** we create a fictitious bivariate *(x,y)* dataset of 500 points where a nonlinear relation between x and y is present.  Also the errors are heteroscedastic; that is, the scatter in *y* depends on the *x* variable.  We will use this dataset here and in the next tutorial.  \n",
        "\n",
        "Here we show several common plots --  scatterplot, box plot, histogram, cumulative distribution function. Multipanel plots are designed using the *par(mfrow)* plot parameter; see *help(par)*. \n",
        "\n",
        "Plots can be downloaded to disk using functions like *pdf()*, *png()*, and *jpeg()* where the user presets the plot window.  But a convenient function is *dev.copy2pdf()* where the current window is downloaded.  "
      ],
      "metadata": {
        "id": "iJ1QXYBS0joC"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-u3IR1XwUnCQ",
        "collapsed": true
      },
      "source": [
        "# VII. Examine, summarize and plot bivariate and univariate data \n",
        "\n",
        "\n",
        "set.seed(1)\n",
        "x <- sample(seq(0.01, 3, length.out=500))\n",
        "y <- 0.5*x + 0.3^(x^2) + rnorm(500, mean=0, sd=(0.05*(1+x^2)))\n",
        "xy <- cbind(x, y)\n",
        "\n",
        "plot(xy, pch=20)\n",
        "summary(x) ; summary(y)   \t# Summarizes properties of an R object"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M67sbdDkUnCQ"
      },
      "source": [
        "par(mfrow=c(1,2))  \t\t# Set up a two-panel figure\n",
        "boxplot(x,  notch=T, main='Boxplot for X')  # See help(boxplot.stats) for explanation\n",
        "boxplot(y,  notch=T, pch=20, cex=0.5, main='Boxplot for Y')\n",
        "dev.copy2pdf(file='box.pdf')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ui3EDPA3UnCQ"
      },
      "source": [
        "par(mfrow=c(1,1))\n",
        "hist(x, breaks=30, main='', xlim=range(x), ylim=c(0,100), \n",
        "     xlab='Yvar', col='royalblue4')\n",
        "hist(y, breaks=30, main='', xlab='Yvar', \n",
        "     col='#ee3b3b70', add=TRUE) # add=TRUE suppresses a new plot"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "loN9vpOCUnCR"
      },
      "source": [
        "plot(ecdf(x), pch=20, cex=0.0, verticals=TRUE, main='',ylab='EDF',xlab='')  \n",
        "plot(ecdf(y), pch=20, cex=0.0, verticals=TRUE, add=T)\n",
        "text(2.0,0.5,\"X\") ; text(1.4,0.8,\"Y\")             # text adds annotation within a plot\n",
        "dev.copy2pdf(file='ecdf.pdf')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "****\n",
        "\n",
        "**Section VIII** gives more practice with vectors and arrays. We show how random samples can be drawn from a dataset. Powerful theorems show that *bootstrap resampling* -- random draws with replacement -- can often give reliable estimates of confidence intervals for statistical quantities derived from the dataset."
      ],
      "metadata": {
        "id": "A8yklMY03dFU"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8wNgDzsLUnCR"
      },
      "source": [
        "# VIII. Arrays, data frames and filtering\n",
        "\n",
        "# Here xy is an `array' of numbers created by `column bind'\n",
        "\n",
        "xy <- cbind(x, y)  ;  str(xy) \t\n",
        "\n",
        "# A data.frame associates names to the columns\n",
        "\n",
        "xy <- as.data.frame(xy)         \n",
        "names(xy) <- c('Xvar', 'Yvar')  \n",
        "\n",
        "# Collect rows where the first column value exceeds 2\n",
        "\n",
        "high_x1 <- xy[xy[,1]>2,]        \n",
        "high_x2 <- subset(xy, xy[,1]>2)\t# Another way to extract rows\n",
        "setequal(high_x1, high_x2)      # test equality of two vectors"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qHuWXQX8UnCR"
      },
      "source": [
        "# Sampling and bootstrapping \n",
        "\n",
        "trials <- sample.int(length(xy[,1]),20) # 20 random rows\n",
        "xy[trials,]\t\n",
        "\n",
        "# 20 bootstrap resamples\n",
        "\n",
        "trials <- sample.int(length(xy[,1]),20, replace=T) \n",
        "xy[trials,]\t\n",
        "\n",
        "# Estimate the standard error of the median of Yvar\n",
        "\n",
        "median(xy[,2]) \n",
        "\n",
        "# Median absolute deviation estimate of median s.e.\n",
        "\n",
        "mad(xy[,2]) / sqrt(500)\t"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "H2wHjvLBUnCS"
      },
      "source": [
        "library(boot)  \t\t\t# The following function in a base-R library\n",
        "med <- function(x,index){ median(x[index]) }\n",
        "\n",
        "# Read help(boot) to understand its output list structure\n",
        "\n",
        "boot(xy[,2], med, R=1000) # Bootstrap estimate of median s.e. \n",
        "hist(boot(xy[,2], med, R=1000)$t, breaks=50, \n",
        "     xlab='Bootstrap median of Yvar')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "aUfbSQAdUnCS"
      },
      "source": [
        "> **Exercise**  Use bootstrap resampling (random sampling with replacement) to estimate uncertainties of a statistic. Create a univariate sample with a weird distribution … maybe sampling from a polynomial or nonlinear function over some interval.  First, calculate the median and a robust measure of its standard error: $1.48*(MAD)$ where MAD is the median absolute deviation and the 1.48 scales it to the standard deviation for a Gaussian distribution.   Second, estimate the uncertainty of the median from a bootstrap resampling.  Give the standard error, 95% confidence interval, and plot a histogram of the bootstrap medians."
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Section IX** gives more practice with bivariate plotting, and shows application of two correlation coefficients that give a probability of relationship between the *x* and *y* variables.  Pearson's correlation coefficient is commonly used but is limited: it assumes a linear relationship and requires homoscedastic Gaussian errors. The nonparametric Kendall's tau coefficient makes no assumptions and is therefore more generally useful.  "
      ],
      "metadata": {
        "id": "lVa4nFvV4YAs"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FcTA6VYIUnCS"
      },
      "source": [
        "# IX. Bivariate plots and correlation tests\n",
        "\n",
        "# Scatterplot. See help(points) for symbol shapes defined by the 'pch' argument.\n",
        "par(mfrow=c(1,2))\n",
        "plot(xy, pch=20, cex=0.5)\t \n",
        "plot(log10(xy), pch=20, cex=0.5, xlab='log(Xvar)', ylab='log(Yvar)')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fYHAd0yDUnCS"
      },
      "source": [
        "length(x[x>2])\t\t# State length of a vector.  Use `dim' for an array or data.frame.\n",
        "# Parametric hypothesis test for bivariate correlation \n",
        "cor.test(x[x>2],y[x>2], method='pearson')\t\n",
        "# Nonparametric hypothesis test for bivariate correlation\n",
        "cor.test(x[x>2],y[x>2], method='kendall')\t"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "iMD_GRF6UnCT"
      },
      "source": [
        "****\n",
        "\n",
        "**Section X: R for advanced astronomical data analysis**\n",
        "\n",
        "R was designed in the 1990s for interactive data exploration and analysis by a single user using a single CPU core with a graphics window.  The compiler was not optimized for speed.  However, around 2012 a new compiler was implemented similar to the one underlying Python.  R and Python thus have roughly similar computational efficiency.  Both are slower than C, C++ and Fortran which are used in both environments for computationally heavy portions of an analysis. Many CRAN packages use C++ under the hood.  \n",
        "\n",
        "In addition to improved compilers, R has also been adapted for use in a multi-core environment, including GPUs.  Two capabilities are now in base-R (Rmpi and snow), but most are in several dozen CRAN packages; see the [CRAN Task View on High Performance and Parallel Computing with R](https://cran.r-project.org/web/views/HighPerformanceComputing.html).  \n",
        "\n",
        "Generally, R has not been viewed as a language appropriate for pipeline analysis of large datasets such as those emerging from astronomical satellites or observatories.  This criticism is valid in the sense that it does not have a large and growing software library of specialized code for astronomical data, such as the ~2000 CRAN packages in the _Bioconductor_ toolbox for genomics or the diverse packages in the _astropy_ project in Python.  \n",
        "\n",
        "However, R can be operated in batch mode, running scripts with the _Rscript_ (or with more flexibility, CRAN package _littler_) command line interface from the operating system prompt.  These commands can be embedded in parallelized scripts for supercomputers using PBS or related scripts.  Using Rscript, we at Penn State have reanalized the full dataset of major surveys from NASA satellites and ground-based observatories involving millions of input and output files.  The motivation for this approach is that the instrument-specific characteristics of the data reduction had already been performed in Python (e.g NASA Levels 2 and 3 data products), and the data were ready for advanced statistical analysis using CRAN packages.  In our case, sophisticated CRAN packages originally designed for econometrics and machine learning play a central role, along with a Fortran code we wrote for a computationally intensive step.  R code for astronomical pipelines can be remarkably brief if one effectively uses the tens-of-thousands of advanced functions available in CRAN.\n",
        "\n",
        "Many statistical functionalities of R can also be run directly from Python and several other languages through bi-language interfaces.  R scripts can be embedded in Python using [rpy2](https://pypi.org/project/rpy2/) and, conversely, Python programs can be embedded in R scripts using CRAN package [reticulate](https://cran.r-project.org/web/packages/reticulate/index.html) or several other options.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "3O-QjliYUnCT"
      },
      "source": [
        "# XI.  Resources for further study of R\n",
        "\n",
        "We give here a miscellaneous collection of commands and resources useful for learning more about R as a programming language.  \n",
        "\n",
        "* There are >700 books with 'R' in the title, most presenting both methodology and code tutorials.  A new book on R is published every ~10 days.  Two are devoted specifically to astronomy:\n",
        "  * [_Modern Statistical Methods for Astronomy with R Applications_](https://doi.org/10.1017/CBO9781139015653), 2012 E. D. Feigelson & G. J. Babu, Cambridge Univ Press\n",
        "  * [_Bayesian Models for Astrophysical Data Using R, JAGS, Python, and Stan_](https://doi.org/10.1017/CBO9781316459515), 2018, J. M. Hilbe, R. S. de Souza & E. E. O. Ishida, Cambridge Univ Press\n",
        "\n",
        "* Two high-quality introductions to R:\n",
        "  * From the R Core Team: https://cran.r-project.org/doc/manuals/R-intro.html\n",
        "  * From Carnegie-Mellon University: http://www.stat.cmu.edu/~cshalizi/statcomp/14/\n",
        "\n",
        "* There are vast informal online learning resources about R programming: \n",
        "  *  [R-bloggers](https://www.r-bloggers.com) aggregates entries from ~1100 blogs.  \n",
        "  * [Stack Overflow] has ~400,000 questions and answers about R programming.\n",
        "  * CRAN packages are often described in two open-access journals: [Journal of Statistical Software](https://www.jstatsoft.org/index) and [The R Journal](https://journal.r-project.org).  These articles often appear as _vignettes_ within the R environment.  \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SIAhSLscUnCU"
      },
      "source": [
        "# A list of the ~30 important CRAN packages embedded in the base-R environment\n",
        "library()\n",
        "\n",
        "# A full list of ~400 functions in R's `base' package\n",
        "library(help = \"base\")\n",
        "\n",
        "# Statistics in base R (~400 functions, tens-of-thousands more in CRAN and elsewhere in R)\n",
        "library(help='stats')\n",
        "\n",
        "# List current contents of your session environment\n",
        "ls()\n",
        "\n",
        "# Programming utilities including:\n",
        "#    Use `source' to bring in external R scripts\n",
        "#    Use `edit' to edit an R object\n",
        "#    Use 'environment' to segregate a collection of objects\n",
        "#    Functions `debug', `trace' and `browser' assist with code testing\n",
        "#    Function 'process.events' allows low-level handling of R commands\n",
        "library(help = 'utils')\n",
        "\n",
        "# Loops:  for( i in 1:100) { ... }\n",
        "# Program flow control:  if/else, ifelse, switch, while, repeat, next, break, stop\n",
        "foo <- 2\n",
        "if(foo == 1) cat('Hello world!') else cat('Do nothing')\n",
        "for(i in 1:10) { cat(' Num = ', i, '\\n') }\n",
        "\n",
        "# Graphics and devices in base R (other packages in CRAN)\n",
        "library(help='graphics')\n",
        "library(help='grDevices')\n",
        "\n",
        "# Parallel computing control in base R \n",
        "# CRAN has dozens of other high performance computing packages\n",
        "library(help='parallel')\n",
        "\n",
        "# Run an R script residing on disk\n",
        "help(source)\n",
        "\n",
        "# Save R objects (or your full environment) onto disk\n",
        "help(save) ; help(load)\n",
        "\n",
        "# Save or load history of R commands\n",
        "help(savehistory)  ;  help(loadhistory)\n",
        "\n",
        "# Connections, pipes, sockets, URLs, clipboard, compression, etc.\n",
        "help(connections)\n",
        " \n",
        "# Interact with host computer \n",
        "Sys.info()\n",
        "system('ls -l')\n",
        "system.time(fft(seq(0,1,length.out=1000000)))\t# A million fast Fourier transforms\n",
        "\n",
        "# Construct composite strings using 'paste'\n",
        "# Extract postions of a string using `substring'\n",
        "band_ir <- 'J'\n",
        "paste('NGC1068',band_ir,'FITS', sep='.')\n",
        "\n",
        "# FITS format reader/writer\n",
        "install.packages('FITSio') ; library(FITSio)\n",
        "\n",
        "# IDL Astro Library translated into R\n",
        "install.packages('astrolibR') ; library(astrolibR)\n",
        "\n",
        "# R/CRAN functions are public domain and can be wrapped from Python\n",
        "# programs using package rpy2.  Example:\n",
        "### pip install rpy2\n",
        "### import rpy2\n",
        "### import rpy2.robjects as robjects\n",
        "### R = robjects.r\n",
        "### ranGauss = R.rnorm(100)\n",
        "### print ranGauss\n",
        "\n",
        "# Python code can be wrapped into R using CRAN package 'reticulate' (among others)\n",
        "# R has similar interfaces to many other languages."
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}