%\documentclass[handout]{beamer} \documentclass{beamer} \usetheme{Rochester} \usecolortheme{wolverine} %\definecolor{UBCblue}{rgb}{0.04706, 0.13725, 0.26667} % UBC Blue (primary) %\usecolortheme[named=UBCblue]{structure} \definecolor{LEEgreen}{rgb}{0.04,0.667,0.137} \usepackage{amsmath} \usepackage{multirow} \usepackage{latexsym} \usepackage{amssymb} \usepackage{pifont} %\usepackage{epic} \usepackage{graphicx} \usepackage{skull} %\usepackage{multicol} %\usepackage{cancel} %\usepackage{pdfpages} \usepackage{listings} \usepackage{mathtools} \usepackage{fancybox} \usepackage{ulem} \DeclarePairedDelimiter\abs{\lvert}{\rvert} \addtolength{\parskip}{1ex} \setbeamercovered{transparent=5} \definecolor{ForestGreen}{cmyk}{0.5, 0, 0.5, 0.5} \definecolor{DarkBlue}{rgb}{0.0, 0.08, 0.45} \definecolor{BrickRed}{cmyk}{0, .89, .94, .28} \definecolor{UBCblue}{rgb}{0.04706, 0.13725, 0.66667} % UBC Blue (primary) \def \hi{\textcolor{BrickRed}} \def \E{{\textrm{E}}} \newcommand{\ex}[2][]{\ensuremath{\mathrm{E}_{#1}\!\left[#2\right]}} \newcommand{\var}[2][]{\ensuremath{\mathrm{Var}_{#1}\left(#2\right)}} \newcommand{\cov}[2][]{\ensuremath{\mathrm{Cov}\left(#2\right)}} \newcommand{\cor}[2][]{\ensuremath{\mathrm{Cor}\left(#2\right)}} \newcommand{\given}{\ensuremath{\;|\;}} \def \R{\mathbb{R}} \def \proof{\noindent{\bf Proof: }} \def \eproof{\qed\vskip .3cm} \def \ss{S} \def \union{\cup} \def \intersect{\cap} \newcommand{\vsp}{\vspace{0.3cm}} \newcommand{\vspl}{\vspace{0.15cm}} \newenvironment{answers}{\color{BrickRed} ~}{\color{black} ~} \newcommand{\pr}[2][]{\ensuremath{\mathrm{Pr}_{#1}\!\left(#2\right)}} \newcommand{\perm}[2]{\operatorname{P}^{#1}_{#2}} \newcommand{\sol}{\setbeamercolor{structure}{fg=green!60!black}} \renewcommand{\vec}[1]{\underline{#1}} \title{MAS2906: Computational Probability and Statistics with R} \author[Lee Fawcett]{Newcastle University \\ \texttt{lee.fawcett@ncl.ac.uk}} \date{Semester 1, 2024/25} %\usepackage[usenames, dvipsnames]{xcolor} \begin{document} %%From web for code/highlighting \definecolor{mygreen}{rgb}{0,0.6,0} \definecolor{mygray}{rgb}{0.5,0.5,0.5} \definecolor{mymauve}{rgb}{0.58,0,0.82} \definecolor{light-gray}{gray}{0.95} \lstset{ backgroundcolor=\color{light-gray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}; should come as last argument basicstyle=\footnotesize, % the size of the fonts that are used for the code breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace breaklines=true, % sets automatic line breaking captionpos=b, % sets the caption-position to bottom commentstyle=\color{mygreen}, % comment style deletekeywords={...}, % if you want to delete keywords from the given language escapeinside={\%*}{*)}, % if you want to add LaTeX within your code extendedchars=false, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8 frame=single, % adds a frame around the code keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible) keywordstyle=\color{blue}, % keyword style language=R, % the language of the code morekeywords={*,...}, % if you want to add more keywords to the set numbers=left, % where to put the line-numbers; possible values are (none, left, right) numbersep=5pt, % how far the line-numbers are from the code numberstyle=\tiny\color{mygray}, % the style that is used for the line-numbers rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here)) showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces' showstringspaces=false, % underline spaces within strings only showtabs=false, % show tabs within strings adding particular underscores stepnumber=1, % the step between two line-numbers. If it's 1, each line will be numbered stringstyle=\color{mymauve}, % string literal style tabsize=2, % sets default tabsize to 2 spaces title=\lstname % show the filename of files included with \lstinputlisting; also try caption instead of title } \begin{frame} \frametitle{\small \textit{\color{purple}{Lecture: 2pm, Monday 23rd September}}} \titlepage \end{frame} \section{General Housekeeping} \begin{frame} \LARGE \textbf{Chapter 1} \textbf{General Housekeeping} \end{frame} \subsection{Lecturer information} \frame{ \frametitle{1.1 Lecturer information} The module leader for MAS2906 is Dr. Lee Fawcett \begin{itemize} \item Lee is based in room 2.07 on the second floor of the Herschel Building \pause \item If you have any questions about the module he can be contacted via email: \texttt{lee.fawcett@newcastle.ac.uk} \pause \item Feel free to post a query to the module Canvas page \pause \item Or you are welcome to schedule a face-to-face/Zoom meeting \end{itemize} } \frame{ \frametitle{1.1 Lecturer information} There will be several other people assisting with this module, primarily in the computer practical sessions. \vspace{0.3cm} They're all really friendly PhD students -- Give them a shout in class if you're stuck with anything! } \subsection{Classes} \frame{ \frametitle{1.2 Classes} In the first half of the module, until enrichment week, there will be a weekly lecture (1 hour) and a weekly computer lab (2 hours). It's vital you attend both! After enrichment week we start to consider more theoretical aspects of the module and so there will be an extra lecture scheduled every week to help with this. \begin{itemize} \item Before enrichment week: \begin{itemize} \item [--]Lecture: Monday 2-3, Herschel LT1 \item [--]Computer lab: Tuesday 1-3, Henry Daysh Computing Cluster 1.14 \underline{or} Wednesday 11-1, Herschel PC \end{itemize} \item After enrichment week: \begin{itemize} \item [--]Lectures: Monday \textbf{2-4}, Herschel LT1 \item [--]Computer lab: Tuesday 1-3, Henry Daysh Computing Cluster 1.14 \underline{or} Wednesday 11-1, Herschel PC \end{itemize} \end{itemize} } \subsection{Assessment} \frame{ \frametitle{1.3 Assessment} No written exam! \begin{itemize} \item Problem Solving Ex. 1 -- Numbas assignment (\textbf{15\%}) \pause \begin{itemize} \item [--] Covers basic R work from the first few weeks of the module \item [--] Completed within Numbas \item [--] Practice material will be provided \item [--] Will be open for a two week period, with a deadline of \textcolor{red}{\bf 4pm, Friday 1 November 2024} \end{itemize} \end{itemize} } \frame{ \frametitle{1.3 Assessment} \begin{itemize} \item Problem Solving Ex. 2 -- Mini Project (\textbf{15\%}) \pause \begin{itemize} \item [--] Covers some of the more theoretical aspects of the course \item [--] Will also require you to use some of your R skills from the first few weeks of the module \item [--] Project material will be released three weeks before the deadline of \textcolor{red}{\bf 4pm, Friday 6 December 2024} \item [--] Submission via Canvas upload \end{itemize} \end{itemize} } \frame{ \frametitle{1.3 Assessment} \begin{itemize} \item Numbas exam (\textbf{70\%}) \pause \begin{itemize} \item [--] Will cover a range of material from across the module \item [--] Will be completed within the Numbas application, but will also require a Canvas upload element \item [--] Sample test questions will be provided \item [--] Will take place during the \textcolor{red}{\bf January assessment period} -- precise details will be released later in the semester \end{itemize} \end{itemize} } \subsection{Late work policy} \frame{ \frametitle{1.4 Late work policy} \begin{center} See Canvas for details! \end{center} } \subsection{Course materials} \frame{ \frametitle{1.5 Course materials} All course materials will be posted to the dedicated Canvas page for MAS2906: \begin{itemize} \item A full copy of these \textcolor{blue}{\bf lecture notes}, \textit{with gaps} \pause \item All \textcolor{blue}{\bf slides} used to support lectures/computer labs after these have taken place \pause \item All \textcolor{blue}{\bf video recordings} \pause \item \textcolor{blue}{\bf Solutions} to assessed exercises \pause \item Any other materials, including additional handouts/worksheets, links to useful Numbas material, specimen exam papers and solutions, $\ldots$ \end{itemize} } \subsection{Books and other resources} \frame{ \frametitle{1.6 Books and other resources} Some useful books include: \pause \begin{itemize} \item W.N. Venables and D.M. Smith: An Introduction To R. Free to download at: {\small \begin{center} {\color{blue} \url{https://cran.r-project.org/doc/manuals/R-intro.pdf}} \end{center}} \item D. Braun and D.J. Murdoch: A First Course In Statistical Programming With R \item S. Ross: A First Course in Probability \item N.A. Weiss, P.T. Holmes and H. Michael: A Course In Probability \end{itemize}} \subsection{General background} \frame{ \frametitle{1.7 General background} In this module we will learn about computational tools for performing probability calculations and statistical analysis. \vsp We will use the R programming language. \vsp R is very widely used in universities, businesses and other organizations due to its many advantages over other similar software packages. \begin{itemize} \item It is available for \textcolor{blue}{\bf free} \item It is `\textcolor{blue}{\bf open source}' and can be used for a very wide range of statistical and scientific applications \item \textcolor{blue}{\bf Libraries} of R functions and data-structures have been created and are readily available on the web \end{itemize} } \frame{ \frametitle{1.7 General background} Learning to use R is important because: \begin{itemize} \item [(i)] It is the foremost tool in \textcolor{blue}{\bf modern computational statistics} \item [(ii)] It is a great way to learn general concepts in \textcolor{blue}{\bf programming} which can be used with other languages or in other contexts \item [(iii)] It can be used to illustrate ideas in \textcolor{blue}{\bf probability and statistics} \end{itemize} \pause In MAS1608 last year we focused mainly on point (ii) above; in this module, we focus mainly on points (i) and (iii). } \frame{ \frametitle{1.7 General background} \textbf{Rough structure of the course} \pause \begin{itemize} \item Weeks 1--3: R Refresher/Intro \item Weeks 4--6: Statistical Methods: Simulation \item Weeks 7--10: Statistical Methods: Transformations/Bivariate distributions \item Week 11: Revision \end{itemize} } \section{R Refresher} \begin{frame} \LARGE \textbf{Chapter 2} \textbf{R Refresher} \end{frame} \frame{ \frametitle{2.1 Preliminaries} R is installed on all University machines. \pause \vsp In this module we will be using \textcolor{blue}{\bf RStudio}, which is an R \textit{Integrated Development Environment}, or IDE. \vsp This just makes R much easier to use and navigate, through the four windows that you see when you open RStudio (text editor, console window, file explorer and graphics display). \pause \vsp You might remember from MAS1608 that both R and RStudio are free, so you can install them on your own computer. See {\small \begin{center} \url{http://www.ncl.ac.uk/maths-physics/students/resources/} \end{center}} or the Stage 2 induction pages on Canvas, for more details.} \frame{ \frametitle{2.1 Preliminaries} The first three weeks of this module are designed to recap the basics from MAS1608. \vsp Specifically, we will review: \begin{itemize} \item [(i)] Data types and vectors in R \item [(ii)] Dataframe manipulation \item [(iii)] Data summaries \item [(iv)] Functions and control statements \end{itemize} } \frame{ \frametitle{2.1 Preliminaries} You might recall the \texttt{movies} dataset from MAS1608. \vsp This was obtained from the internet movie database (IMDB): \url{http://imdb.com/}, which is a website devoted to collecting movie data supplied by studios and fans. \vsp More information can be found online: \begin{center} \url{http://imdb.com/help/show_leaf?about} \end{center}\pause This is a great dataset to use to demonstrate R's dataframe manipulation and querying capabilities -- it is large, has many variables, and is interesting. Everyone loves movies! } \frame{ \frametitle{2.1 Preliminaries} \input movietableforslides } \begin{frame}[fragile] \frametitle{2.1 Preliminaries} We have written an add-on R package for the module which will allow easy access to many interesting datasets. \vsp Installing the package is straightforward: % Turn off + signs in R listings % Other sweave tricks %\DefineVerbatimEnvironment{Sinput}{Verbatim} {xleftmargin=2em} %\DefineVerbatimEnvironment{Soutput}{Verbatim}{xleftmargin=2em} %\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em} %\fvset{listparameters={\setlength{\topsep}{0pt}}} %\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}} \begin{lstlisting} install.packages("ncldata", repos = "http://R-Forge.R-project.org", type = "source") \end{lstlisting} \pause To load the package, type: \begin{lstlisting} library(ncldata) \end{lstlisting} Then, to specifically load the movies dataset, type: \begin{lstlisting} data(movies) \end{lstlisting} \end{frame} \subsection{Data types in R} \begin{frame}[fragile] \frametitle{2.2 Data types in R} R has a variety of data types. For example: \begin{lstlisting} v = TRUE #%* \textcolor{ForestGreen}{`Logical' or `boolean' - TRUE or FALSE} *) w = "fred" #%* \textcolor{ForestGreen}{`Character' or `string'} *) x = 5.0 #%* \textcolor{ForestGreen}{`Double' or `numerical'} *) \end{lstlisting} \vspace{-0.5cm} There are also some `special' data types: \begin{lstlisting} y = 5/0 #%* \textcolor{ForestGreen}{Returns \texttt{Inf} }*) z = y-y #%* \textcolor{ForestGreen}{Returns \texttt{NaN} }*) \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.2 Data types in R} Another important data type is \texttt{NA}, short for `Not Available' and used for missing values. \vsp A summary of the data types in R is given in the Table below. \vsp {\small %\begin{table}[!h] \begin{tabular}{l|l l l l} \hline Type &Example 1 &Example 2 &Example 3 &Example 4 \\ \hline Doubles &2 &3.1242 &-45.6 &4e-10 \\ Logicals/Boolean&TRUE &FALSE & & \\ Characters &``Fred'' &``TRUE'' &``Male'' &``2'' \\ Infinity &Inf &5/0 & & \\ Other &NaN &NA & & \\ \hline \end{tabular} %\caption{A summary of some data types in R.} %\label{table:datatypes} %\end{table} } \end{frame} \begin{frame}[fragile] \frametitle{\textcolor{purple}{Quiz \#1}} In the code below, what is the value of x? \begin{lstlisting} x = 5 x = 2*x + 1 x + 5 \end{lstlisting} \includegraphics[width=1cm]{ironman.jpg} \begin{itemize} \item [1.] 5 \item [2.] 11 \item [3.] 16 \item [4.] 22 \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{\textcolor{purple}{Quiz \#1}} In the code below, what is the value of x? \begin{lstlisting} x = 5 x = 2*x + 1 x + 5 \end{lstlisting} \includegraphics[width=1cm]{ironman.jpg} \begin{itemize} \item [1.] 5 \item [2.] \textcolor{red}{\bf 11} \item [3.] 16 \item [4.] 22 \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{\textcolor{purple}{Quiz \#2}} In the code below, what type of variable is x? \begin{lstlisting} x = 5 x = 2*x + 1 x + 5 \end{lstlisting} \includegraphics[width=1cm]{ironman.jpg} \begin{itemize} \item [1.] Character \item [2.] Double \item [3.] Logical \item [4.] NA \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{\textcolor{purple}{Quiz \#2}} In the code below, what type of variable is x? \begin{lstlisting} x = 5 x = 2*x + 1 x + 5 \end{lstlisting} \includegraphics[width=1cm]{ironman.jpg} \begin{itemize} \item [1.] Character \item [2.] \textcolor{red}{\bf Double} \item [3.] Logical \item [4.] NA \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{\textcolor{purple}{Quiz \#3}} In the code below, what type of variable is y? \begin{lstlisting} x = 0 y = 2/0 y = "Fred" \end{lstlisting} \includegraphics[width=1cm]{ironman.jpg} \begin{itemize} \item [1.] NaN \item [2.] Double \item [3.] NA \item [4.] Character \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{\textcolor{purple}{Quiz \#3}} In the code below, what type of variable is y? \begin{lstlisting} x = 0 y = 2/0 y = "Fred" \end{lstlisting} \includegraphics[width=1cm]{ironman.jpg} \begin{itemize} \item [1.] NaN \item [2.] Double \item [3.] NA \item [4.] \textcolor{red}{\bf Character} \end{itemize} \end{frame} \subsection{Vectors} \begin{frame}[fragile] \frametitle{2.3 Vectors} Vectors are the most basic of all data structures, but are used in almost all R code. \vsp An R vector contains $n$ values of the same type, where $n$ can be zero. \vsp For example: \begin{lstlisting} c(0, 1, 2, 3, 4, 5) [1] 0 1 2 3 4 5 (first_vec = c(0, 1, 2, 3, 4, 5)) [1] 0 1 2 3 4 5 (second_vec = c("Male", "Female", "Male")) [1] "Male" "Female" "Male" x = 5.2 y = "Fred" \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.3 Vectors} There are special functions in R that can be used to determine the variable type. For example: \begin{lstlisting} is.double(x) [1] TRUE is.character(x) [1] FALSE is.character(y) [1] TRUE is.vector(first_vec) [1] TRUE \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.3 Vectors} And of course, we can use the \texttt{length} command to determine the length of a vector: \begin{lstlisting} length(first_vec) [1] 6 length(second_vec) [1] 3 \end{lstlisting} \pause And you can use the \texttt{seq} command to generate sequence vectors: \begin{lstlisting} (x1 = seq(1, 6)) [1] 1 2 3 4 5 6 (x2 = seq(-4, 4, by=2)) [1] -4 -2 0 2 4 (x3 = seq(3, -2, by=-0.5)) [1] 3.0 2.5 2.0 1.5 1.0 0.5 0.0 -0.5 -1.0 -1.5 [11] -2.0 \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.3 Vectors} Other useful commands for vectors: \begin{lstlisting} rev(x1) #%* \textcolor{ForestGreen}{Reverses the order of the vector}*) [1] 6 5 4 3 2 1 sort(x3) #%* \textcolor{ForestGreen}{Sorts the vector into ascending order}*) [1] -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 sort(second_vec) #%* \textcolor{ForestGreen}{Can be applied to characters: alphabetical order}*) [1] "Female" "Male" "Male" sum(x2); sum(x3); sum(second_vec) [1] 0 [1] 5.5 Error in sum(second_vec) : invalid 'type' (character) of argument unique(second_vec) #%* \textcolor{ForestGreen}{Lists all unique elements of the vector}*) [1] "Male" "Female" which(x2==2) #%* \textcolor{ForestGreen}{Gives index/indices of particular elements}*) [1] 4 \end{lstlisting} \end{frame} \subsubsection{Vector operations} \begin{frame}[fragile] \frametitle{2.3.1 Vector operations} When our data are in a vector structure, we can apply standard operations to the entire vector. For example: \begin{lstlisting} (x = seq(-4, 4)) [1] -4 -3 -2 -1 0 1 2 3 4 x*x [1] 16 9 4 2 1 0 1 2 4 9 16 x-5 [1] -9 -8 -7 -6 -5 -4 -3 -2 -1 x + x [1] -8 -6 -4 -2 0 2 4 6 8 \end{lstlisting} \end{frame} \subsubsection{Extracting elements from vectors}\label{section:vectorextraction} \begin{frame}[fragile] \frametitle{2.3.2 Extracting elements from vectors} R has a number of useful methods that we can use to extract subsets of our data. \vsp For example, square brackets: \begin{lstlisting} first_vec[2] [1] 1 second_vec[2:3] [1] "Female" "Male" first_vec[4:2] [1] 3 2 1 (z = 1:10) [1] 1 2 3 4 5 6 7 8 9 10 z[z<3 | z>=7] #%* \textcolor{ForestGreen}{Values of z that are less than 3 \textbf{or} at least 7}*) [1] 1 2 7 8 9 10 \end{lstlisting} \vspace{-0.5cm} \pause We will look at extraction from two-dimensional R objects shortly. \end{frame} \begin{frame}[fragile] \frametitle{2.3.2 Extracting elements from vectors} Another useful way of extracting elements from vectors is to use logical vectors. \vsp Recall that the logical data type in R can take only the values \texttt{TRUE} and \texttt{FALSE}. Hence: \begin{lstlisting} A = TRUE B = FALSE !A #%* \textcolor{ForestGreen}{Read as ``Not A"}*) [1] FALSE !B [1] TRUE A & B #%* \textcolor{ForestGreen}{Need ``both A \textbf{and} B"}*) [1] FALSE A | B #%* \textcolor{ForestGreen}{Need ``either A \textbf{or} B"}*) [1] TRUE \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.3.2 Extracting elements from vectors} Then we could declare a vector of such logicals to help us extract from another vector: \begin{lstlisting} (logic1 = c(TRUE, FALSE, TRUE, FALSE)) [1] TRUE FALSE TRUE FALSE (vec1 = seq(3, 6)) [1] 3 4 5 6 vec1[logic1] #%* \textcolor{ForestGreen}{Will extract the elements in vec1 where logic1 is true}*) [1] 3 5 \end{lstlisting} \pause We can also sum a logical vector. For example: \begin{lstlisting} sum(logic1) [1] 2 sum(vec1[logic1]) [1] 8 \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{Example 2.1} If you've not already got RStudio installed on the machine you're using, \textbf{you must do so now} by following the advice on the Canvas Community or our student resources webpage: {\small \begin{center} \url{http://www.ncl.ac.uk/maths-physics/students/resources/} \end{center} } From this link, click on the link for our ``Computing page", log-in using your University details, and then click on the link for ``R". \begin{itemize} \item Here you will find some basic resources aimed at students who have never used R before \item Scroll down, and you'll find information about how to install R on your home PC \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{Example 2.1} \includegraphics[width=1cm]{ironman.jpg} In RStudio, commit the following lines of code. Both lines read in files stored on a webpage. \textbf{Do not copy-and-paste this code -- you must type it in carefully!} \begin{lstlisting} x = scan("https://www.mas.ncl.ac.uk/~nlf8/teaching/mas2906/x.txt") y = as.logical(scan("https://www.mas.ncl.ac.uk/~nlf8/teaching/mas2906/y.txt", what="logical")) \end{lstlisting} \vspace{-1cm} {\tiny \begin{center} \begin{tabular}{l l} Questions about x &Questions about y \\ \hline (a) What type of data is in the vector x? &(a) What type of data is in the vector y? \\ (b) How many values are in the vector x? &(b) What is the first value of y?\\ (c) What is the smallest value in the vector x? &(c) What is the last value of y?\\ (d) What is the value of x in position 3500? &(d) How many times does TRUE occur? \\ (e) What is the 97th smallest value in x? &(e) When is the 472nd occurrence of TRUE? \\ (f) Find $\sum_{i=5612}^{8664}\text{x}_{i}$. &(f) How many occurrences of TRUE are there \\ & \textcolor{white}{LL} between positions 504 and 593 (inclusive)? \\ \hline \end{tabular} \end{center} } \end{frame} \subsection{Dataframes}\label{sec:dataframes} \begin{frame} \frametitle{2.4 Dataframes} A dataframe is a special kind of object in R. \vsp We use dataframes for storing and managing datasets that have a rectangular structure. \vsp Typically the rows correspond to \textcolor{blue}{\bf cases} and the columns \textcolor{blue}{\bf variables}. \vsp The crucial difference between a dataframe and a matrix is that all values in a matrix must be of the same type, whereas in a dataframe we can have responses for cases on variables of different types. \end{frame} \begin{frame} \frametitle{2.4 Dataframes} For example, we might have survey results on 10 individuals, who have each been asked 5 questions. \begin{itemize} \item This would give a dataframe with 10 rows and 5 columns \item We call this a ``\textcolor{blue}{\bf 10 by 5 dataframe}" (generally, we refer to $r \times c$ dataframes) \item The variable types need not all be the same: \begin{itemize} \item [--] The first column could be a vector of doubles recording the individuals' height in metres \item [--] The second column could be a vector of logicals taking the value \texttt{TRUE} if they like Marmite \item [--] The third column could be a vector of characters taking the name of their favourite book, etc. \end{itemize} \end{itemize} \end{frame} \subsubsection{Creating and querying your own dataframe} \begin{frame}[fragile] \frametitle{2.4.1 Creating and querying your own dataframe} You can construct your own dataframe using R's \texttt{data.frame} command. For example: \begin{lstlisting} height_m = c(1.81, 1.75, 1.58) weight_kg = c(NA, NA, 70.3) marmite = c(TRUE, FALSE, FALSE) book = c("Midnight in Chernobyl", "To Kill a Mockingbird", "Life of Pi") (small.survey = data.frame(height_m, weight_kg, marmite, book)) height_m weight_kg marmite book 1 1.81 NA TRUE Midnight in Chernobyl 2 1.75 NA FALSE To Kill a Mockingbird 3 1.58 70.3 FALSE Life of Pi \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.4.1 Creating and querying your own dataframe} If you had a large dataframe and wanted to quickly know its size, you could use the \texttt{dim} command. For example: \begin{lstlisting} dim(small.survey) [1] 3 4 \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.4.1 Creating and querying your own dataframe} We can extract elements from a dataframe using square brackets (just as we demonstrated in Section 2.3.2 for vector extraction). \vsp If the columns have labels, we can also use the \$ command. \begin{lstlisting} small.survey[2, 4] #%* \textcolor{ForestGreen}{Extract row 2, column 4}*) [1] To Kill a Mockingbird small.survey[ , 4] #%* \textcolor{ForestGreen}{Extract all rows from column 4}*) [1] Midnight in Chernobyl To Kill a Mockingbird Life of Pi small.survey$book #%* \textcolor{ForestGreen}{Does exactly the same as line 4!}*) [1] Midnight in Chernobyl To Kill a Mockingbird Life of Pi small.survey[3, ] #%* \textcolor{ForestGreen}{Extract all columns from row 3}*) height_m weight_kg marmite book 3 1.58 70.3 FALSE Life of Pi \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{2.4.1 Creating and querying your own dataframe} Though used less, the \texttt{subset} command can be applied to extract information from selected columns only. For example: \begin{lstlisting} (small.survey2 = subset(small.survey, select = 3:4)) marmite book 1 TRUE Midnight in Chernobyl 2 FALSE To Kill a Mockingbird 3 FALSE Life of Pi \end{lstlisting} \pause And what if you wanted the height and favourite book of respondents who didn't like Marmite? \begin{lstlisting} subset(small.survey,marmite == FALSE,select=c(1,4)) height_m book 2 1.75 To Kill a Mockingbird 3 1.58 Life of Pi \end{lstlisting} \end{frame} \subsubsection{The movies dataframe}\label{sec:moviesdataframe} \begin{frame}[fragile] \frametitle{2.4.2 The movies dataframe} The Movies dataset is available from within the \texttt{ncldata} package. \pause \vsp Recall that you can access this in the following way: \begin{lstlisting} install.packages("ncldata", repos="http://R-Forge.R-Project.org", type="source") library(ncldata) data(movies) \end{lstlisting} \pause We will now explore the movies dataset in RStudio. \end{frame} \begin{frame} \frametitle{Example 2.2} \begin{tabular}{c l} \multirow{3}{*}{\includegraphics[width=1cm]{ironman.jpg}} &Follow the steps at the start of Section 2.4.2 to load the \\ &\texttt{ncldata} package and to read in the \texttt{movies} dataframe. \\ & Then answer the following questions. \end{tabular} \begin{itemize} \item [(a)] How many people submitted votes for the films ``\textit{10 Things I Hate About You}", ``\textit{Midnight in the Garden of Good and Evil}" and ``\textit{Zoolander}"? \item [(b)] How many films are classified as ``Action"? How many are classified as ``Romance"? How many ``Short" films are there in the database? \item [(c)] How many films were made between the Years 1970 and 2000 (inclusive)? \item [(d)] How many movies are longer than two hours in duration? \item [(e)] How many movies made since the year 2005 (not including 2005) are longer than three hours in duration? \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{Example 2.3} \begin{tabular}{c l} \multirow{3}{*}{\includegraphics[width=1cm]{ironman.jpg}} & Load and explore the \texttt{music} dataset from the \texttt{ncldata} \\ &package by typing the following lines of code:\footnote{This will only work if you have correctly installed the \texttt{ncldata} package and read the package into your working RStudio space using the \texttt{library} command; see example 2.2 and Section 2.4.2)} \\ \end{tabular} \begin{lstlisting} data(music) #%* \textcolor{ForestGreen}{What's this dataset all about? Let's find out}*) ?music \end{lstlisting} \end{frame} \begin{frame} \frametitle{Example 2.3} \begin{itemize} \item [(a)] How many tracks do we have data for, and on how many variables? \item [(b)] How many tracks do we have from \textit{The Beatles} in this dataframe? What about \textit{Lady GaGa}? \item [(c)] How many tracks have a tempo of at least 130 beats per minute? \item [(d)] Of the tracks with a tempo of at least 130 beats per minute, how many have been released since the year 2000 (including the year 2000)? \item [(e)] What is the most popular key for tracks in this database? \end{itemize} \end{frame} \begin{frame} \frametitle{All done! Now what?} \begin{itemize} \item Computer practical -- \textcolor{blue}{\bf tomorrow (1-3)} or \textcolor{blue}{\bf Wednesday (11-1)} \vspl \item Complete \textcolor{blue}{\bf Examples 2.1 and 2.2} in your own time (video solutions to follow) \pause \vspl \item Next lecture: \textcolor{blue}{\bf Monday 30th Setember, 2-3pm} \end{itemize} \end{frame} \end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OLD STUFF \begin{frame}\frametitle{MAS2602 Arrangements} \begin{itemize} \item Statistics classes run in teaching weeks 5--9 \item Lecturer is Dr. Lee Fawcett (\textcolor{blue}{\texttt{lee.fawcett@ncl.ac.uk}}) \item Schedule: 4 lectures, 4 practicals, 1 revision class, 1 class test \item Office hours: Tuesdays 3-4; room 2.07 Herschel \\ (also MAS3902: Tuesdays 2-3 and Thursdays 3-4) \end{itemize} \end{frame} \begin{frame} \frametitle{Schedule} %\vspace{-1.1cm} {\small \begin{center} \begin{tabular}{llll} Week 5 &Mon 28 Oct &11--12 &Lecture (Herschel LT1) \\ &Fri 1 Nov &10--12 &Practical (Herschel PC) \\ & & & \\ Week 6 &Mon 4 Nov &11--12 &Lecture (Herschel LT1) \\ &Wed 6 Nov &11--1 &Practical (Herschel PC) \\ & & & \\ Week 7 &Mon 11 Nov &11--12 &Lecture (Herschel LT1) \\ &Thu 14 Nov &11--1 &Practical (Herschel PC) \\ & & & \\ Week 8 &Mon 18 Nov &11--12 &Lecture (Herschel LT1) \\ &Thu 21 Nov &11--1 &Practical (Herschel PC) \\ &\textcolor{red}{\it Thu 21 Nov} &\textcolor{red}{\it 3pm} &\textcolor{red}{\it Assignment due} \\ &Fri 22 Nov &1--2 &Revision (Herschel LT2) \\ & & & \\ Week 9 &\textcolor{red}{\it Tue 26 Nov} &\textcolor{red}{\it 9--10} &\textcolor{red}{\it Test (Herschel PC)} \\ \end{tabular} \end{center} } \end{frame} \begin{frame} \frametitle{Assessment} \textbf{Assignment (a.k.a.\ ``mini-project'')} \begin{itemize} \item Due \textcolor{red}{\bf Thursday 21 November, 3pm} \item Worth $10\%$ of the module marks \end{itemize} \textbf{Class test} \begin{itemize} \item \textcolor{red}{\bf Tuesday 26 November, 9am} -- one hour long \item Worth $30\%$ of module marks \item Open-book test \item You will write one or two short R programs during the test \end{itemize} \end{frame} \begin{frame} \frametitle{Help and Support} Help available: \begin{itemize} \item Office hours \item Demonstrators in practical sessions \item Books -- see recommendations in booklet \item Email Lee, or just pop in! \item Blackboard and dedicated webpage \end{itemize} \end{frame} \begin{frame} \frametitle{Late work policy} \begin{center} In this module, deadline extensions can be requested for the final project (by means of submitting a PEC form), and work submitted within 7 days of the deadline without good reason will be marked for reduced credit. This module also contains tests worth more than $10\%$ for which rescheduling can be requested (by means of submitting a PEC form). There are mini-projects (worth $10\%$ each) for which it is not possible to extend deadlines and for which no late work can be accepted. For details of the policy (including procedures in the event of illness etc.) please look at the School web site: http://www.ncl.ac.uk/maths/students/teaching/homework/ For problems with deadlines, speak to your personal tutor and prepare a PEC form \end{center} \end{frame} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{frame} \begin{Large} \begin{center} {\bf Lecture 1: Introduction and Simulation of Random Variables} \end{center} \end{Large} \end{frame} \begin{frame}\frametitle{Introduction} In this part of the module we will do statistics with R: \begin{itemize} \item R is the foremost tool in modern computational statistics \item Using R teaches general concepts in programming \item It can be used to illustrate mathematical ideas in probability and statistics \end{itemize} \pause \textbf{Today's lecture:} simulating random variables \begin{enumerate} \item Simulating random variables seen in MAS1604 \item Using this to simulate more complicated probability models \end{enumerate} \pause \textbf{Friday's practical:} \begin{enumerate} \item Revision of R from MAS1608 \item Putting today's material into practice \end{enumerate} \end{frame} \begin{frame}\frametitle{The binomial distribution} If $X\sim Bin(n,p)$ then $X$ has PMF\ (probability mass function) given by \textcolor{BrickRed}{ \begin{equation*} p_X(k) = \binom{n}{k}p^k(1-p)^{n-k},\quad\text{for\ }k=0,1,\ldots,n. \end{equation*}} There is no closed formula for the CDF (cumulative distribution function). \pause R commands: \begin{itemize} \item \texttt{dbinom} -- calculate PMF \item \texttt{pbinom} -- calculate CDF \item \texttt{rbinom} -- generate random sample \end{itemize} \end{frame} \begin{frame}[fragile]\frametitle{R commands for the binomial distribution} \begin{lstlisting} dbinom(5, 10, 0.7) #%* \textcolor{ForestGreen}{$\pr{X=5}, X\sim Bin(10,0.7)$} *) pbinom(4, 10, 0.7) #%* \textcolor{ForestGreen}{$\pr{X\leq 4}, X\sim Bin(10,0.7)$} *) rbinom(50, 10, 0.7) #%* \textcolor{ForestGreen}{Sample a $Bin(10,0.7)$ distribution 50 times} *) \end{lstlisting} \end{frame} \begin{frame}[fragile]\frametitle{Creating a bar plot} \begin{lstlisting} x = rbinom(50, 10, 0.7) plot(table(x), xlim = c(0,10), xaxt = 'n', xlab = 'x', ylab = 'frequency') axis(1, at = seq(0, 10)) \end{lstlisting} \vspace{-1.5cm} \begin{center} \includegraphics[width=8cm]{../notes/fig1-1} \end{center} \end{frame} \begin{frame}\frametitle{The geometric distribution} If $Y\sim Geom(p)$ then $Y$ has PMF and CDF given by \textcolor{BrickRed}{ \begin{align*} p_Y(k) &= (1-p)^{k-1}p,\\ F_Y(k) &= 1-(1-p)^k,\quad\text{for\ }k=1,2,\ldots. \end{align*}} Note that $Y$ takes values in $1,2,3,\ldots$: \begin{itemize} \item R uses a slightly different definition \item We use the definition that $Y$ is the number of Bernoulli trials with up to and including first success \item R counts number of trials up to, but \hi{not} including the first success, so in R geometric random variables take values $0,1,2,\ldots$. \end{itemize} \end{frame} \begin{frame}[fragile]\frametitle{R commands for the geometric distribution} Adjust the arguments to account for different definition: \begin{lstlisting} dgeom(4, 0.2) #%* \textcolor{ForestGreen}{$\pr{Y=5}, Y\sim Geom(0.2)$} *) pgeom(2, 0.2) #%* \textcolor{ForestGreen}{$\pr{Y\leq 3}, Y\sim Geom(0.2)$} *) 1 + rgeom(100, 0.2) #%* \textcolor{ForestGreen}{Sample a $Geom(0.2)$ distribution 100 times} *) \end{lstlisting} Here's a function to replace \texttt{dgeom} with our definition of the geometric distribution: \begin{lstlisting} mydgeom = function(x, p) { dgeom(x-1, p)} \end{lstlisting} \end{frame} \begin{frame}\frametitle{The Poisson distribution} If $Z\sim Po(\lambda)$ then $Z$ has PMF given by \textcolor{BrickRed}{ \begin{equation*} p_Z(k) = \frac{\lambda^k}{k!}e^{-\lambda},\quad\text{for\ }k=0,1,2,\ldots. \end{equation*}} There is no closed formula for the CDF (cumulative distribution function). R commands: \begin{itemize} \item \texttt{dpois} -- calculate PMF \item \texttt{ppois} -- calculate CDF \item \texttt{rpois} -- generate random sample \end{itemize} \end{frame} \begin{frame}[fragile]\frametitle{R commands for the Poisson distribution} \begin{lstlisting} dpois(5, 3.5) #%* \textcolor{ForestGreen}{$\pr{Z=5}, Z\sim Po(3.5)$} *) ppois(2, 3.5) #%* \textcolor{ForestGreen}{$\pr{Z\leq 2}, Z\sim Po(3.5)$} *) rpois(100, 3.5) #%* \textcolor{ForestGreen}{Sample a $Po(3.5)$ distribution $100$ times} *) \end{lstlisting} \vspace{-0.5cm} \end{frame} \begin{frame}\frametitle{Summary} \begin{center} \begin{tabular}{c|ccc} Distribution & Binomial & Poisson & Geometric\ $\skull$\\ \hline PMF & \texttt{dbinom(...)} & \texttt{dpois(...)} & \texttt{dgeom(...)} \\ CDF & \texttt{pbinom(...)} & \texttt{ppois(...)} & \texttt{pgeom(...)} \\ sample & \texttt{rbinom(...)} & \texttt{rpois(...)} & \texttt{rgeom(...)} \\ \end{tabular} \end{center} \end{frame} \begin{frame}\frametitle{Continuous random variables} R has functions for the uniform, exponential and normal distributions: \begin{center} \begin{tabular}{c|ccc} Distribution & Uniform & Exponential & Normal\\ \hline PDF & \texttt{dunif(...)} & \texttt{dexp(...)} & \texttt{dnorm(...)} \\ CDF & \texttt{punif(...)} & \texttt{pexp(...)} & \texttt{pnorm(...)} \\ quantile & \texttt{qunif(...)} & \texttt{qexp(...)} & \texttt{qnorm(...)} \\ sample & \texttt{runif(...)} & \texttt{rexp(...)} & \texttt{rnorm(...)} \\ \end{tabular} \end{center} \end{frame} \begin{frame}[fragile]\frametitle{Continuous random variables -- R examples} \begin{lstlisting} dunif(3, 2, 5) #%* \textcolor{ForestGreen}{$f_{X}(3), X\sim U(2, 5)$} *) pexp(1.5, 5) #%* \textcolor{ForestGreen}{$F_{Y}(1.5), Y\sim Exp(5)$} *) rnorm(10, 0, 2) #%* \textcolor{ForestGreen}{Sample a $N(0, 2^2)$ distribution 10 times} *) \end{lstlisting} For the standard uniform ($U(0,1)$) and standard normal ($N(0,1)$) distributions you don't need to provide the parameters $a=0$, $b=1$ and $\mu=0$, $\sigma=1$ respectively. For example: \begin{lstlisting} runif(20) #%* \textcolor{ForestGreen}{Samples a $U(0,1)$ distribution 20 times} *) pnorm(1.96) #%* \textcolor{ForestGreen}{$\pr{Z<1.96}, Z\sim N(0,1)$} *) [1] 0.9750021 \end{lstlisting} \end{frame} \begin{frame}[fragile]\frametitle{Quantiles} The quantile functions \texttt{qunif, qexp, qnorm} solve equations like \begin{equation*} F_X(\alpha) = p \end{equation*} for $\alpha$ given a probability $p$. For example: \begin{lstlisting} qnorm(0.9750021) [1] 1.96 \end{lstlisting} \end{frame} \begin{frame}\frametitle{Quantiles} \begin{example} Suppose annual maximum wave heights observed off the coast at a flood-prone town are assumed Normally distributed, with mean 2 metres and standard deviation 0.5 metres. \begin{itemize} \item [(a)] Write down the R command to find the probability that the largest wave height next year will exceed 3.25 metres. \item [(b)] Write down the R command to estimate the height of a new sea wall such that we might expect the town to be flooded, on average, once per century. Why might our modelling assumption be invalid? \end{itemize} \end{example} \end{frame} \frame{ \frametitle{Example 1.1: Solution (a)} \includegraphics[width=\textwidth, height=\textheight]{seasurge} } \begin{frame}[fragile] \frametitle{Example 1.1: Solution (a)} Could work out the old-fashioned way: \pause \begin{eqnarray*} \pr{X>3.25} &=&\pr{Z>\frac{3.25-2}{0.5}} \\ & & \\ &=& \pr{Z>2.5} \\ & & \\ &=& 1-\pr{Z \leq 2.5} = 1-0.994 = 0.006. \end{eqnarray*} \pause Or could just use R: \begin{lstlisting} 1-pnorm(3.25, 2, 0.5) [1] 0.006209665 \end{lstlisting} \end{frame} \begin{frame}[fragile] \frametitle{Example 1.1: Solution (b)} \includegraphics[width=\textwidth, height=\textheight]{seasurge2} \end{frame} \begin{frame}[fragile] \frametitle{Example 1.1: Solution (b)} Similarly: \begin{lstlisting} qnorm(0.99, 2, 0.5) [1] 3.163174 \end{lstlisting} \end{frame} \begin{frame}\frametitle{A more advanced model} \begin{itemize} \item Number of arrivals per day at an IT help-desk is modelled using a Poisson distribution \item Mean of the Poisson distribution might vary from day-to-day \item Suppose the number of arrivals $X\sim Po(\Lambda)$ \item $\Lambda$ is itself a random variable, with $\Lambda\sim Exp(c)$ for a constant $c=0.05$ \end{itemize} What can we say about the distribution of $X$? \begin{enumerate} \item What are the expectation and variance of $X$? \item What is $\pr{X>30}$? \end{enumerate} \end{frame} \begin{frame}[fragile]\frametitle{Arrival model -- R code} \begin{lstlisting} arrivals = function(n, c = 0.05) { x = vector(mode = 'numeric',length = n) for (i in 1:n) { lambda = rexp(1, rate = c) x[i] = rpois(1, lambda) } return(x) } \end{lstlisting} \end{frame} \begin{frame}[fragile]\frametitle{Arrival model -- R code} Bar plot: \begin{lstlisting} x = arrivals(500) plot(table(x), xlim = c(0, 50), xaxt = 'n', xlab = 'x', ylab = 'frequency') axis(1, at = seq(0, 50, 5)) \end{lstlisting} \vspace{-1.5cm} \begin{center} \includegraphics[width=8cm]{../notes/fig1-2} \end{center} \end{frame} \begin{frame}\frametitle{Using simulated samples} \begin{itemize} \item The \hi{sample mean} is an approximation to the distribution mean \item The same applies for the sample variance \item To calculate $\pr{X> 30}$ we count the proportion of times this occurs in the sample \end{itemize} We expect the approximation to improve \hi{as we increase the sample size} \end{frame} \begin{frame}[fragile]\frametitle{Using simulated samples} \begin{lstlisting} mean(x) [1] 19.618 var(x) [1] 382.1764 sum(x>30)/500 [1] 0.202 \end{lstlisting} \end{frame} \begin{frame}\frametitle{Mixtures of normal distributions} Suppose \begin{itemize} \item $\mu_1,\mu_2,\sigma_1>0,\sigma_2>0$ are fixed constants, \item $w_1,w_2$ are positive constants with $w_1+w_2=1$. \end{itemize} Consider the following function: \begin{equation*} f(x) = w_1f_1(x)+w_2f_2(x) \end{equation*} where $f_1$ and $f_2$ are the density functions for $Z_1\sim N(\mu_1,\sigma_1^2)$ and $Z_2\sim N(\mu_2,\sigma_2^2)$ respectively. Check that $f(x)$ represents a valid probability density function. \end{frame} \begin{frame}\frametitle{Mixtures of normal distributions} \begin{answers} First note that $f(x) \geq 0$ everywhere. Also \begin{eqnarray*} \int_{-\infty}^{\infty} f(x) dx &=& w_{1}\int_{-\infty}^{\infty}f_{1}(x) dx + w_{2} \int_{-\infty}^{\infty}f_{2}(x)dx \pause \\ & & \\ &=& w_1 + w_2 = 1. \end{eqnarray*} \end{answers} \end{frame} \begin{frame}\frametitle{Mixtures of normal distributions} If a random variable $X$ has PDF corresponding to $f(x)$ we say it is a \hi{mixture} of normal distributions $N(\mu_1,\sigma_1^2)$ and $N(\mu_2,\sigma_2^2)$ with weights $w_1,w_2$. Note that $X$ is \hi{not} the sum of two normal random variables i.e. \begin{equation*} X \neq w_1Z_1+w_2Z_2\quad\text{where\ }Z_i\sim N(\mu_i,\sigma_i^2),\ i=1,2. \end{equation*} \end{frame} \begin{frame}\frametitle{Example} For example: suppose $\mu_1=3, \sigma_1=1$ and $\mu_2=6, \sigma_2=2$ with $w_1=w_2=1/2$. \vspace{-1.7cm} \begin{center} \includegraphics[width=9.5cm]{../notes/fig1-3} \end{center} \end{frame} \begin{frame}\frametitle{Sampling from a mixture distribution} \begin{enumerate} \item Sample a random variable $J\in\{1,2\}$ such that $\pr{J=1}=w_1$ and $\pr{J=2}=w_2$. \item The random variable $J$ tells you which \hi{component} of the mixture to sample $X$ from. \item If $J=1$ then sample $X$ from $N(\mu_1,\sigma_1^2)$, but if $J=2$ then sample $X$ from $N(\mu_2,\sigma_2^2)$. \end{enumerate} \end{frame} \begin{frame}[fragile]\frametitle{Sampling from a mixture distribution} \begin{lstlisting} normal.mixture = function(n, mu1, sig1, w1, mu2, sig2, w2) { p = c(w1, w2) x = vector(mode = 'numeric', length = n) for (i in 1:n) { j = sample(c(1, 2), 1, prob = p) if (j==1) { x[i] = rnorm(1, mu1, sig1) } else { x[i] = rnorm(1, mu2, sig2) } } return(x) } \end{lstlisting} \end{frame} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{frame} \begin{Large} \begin{center} {\bf Lecture 2: Inverse CDF method} \end{center} \end{Large} \end{frame} \begin{frame}\frametitle{Today's lecture} In this lecture we look at the \textcolor{purple}{\bf inverse CDF method} for simulating from continuous random variables. \end{frame} \begin{frame}\frametitle{The inverse CDF method} This is a method for simulating \hi{univariate continuous} random variables \begin{itemize} \item Let $U\sim U(0,1)$ \item Suppose $F(x)$ is a well-defined CDF which is invertible \item Then the random variable $X=F^{-1}(U)$ has CDF $F(x)$ \end{itemize} \end{frame} \begin{frame} \frametitle{The inverse CDF method} \begin{answers} \[ X = F^{-1}(U). \pause \] So \begin{eqnarray*} \pr{X \leq x} &=& \pr{F^{-1}(U) \leq x} \\ \pause & & \\ &=&\pr{F(F^{-1}(U)) \leq F(x)} \\ \pause & & \\ & = &\pr{U \leq F(x)}\\ \pause & &\\ &=&F(x), \end{eqnarray*} since $0 \leq F(x) \leq 1$ (see sketch on board). \end{answers} \end{frame} \begin{frame}\frametitle{The inverse CDF method} \begin{example} Suppose we obtain two observations from a $U(0,1)$ distribution: $0.1$ and $0.85$. Use these values to obtain two observations from $X\sim Exp(2)$. \end{example} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} \noindent If $X \sim Exp(2)$ then $F_{X}(x) = 1-e^{-2x}$. \pause \vspace{0.3cm} If $y = 1-e^{-2x}$, then \begin{eqnarray*} e^{-2x} &=& 1-y \\ & & \\ -2x &=& \log (1-y) \\ & & \\ x &=& -\frac{1}{2} \log(1-y). \end{eqnarray*} \pause Plug in $y=0.1$ and $y=0.85$; gives realized values of $0.0527$ and $0.949$. \end{answers} \end{frame} \begin{frame}\frametitle{Example} \begin{itemize} \item Patients arriving at a doctor's surgery have appointment times $15$ minutes apart. \item The random variable $L$ measures how late a patient is for his/her appointment, with negative values denoting early arrivals. \end{itemize} The PDF of $L$ is \begin{equation*} f_L(\ell) = \begin{cases} 0.0025(\ell+20),\ & -20\leq \ell \leq 0,\\ 0.05e^{-\ell/10},& \ell>0. \end{cases} \end{equation*} \end{frame} \begin{frame}\frametitle{Questions} \begin{enumerate}[a.] \item Sketch $f_L(\ell)$. \item Show that the CDF of $L$ is \begin{equation*} F_L(\ell) = \begin{cases} 0,&\ell<-20,\\ 0.00125(\ell+20)^2, & -20\leq \ell \leq 0,\\ 1-\frac{e^{-\ell/10}}{2}, & \ell>0. \end{cases} \end{equation*} \item Suppose that $0.205$ and $0.713$ are two realized values of $U\sim U(0,1)$. Show that the corresponding realized values of $L$ are $-7.19$ and $5.55$ under the inverse CDF simulation method. \end{enumerate} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} \textbf{First, the PDF:} \vspace{-1.4cm} \begin{center} \includegraphics[width=8.25cm]{Ex24.pdf} \end{center} \end{answers} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} \noindent \textbf{Now to find the CDF:} \pause $L$ lies between $-20$ and $+\infty$, so $F_{L}(\ell) = 0$ when $\ell < -20$. \pause \vspace{0.3cm} \noindent For $-20 \leq \ell \leq 0$: \pause \begin{eqnarray*} F_{L}(\ell) &=& \int_{-20}^{\ell} f_{L}(s) ds \pause \\ & & \\ &=& \int_{-20}^{\ell} 0.0025(s+20) ds \\ \pause & & \\ &=&\left[\frac{0.0025}{2}(s+20)^{2}\right]^{\ell}_{-20} \pause\\ & & \\ &=& 0.00125(\ell+20)^{2} \quad \textsf{for $-20 \leq \ell \leq 0$}. \end{eqnarray*} Also, $F_{L}(0)=0.5$. \end{answers} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} Now \begin{eqnarray*} F_{L}(\ell) &=& \int_{-20}^{0}f_{L}(s) ds + \int_{0}^{\ell}f_{L}(s) ds \quad \textsf{for $\ell>0$} \pause \\ && \\ &=& 0.5 + \int_{0}^{\ell}0.05e^{-s/10} ds \pause \\ && \\ &=& 0.5 + \left[0.05 \times -10 e^{-s/10}\right]_{0}^{\ell} \pause \\ && \\ &=&0.5 +0.5 - 0.5 e^{-\ell/10} \pause \\ & & \\ &=& 1- \frac{e^{-\ell/10}}{2}, \quad \textsf{as required.} \end{eqnarray*} \end{answers} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} \noindent \textbf{Now for some realized values:} \pause \vspace{0.3cm} $u=0.205 <0.5$, so realized value $\ell$ lies between $-20$ and $0$. Thus \pause \[ F_{L}(\ell) = 0.00125(\ell+20)^2. \] \pause Find the inverse, and evaluate at $u=0.205$: \pause \begin{eqnarray*} y &=&0.00125(\ell+20)^2, \quad \textsf{so} \pause\\ \ell &=&\pm\sqrt{\frac{y}{0.00125}}-20, \end{eqnarray*}\pause giving \[ \ell = \sqrt{\frac{0.205}{0.00125}}-20 = -7.2. \] \end{answers} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} Similarly, $u=0.713 >0.5$ so realized value $>0$. Thus\pause \[ F_{L}(\ell) = 1- \frac{e^{-\ell/10}}{2}. \]\pause Find the inverse, and evaluate at $u=0.713$: \pause \begin{eqnarray*} y &=& 1-\frac{e^{-\ell/10}}{2} \quad \textsf{so} \pause \\ \ell &=& -10\log(2-2y), \end{eqnarray*} giving \[ \ell = -10\log(2-2\times 0.713) = 5.55. \] \end{answers} \end{frame} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{frame} \begin{Large} \begin{center} {\bf Lecture 3: Transformation of \\ Random Variables} \end{center} \end{Large} \end{frame} \begin{frame} \frametitle{Recap Quiz} \begin{itemize} \item [1.] Write down R code that simulates 100 observations from the $U(0,1)$ distribution and stores the output in the vector \texttt{U}. \item [2.] Suppose $X$ has CDF \[ F_{X}(x) = 1- \frac{1}{x}, \quad x\geq 1. \] Find $F^{-1}_{X}(x)$. \item [3. ] Write an R function called \texttt{inv.cdf} to return the inverse of $F_{X}(x)$. \item [4. ] Write down a single line of R code to generate 100 realizations of the random variable $X$. \item [5.] Write down a single line of R code to produce a histogram of your realizations of $X$, coloured yellow and with a suitable title. \end{itemize} \end{frame} \begin{frame}\frametitle{Transforming a normal random variable} \begin{itemize} \item Consider how we might model a \hi{continuous positive} random quantity e.g.\ blood pressure or height \item The Normal distribution is not suitable! \end{itemize} One way to obtain a positive random variable is to define \begin{equation*} Y = \exp(X),\quad\text{where\ }X\sim N(\mu,\sigma^2). \end{equation*} Then $Y$ is a random variable which is a \hi{transformation} of $X$. \end{frame} \begin{frame}[fragile]\frametitle{Simulation to explore transformation} \begin{lstlisting} y = exp(rnorm(250, mean = 1.5, sd = 0.3)) hist(y, xlab = 'Y', ylab = 'frequency') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{../notes/fig3-1} \end{center} \end{frame} \begin{frame}\frametitle{Calculating PDFs} \begin{itemize} \item This looks like a useful distribution -- what is its PDF? \item $Y=\exp(X)$ so $X = \log(Y)$ \item So $f_Y(y)=f_X(\log(y))$ is a good guess, but is \hi{NOT} the right answer! \end{itemize} \vspace{1cm} \begin{example} If $X\sim U(0,1)$, what is the PDF of $Y=X^3$? \end{example} \end{frame} \begin{frame} \frametitle{Solution} \begin{answers} \begin{eqnarray*} F_{Y}(y) &=&\pr{Y \leq y} \pause\\ &=&\pr{X^{3} \leq y} \pause\\ &=& \pr{X \leq y^{1/3}} \\ &=& F_{X}(y^{1/3}). \pause \end{eqnarray*} Therefore \begin{eqnarray*} f_{Y}(y) &=& \frac{d}{dy} F_{X}(y^{1/3}) \pause \\ &=&f_{X}(y^{1/3}) \times \frac{1}{3}y^{-2/3} \quad \textsf{by chain rule} \pause \\ &=& \frac{1}{3}y^{-2/3}, \quad 0\leq y \leq 1 \textsf{ (and zero otherwise)}. \end{eqnarray*} \end{answers} \end{frame} \begin{frame}[fragile]\frametitle{Checking via simulation} To confirm this result, we can simulate $Y$: \begin{lstlisting} y = (runif(500))^3 hist(y, xlab = 'Y', ylab = 'frequency') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{../notes/fig3-2} \end{center} \end{frame} \begin{frame}\frametitle{Further examples} \begin{example} If $X\sim N(0,1)$, what is the PDF of $Y=X^2$? \end{example} \pause \begin{answers} \[ F_{Y}(y) =\pr{Y \leq y} =\pr{X \leq y^{1/2}} = F_{X}(y^{1/2}). \] \pause Then \begin{eqnarray*} f_{Y}(y) &=& \frac{d}{dy} F_{X}(y^{1/2}) \pause \\ &=&f_{X}(y^{1/2}) \times \frac{1}{2}y^{-1/2} \quad \textsf{by chain rule} \pause \\ &=&\frac{1}{\sqrt{2\pi}}e^{-y/2} \times \frac{1}{2}y^{-1/2} \pause \\ &=&\frac{1}{2 \sqrt{2 \pi y}}e^{-y/2} \quad \textsf{for $y \geq 0$, zero otherwise.} \end{eqnarray*} \end{answers} \end{frame} \begin{frame}\frametitle{Further examples} \begin{example} If $X\sim N(\mu,\sigma^2)$, what is the PDF of $Y=\exp(X)$? \\ This is the example used at the start of the lecture. We say $Y$ has the \hi{log-normal} distribution and write $Y\sim LN(\mu,\sigma^2)$. \end{example} \pause \begin{answers} \begin{eqnarray*} F_{Y}(y) &=& F_{X}(\log y) \quad \textsf{by usual argument with CDF} \pause \\ f_{Y}(y) &=&\frac{d}{dy}F_{X}(\log y) \pause \\ &=&f_{X}(\log y) \times \frac{1}{y} \quad \textsf{for $y \geq 0$} \pause \\ &=&\frac{1}{\sigma y \sqrt{2 \pi}}\exp\left\{-\frac{1}{2}\left(\frac{(\log y - \mu)^{2}}{\sigma^{2}}\right)\right\}. \end{eqnarray*} \end{answers} \end{frame} \begin{frame}\frametitle{The univariate transformation method} If $X$ has pdf $f_X(x)$ and $h(x)$ is a monotonically increasing or monotonically decreasing function then \begin{itemize} \item $h$ is invertible and \item the PDF of $Y=h(X)$ is \begin{equation*} f_Y(y) = f_X(h^{-1}(y))\times\left| \frac{d\left( h^{-1}(y) \right)}{dy} \right| \end{equation*} when $y$ is in the range of $h(x)$ and zero otherwise. \end{itemize} This result allows direct calculation of $f_Y(y)$ and is called the \hi{univariate transformation method}. \begin{center} \textcolor{BrickRed}{\bf See whiteboard for Proof.} \end{center} \end{frame} \begin{frame}\frametitle{Some Examples} \begin{example} When $X\sim U(0,1)$, use the transformation method to find the PDF of $Y=X^3$. \end{example} \pause \begin{answers} \[ h(x) =x^{3}, \quad h^{-1}(y) = y^{1/3}. \] \pause \begin{eqnarray*} f_{Y}(y) &=& 1 \times \frac{d}{dy}y^{1/3} \pause \\ &=& \frac{1}{3}y^{-2/3}, \quad 0\leq y \leq 1. \end{eqnarray*} \end{answers} \end{frame} \begin{frame}\frametitle{Some examples} \begin{example} When $X\sim U(-1,1)$ and $Y=\exp(X)$, use the transformation method to find the PDF of $Y$. \end{example} \pause \begin{answers} \[ h(x) = e^{x}, \quad h^{-1}(y) = \log (y). \]\pause \begin{eqnarray*} f_{Y}(y) &=& \frac{1}{2} \times \abs[\Big]{\frac{d}{dy} \log (y)} \pause \\ &=& \frac{1}{2|y|} = \frac{1}{2y}, \quad e^{-1} \leq y \leq e^{1}, \textsf{ zero otherwise.} \end{eqnarray*} \end{answers} \end{frame} \begin{frame}\frametitle{Some examples} \begin{example} Use the transformation method to find the PDF of the log-normal distribution. \end{example} \pause \begin{answers} \[ h(x) = e^{x}, \quad h^{-1}(y) = \log (y). \] \pause \begin{eqnarray*} f_{Y}(y) &=& f_{X}(\log (y)) \times \frac{1}{y} \qquad \textsf{for $y>0$} \pause \\ &=& \frac{1}{y} \times \frac{1}{\sigma \sqrt{2 \pi}} \exp\left[-\frac{1}{2} \left(\frac{\log (y) - \mu}{\sigma}\right)^{2}\right]. \end{eqnarray*} \end{answers} \end{frame} \begin{frame}\frametitle{Cauchy example} \begin{example} Suppose $X\sim U(-\frac{\pi}{2}, \frac{\pi}{2})$ and $Y=\tan(X)$. \\ Use simulation to study the distribution of $Y$ and explicitly compute the PDF. \end{example} \end{frame} \begin{frame}[fragile]\frametitle{Solution} \begin{lstlisting} y = tan(runif(100, -pi/2, pi/2)) hist(y, xlab = 'y', ylab = 'frequency', main = '') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{../notes/fig3-3} \end{center} \end{frame} \begin{frame}\frametitle{Solution} \begin{answers} This distribution looks odd! Calculate the PDF: \pause \[ h(x) = \tan(x), \quad h^{-1}(y) = \tan^{-1}(y). \] \pause \begin{eqnarray*} f_{Y}(y) &=& \frac{1}{\pi} \frac{d}{dy}\left(\tan^{-1} y \right) \quad -\infty < y < \infty, \pause \\ &=&\frac{1}{\pi(1+y^{2})} \quad -\infty