\documentclass[12pt,a4paper]{article} \usepackage{amsthm,amsmath,amsfonts} \usepackage{latexsym} \usepackage{fancyhdr} \usepackage{enumerate} %\usepackage{marginnote} \usepackage{graphicx} \usepackage{skull} \usepackage{fancybox} \usepackage{hyperref} \usepackage[T1]{fontenc} %\usepackage{xcolor} %%From web for code highlighting:: \usepackage{listings} \usepackage{cmbright} \addtolength{\parskip}{1ex} \baselineskip = 6mm %\parskip = 2mm \textheight 23.6cm \textwidth 16.6cm \oddsidemargin 0.0cm \evensidemargin 0.0cm \topmargin -1.4cm %\headsep -.60in \setlength{\parindent}{0cm} \pagestyle{fancy} \lhead{} \chead{} \rhead{} \lfoot{} \cfoot{\thepage} \rfoot{} \rhead{} \newcommand{\pr}[2][]{\ensuremath{\mathrm{Pr}_{#1}\!\left(#2\right)}} \newenvironment{myhighlight}{\textit}{} \renewenvironment{emph}{\textbf}{} \newcommand{\Exp}[2][]{\ensuremath{\mathrm{Exp}\left(#2\right)}} \newcommand{\Bin}[2][]{\ensuremath{\mathrm{Bin}\left(#2\right)}} \newcommand{\Po}[2][]{\ensuremath{\mathrm{Po}\left(#2\right)}} \newcommand{\cov}[2][]{\ensuremath{\mathrm{Cov}\left(#2\right)}} \newcommand{\cor}[2][]{\ensuremath{\mathrm{Cor}\left(#2\right)}} \newcommand{\diag}[2][]{\ensuremath{\mathrm{diag}\left(#2\right)}} \newcommand{\given}{\ensuremath{\ |\ }} \newcommand{\stationary}{\ensuremath{\vec{\pi}}} \newcommand{\identity}{\ensuremath{\vec{I}}} \renewcommand{\vec}[1]{\underline{#1}} \newcommand{\remove}[1]{} %\newcommand{\cover}[2]{\hi{#1}} \newcommand{\cover}[2]{#2} \newcommand{\hi}[2][]{\textbf{#2}} \newcommand{\bigfont}{\fontsize{1.45cm}{1.5cm}\selectfont} \newtheoremstyle{exmp}% hnamei {3pt}% hSpace abovei {3pt}% hSpace belowi {}% hBody fonti {}% hIndent amounti1 {\bf}% hTheorem head fonti {:}% hPunctuation after theorem headi {.5em}% hSpace after theorem headi2 {}% hTheorem head spec (can be left empty, meaning `normal')i \theoremstyle{exmp} \newtheorem{example}{\textcolor{purple}{\textbf{Example}}}[section] %\def \var{{\textrm{Var}}} \newcommand{\E}[2][]{\ensuremath{\mathrm{E}\left[#2\right]}} \def \R{\mathbb{R}} \def \Z{\mathbb{Z}} \def \proof{\noindent{\bf Proof: }} \def \eproof{\qed\vskip .3cm} \def \union{\cup} \def \intersect{\cap} \newcommand{\ex}[2][]{\ensuremath{\mathrm{E}_{#1}\!\left[#2\right]}} \newcommand{\var}[2][]{\ensuremath{\mathrm{Var}_{#1}\left(#2\right)}} %\renewcommand{\det}[1][]{\ensuremath{\mathrm{det}\left(#1\right)}} \newtheorem{defn}{Definition} %\newtheorem*{example}{Example} \newtheorem{result}{Result}[section] \newtheorem{algorithm}{Algorithm}[section] \title{\includegraphics[width=0.9\textwidth]{crestBW.jpg} \\ \vspace{3cm} {\bigfont \bf MAS2906} \\ \vspace{1cm} {\bigfont Computational Probability and Statistics with R } \\ \vspace{3cm} } \date{\LARGE Semester 1, 2023/24} \author{\LARGE Dr. Lee Fawcett} \usepackage{Sweave} \usepackage{hyperref} \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} \definecolor{LEEgreen}{rgb}{0.04,0.667,0.137} \definecolor{LEEred}{rgb}{0.44,0.15,0.137} \definecolor{UBCblue}{rgb}{0.04706, 0.13725, 0.66667} % UBC Blue (primary) \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 } %\SweaveOpts{concordance=TRUE} \maketitle \pagenumbering{gobble} \setcounter{page}{-3} % to get contents page labelled (i) \newpage \textcolor{white}{Lee} \newpage \setcounter{page}{1} \pagenumbering{arabic} \lhead{MAS2906} \rhead{Chapter 1: General Housekeeping} %\rhead{Information for Statistics Classes (Weeks 5--9)} \section{General Housekeeping} \subsection{Lecturer information} The module leader for MAS2906 is Dr. Lee Fawcett. If you have any questions about the module he can be contacted via email at \texttt{lee.fawcett@newcastle.ac.uk}. Lee is based in room 2.07 on the second floor of the Herschel Building. There will be several other people on the MAS2906 team assisting with this module, primarily in the computer practical sessions. They're all really friendly PhD students -- give them a shout in class if you're stuck with anything! \subsection{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: Wednesday 10-11, Herschel LT2 \item Computer lab: Thursday 9-11, Herschel Computing Cluster \end{itemize} \item After enrichment week: \begin{itemize} \item Lectures: Tuesday 10-11, Herschel LT1; Wednesday 10-11, Herschel LT2 \item Computer lab: Thursday 9-11, Herschel Computing Cluster \end{itemize} \end{itemize} Full details can be found in the \href{https://ncl.instructure.com/courses/51393/pages/module-schedule}{Module Schedule} on Canvas. Occasionally, there will be additional tasks for you to complete during the week, and of course there will be the formal assessments to complete! \subsection{Assessment} MAS2906 has no written exam; you will be assessed via a Numbas test in January, and two homework assignments: \begin{itemize} \item [--] Problem-Solving Exercises 1 -- Numbas assignment (\textbf{15\%}): Covers basic R work from the first few weeks of the module. The assignment will be entirely completed within the Numbas application. Practice material will be provided. The assignment will be open for a two week period, with a deadline of \textcolor{red}{\bf 4pm, Friday 3 November 2023} \item [--] Problem-Solving Exercises 2 -- Mini Project (\textbf{15\%}): The project will cover some of the more theoretical aspects of the course, but will also require you to use some of your R skills from the first few weeks of the module. The project material will be released three weeks before the deadline of \textcolor{red}{\bf 4pm, Friday 8 December 2023}. Submission will be via Canvas upload \item [--] Numbas test (\textbf{70\%}): This test, during the formal University examination period, will cover a range of material from across the module. The test will be completed within the Numbas application, but will also require a Canvas upload element. Sample test questions will be provided, and dedicated time will be made for revision in class. The test will take place during the \textcolor{red}{\bf January assessment period} -- precise details will be released later in the semester. \end{itemize} \subsection{Late work policy} In this module, deadline extensions can be requested for the 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). For details of the policy (including procedures in the event of illness etc.) please consult the Mathematics, Statistics \& Physics Community pages on Canvas, under: Assessment Information, Late Work and Missed Assessments. \subsection{Course materials} All course materials will be posted to the dedicated MAS2906 Canvas page -- where you will find: \begin{itemize} \item A full copy of these lecture notes, with gaps \item A copy of the lecture notes with the gaps filled in, after the associated lectures/computer labs \item All slides used to support lectures/computer labs, after these classes have taken place \item All video recordings \item Solutions to assessed exercises \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} At various points in the module you will be directed to additional reading, should you want extra examples or just another source of information relating to the material presented in the notes. This will primarily be text books (many of which are available to view online), but perhaps other sources (e.g. journal publications, magazine articles etc.). \\ \\ Some useful books include: \begin{itemize} \item W.N. Venables and D.M. Smith: An Introduction To R. Free to download at: \begin{center} \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} In MAS2906 we will learn about computational tools for performing probability calculations and statistical analysis. We will use the R programming language. R is very widely used in universities, businesses and other organizations due to its many advantages over other similar software packages. It is available free of charge, it is `open source' and can be used for a very wide range of statistical and scientific applications. Libraries of R functions and data-structures have been created and are readily available on the web, enabling statisticians to perform a very wide range of tasks but with a minimum of additional programming. Learning to use R is important because: \begin{itemize} \item [(i)] It is the foremost tool in modern computational statistics \item [(ii)] It is a great way to learn general concepts in programming which can be used with other languages or in other contexts \item [(iii)] It can be used to illustrate ideas in probability and statistics \end{itemize} In MAS1608 last year we focused mainly on point (ii) above; in this module, we focus mainly on points (i) and (iii). \\ \\ In the first 3 weeks of the module we will go over the basics of programming in R -- this will also hep you with your compuing in other statistics modules. For the rest of semester 1, we will present a mixture of background theory in probability and statistics, together with elements of the R programming language. Thus, from around the fourth week of the course, activities will contain both some statistical theory and some R programming concepts. It is vital that, if possible, you attend your lectures and practical sessions. This is a great opportunity to ask us questions face-to-face, seek help from other students, and of course prepare for the two tests. It is even more important that you follow instructions on Canvas and complete the other activities by the dates suggested -- if you are told that some activities should be completed before the next computer practical, then make sure you do exactly that -- the work set in the computer practical class will assume you have completed all work to-date! % 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}} \newpage \lhead{MAS2906} \rhead{Chapter 2: R Refresher} \section{R Refresher} \subsection{Preliminaries} R is installed on all University machines. In this module we will be using RStudio, which is an R \textit{Integrated Development Environment}, or IDE. 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). You might remember from MAS1608 that both R and RStudio are free, so you can install them on your own computer. Versions of R and RStudio are available for Windows, Apple Mac and Linux. See the \href{https://ncl.instructure.com/courses/32642/pages/r-computing-resources?module_item_id=1231374}{Computing Resources} section of the Maths, Stats \& Physics Canvas Community for more information. \\ \\ We will assume a limited working knowledge of R and RStudio. The first three weeks are designed to recap the basics from MAS1608. Specifically, we will review (i) data types and vectors in R; (ii) dataframe manipulation; (iii) data summaries; (iv) functions and control statements. Everything we do for the remainder of semester 1 will then assume you have a good grasp of these basics, so please work hard to make sure you're up-to-speed with everything as we progress through this Section of the notes! And of course, get in touch if you're not sure about anything -- either by email or via the Canvas discussion board. You might recall the \texttt{movies} dataset from MAS1608. 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. It claims to be the biggest movie database on the web and is run by \textit{Amazon}. More information can be found online: \begin{center} \url{http://imdb.com/help/show_leaf?about} \end{center} This includes information about the data collection process. 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! See Table \ref{table:movies} for an extract of the movies dataset. \input movietable We have written an add-on R package for the module which will allow easy access to many interesting datasets (including the movies dataset). Installing the package is straightforward, and will make loading datasets used in the notes really easy: % 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} \vspace{-0.5cm} To load the package, type: \begin{lstlisting} library(ncldata) \end{lstlisting} \vspace{-0.5cm} Then, to specifically load the movies dataset, type: \begin{lstlisting} data(movies) \end{lstlisting} \vspace{-0.5cm} Alternatively, some datasets can be downloaded directly from the course Canvas page, and some from web sources; examples of how to do this will be given later in the notes. \subsection{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} Another important data type is \texttt{NA}, short for `Not Available' and used for missing values. A summary of the data types in R is given in Table \ref{table:datatypes} below. Note the way R presents doubles in standard form (Example 4); notice also that \underline{anything} declared inside inverted commas will be assumed by R to be of Character type. {\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}} \subsection{Vectors} Vectors are the most basic of all data structures, but are used in almost all R code. An R vector contains $n$ values of the same type, where $n$ can be zero. 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} \vspace{-0.5cm} Some notes on the above code listing: \begin{itemize} \item In lines 1-2 we see a vector of length $n=6$, containing the values $\{0, 1, \ldots, 5\}$. Note that this vector has not been assigned for future use \item In lines 3-4, we see the same vector as in lines 1-2, but this vector has now been assigned to the object ``first\_vec" \item The outer round brackets used in line 3 are not necessary - they are used in these notes and slides to ensure the resulting vector is printed to the notes/screen, for your convenience so you can see what the code has produced \item ``first\_vec'' is a vector of doubles \item ``second\_vec'' is a vector of $n = 3$ characters \item Both x and y are vectors too -- each of length $n=1$ (x is a vector containing $n=1$ double; y is a vector containing $n=1$ character) \end{itemize} 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} 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} And you can use the \texttt{seq} command to generate sequence vectors very quickly: \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 -2.0 \end{lstlisting} The following code listing demonstrates some other built-in R commands that can be very useful when applied to a vector -- \texttt{rev}, \texttt{sort}, \texttt{sum}, \texttt{unique} and \texttt{which}. Make sure you understand what each of these is doing -- especially \texttt{unique} and \texttt{which}! \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 for 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 which(x3<=1) [1] 5 6 7 8 9 10 11 \end{lstlisting} \subsubsection{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} \subsubsection{Extracting elements from vectors}\label{section:vectorextraction} R has a number of useful methods that we can use to extract subsets of our data. For example, we could use square brackets to extract particular elements: \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}{We only want values of z that are less than 3 \textbf{or} at least 7}*) [1] 1 2 7 8 9 10 \end{lstlisting} We will come back to extraction from two-dimensional R objects (e.g. dataframes) shortly. %\\ %\\ %\textbf{Always use round brackets to implement a function, and square brackets for extraction.} Another useful way of extracting elements from vectors is to use logical vectors. 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} 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 only extract the elements in vec1 where logic1 is true}*) [1] 3 5 \end{lstlisting} We can also sum a logical vector. For example: \begin{lstlisting} sum(logic1) [1] 2 sum(vec1[logic1]) [1] 8 \end{lstlisting} \newpage \begin{example} 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 \href{https://ncl.instructure.com/courses/32642/pages/r-computing-resources?module_item_id=1231374}{Computing Resources} section of the Maths, Stats \& Physics Canvas Community. Here you will find some basic resources aimed at students who have never used R before; scroll down, and you'll find information about how to install R on your home PC. 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("http://www.mas.ncl.ac.uk/~nlf8/teaching/mas2906/x.txt") y = as.logical(scan("http://www.mas.ncl.ac.uk/~nlf8/teaching/mas2906/y.txt", what="logical")) \end{lstlisting} \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} \begin{center} \begin{tabular}{l l} Your answers about x &Your answers about y \\ \hline & \\ (a) &(a) \\ & \\ (b) &(b)\\ & \\ (c) \textcolor{white}{What is the smallest value in the vector x?}&(c) \\ & \\ (d) &(d) \\ & \\ (e) &(e) \textcolor{white}{When is the 472nd occurrence of TRUE?}\\ & \\ (f) &(f) \\ & \\ \end{tabular} \end{center} \end{example} \newpage \subsection{Dataframes}\label{sec:dataframes} A dataframe is a special kind of object in R. We use dataframes for storing and managing datasets that have a rectangular structure. Typically the rows correspond to \textit{cases} and the columns \textit{variables}. 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. For example, we might have survey results on 10 individuals, who have each been asked 5 questions. This would give a dataframe with 10 rows and 5 columns (a so-called ``10 by 5 dataframe"; generally, we refer to $r \times c$ dataframes). The variable types need not all be the same -- the first column could be a vector of doubles recording the individuals' height in metres; the second column could be a vector of logicals taking the value \texttt{TRUE} if they like Marmite; the third column could be a vector of characters taking the name of their favourite book, etc. \subsubsection{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} \vspace{-0.25cm} So, case 2 in our survey is 1.75 metres tall, does \textit{not} like Marmite and their favourite book is ``To Kill a Mockingbird". We do not know the weights of cases 1 and 2 (missing values). Here we have 3 cases and 4 variables, so we have a $3 \times 4$ dataframe. If you had a large dataframe and wanted to quickly know its size, you could use the \texttt{dim} command (short for ``dimension"). For example: \begin{lstlisting} dim(small.survey) [1] 3 4 \end{lstlisting} \vspace{-0.25cm} We can extract elements from a dataframe using square brackets -- just as we demonstrated in Section \ref{section:vectorextraction} for extraction from a vector. However, if the columns in our dataframe have labels, we can also use the \$ command to access columns directly. For example: \begin{lstlisting} small.survey[2, 4] #%* \textcolor{ForestGreen}{Will extract the element from the 2nd row, 4th column}*) [1] To Kill a Mockingbird small.survey[ , 4] #%* \textcolor{ForestGreen}{Will extract the elements from \textit{all} rows in the 4th column}*) [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$book[2] #%* \textcolor{ForestGreen}{Does exactly the same as line 1!}*) [1] To Kill a Mockingbird small.survey[3, ] #%* \textcolor{ForestGreen}{Will extract the elements from \textit{all} columns in the 3rd row}*) height_m weight_kg marmite book 3 1.58 70.3 FALSE Life of Pi \end{lstlisting} Though used less, the \texttt{subset} command can be applied to dataframes 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} 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} \subsubsection{The movies dataframe}\label{sec:moviesdataframe} The Movies dataset is available from within the \texttt{ncldata} package. 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} \vspace{-0.7cm} \begin{itemize} \item Line 1 installs the \texttt{ncldata} package to your machine. You will only need to do this once, and in sessions in the computer cluster this has been done for you already. You will make use of this package for other datasets in Statistics modules, including MAS2903 \item Line 3 reads the package into your current working RStudio space -- each time you need to access data from within the package you should use this line of code \item Line 4 then pulls the \texttt{movies} dataset from the \texttt{ncldata} package, ready for you to use \end{itemize} You can now explore the \texttt{movies} dataset: \begin{lstlisting} dim(movies) #%* \textcolor{ForestGreen}{Looks like there are 4988 movies and 23 variables}*) [1] 4988 23 head(movies) #%* \textcolor{ForestGreen}{Dataframe too large to view; let's look at the first six rows}*) Title Year Rating Votes Length Budget BudgetCurrency 1 10,000 BC 2008 5.1 115607 109 1.05e+08 USD 2 101 Dalmatians (I) 1961 7.2 126498 79 NA 3 101 Dalmatians (II) 1996 5.7 84527 103 7.50e+07 USD 4 102 Dalmatians 2000 4.9 29328 100 8.50e+07 USD 5 10 Cloverfield Lane 2016 7.2 226726 104 1.50e+07 USD 6 10 Things I Hate About You 1999 7.2 249949 97 3.00e+07 USD mpaa Action Adventure Animation Biography Comedy Crime Documentary Drama 1 PG-13 TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE 2 FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE 3 FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE 4 FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE 5 PG-13 FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE 6 PG-13 FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE Fantasy Horror Mystery Romance SciFi Short Thriller 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE 3 FALSE FALSE FALSE FALSE FALSE FALSE FALSE 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE 5 FALSE TRUE TRUE FALSE FALSE FALSE FALSE 6 FALSE FALSE FALSE TRUE FALSE FALSE FALSE \end{lstlisting} From the above R code, we can see that the variables in the \texttt{movies} dataframe are: \begin{itemize} \item \textbf{Title} of the movie \item The \textbf{Year} the movie was released \item The movie's average \textbf{Rating} (out of 10), based on the number of \textbf{Votes} submitted \item The \textbf{Length} of the movie (in minutes) \item The movie \textbf{Budget} (this is in USD, as shown by \textbf{BudgetCurrency}, unless otherwise stated) \item The movie's rating as suggested by the Motion Picture Association of America (\textbf{mpaa}) \item Various logical variables indicating the movie's genre (e.g. \textbf{Animation}, \textbf{Horror} etc.) \end{itemize} For example, the movie ``10 Cloverfield Lane" was released in 2016; it has an average rating of 7.2/10 from 226, 726 votes submitted to the Internet Movie Database (IMDB); the movie length is one hour 44 minutes; the budget was around US\$ 15M; it was rated PG-13 by the Motion Picture Association of America; and it has been classified as a drama/horror/mystery movie. \subsubsection*{Examples of querying the movies dataframe} \begin{itemize} \item [1.] Note that there are some missing values -- for example, the budget for the original version of 101 Dalmatians, released in 1961, is unknown. What if we wanted to know how many movies had known budgets? \item [2.] What if we wanted to create a new dataset with only those movies made since the year 1998? \item [3.] Or (harder) movies that cost more than US\$ 100,000 but are not R rated? \end{itemize} For the first query: \begin{lstlisting} x1 = is.na(movies$Budget) #%* \textcolor{ForestGreen}{Will return TRUE if the budget is NA}*) head(x1) #%* \textcolor{ForestGreen}{Compare to dataset above -- looks good}*) [1] FALSE TRUE FALSE FALSE FALSE FALSE x2 = !is.na(movies$Budget) #%* \textcolor{ForestGreen}{Flip things round: NA = FALSE}*) head(x2) [1] TRUE FALSE TRUE TRUE TRUE TRUE budget_known = movies$Budget[x2] #%* \textcolor{ForestGreen}{Apply logical vector to the movie budgets}*) head(budget_known) [1] 1.05e+08 7.50e+07 8.50e+07 1.50e+07 3.00e+07 6.00e+06 length(budget_known) #%* \textcolor{ForestGreen}{This will tell us how many movies have known budgets}*) [1] 3990 length(movies$Budget[!is.na(movies$Budget)]) #%* \textcolor{ForestGreen}{A bit more concise!}*) [1] 3990 \end{lstlisting} \newpage Second query: \begin{lstlisting} movies2 = subset(movies, Year>=1998) dim(movies2) [1] 3478 23 \end{lstlisting} And the last query: \begin{lstlisting} movies3 = movies[movies$Budget>100000 & movies$mpaa != "R", ] dim(movies3) [1] 3174 23 \end{lstlisting} \newpage \begin{example}\label{example:movies} Follow the steps at the start of Section \ref{sec:moviesdataframe} to load the \texttt{ncldata} package and to read in the \texttt{movies} dataframe. Then answer the following questions. \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}"? \textbf{Answers:}\underline{\textcolor{white}{Answers to be inserted here!!}} \item [(b)] How many films are classified as ``Action"? How many are classified as ``Romance"? How many ``Short" films are there in the database? \textbf{Answers:}\underline{\textcolor{white}{Answers to be inserted here!!}} \item [(c)] How many films were made between the Years 1970 and 2000 (inclusive)? \textbf{Answer:}\underline{\textcolor{white}{Answers}} \item [(d)] How many movies are longer than two hours in duration? \textbf{Answer:}\underline{\textcolor{white}{Answers inserted!}} \item [(e)] How many movies made since the year 2005 (not including 2005) are longer than three hours in duration? \textbf{Answer:}\underline{\textcolor{white}{Answers to be inserted here!!}} \end{itemize} \end{example} \vspace{1cm} \begin{example}\label{example:music} Load and explore the \texttt{music} dataset from the \texttt{ncldata} package by typing the following lines of code (Note: 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 \ref{example:movies} and Section \ref{sec:moviesdataframe}): \begin{lstlisting} data(music) #%* \textcolor{ForestGreen}{What's this dataset all about? Let's find out}*) ?music \end{lstlisting} \vspace{-0.7cm} \begin{itemize} \item [(a)] How many tracks do we have data for, and on how many variables? \textbf{Answers:}\underline{\textcolor{white}{Answers to go h}} \item [(b)] How many tracks do we have from \textit{The Beatles} in this dataframe? \textbf{Answer:}\underline{\textcolor{white}{Answers to go he} }\\ What about \textit{Lady GaGa}? \textbf{Answer:}\underline{\textcolor{white}{Answers to be inserted here!!}} \item [(c)] How many tracks have a tempo of at least 130 beats per minute? \textbf{Answer:}\underline{\textcolor{white}{Answers to go here}} \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)? \textbf{Answer:}\underline{\textcolor{white}{Answers to be inserted here!!}} \item [(e)] What is the most popular key for tracks in this database? \textbf{Answer:}\underline{\textcolor{white}{Answers to be inserted her}} \end{itemize} \end{example} \newpage \subsection{Graphical and numerical summaries} R has very powerful built-in functions for summarising datasets both graphically and numerically. Most of these are fairly intuitive, so in this section we provide only a brief review of material from stage 1. For example, we might want to explore the distribution of budgets in the \texttt{movies} dataframe. As we know, there are some missing values in the dataset, but most of the built-in functions for data analysis in R have an argument we can invoke to ignore these: \begin{lstlisting} #%* \textcolor{ForestGreen}{To find the mean of the movie budgets but ignoring the missing values}*) mean(movies$Budget, na.rm = TRUE) [1] 41245023 #%* \textcolor{ForestGreen}{And the standard deviation}*) sd(movies$Budget, na.rm = TRUE) [1] 68612118 #%* \textcolor{ForestGreen}{How about a range of numerical summaries with a single function?}*) summary(movies$Budget, na.rm = TRUE) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.000e+00 9.000e+06 2.500e+07 4.125e+07 5.000e+07 1.900e+09 998 \end{lstlisting} From the listing above we can make a few observations: \begin{itemize} \item The mean budget is larger than the median: $\sim$US\$ 41M, compared to US\$ 25M. This indicates a positively skewed dataset, as the mean is being `pulled up' by some very large budgets (of around USD\$1.9Bn?! Which movie is this?! Must be amazing...) \item The standard deviation seems extremely large; again, this could be due to a small number of extremely large budgets... \item ... or perhaps some unrealistically \textit{small} budgets (of US\$ 0?!) \end{itemize} \begin{lstlisting} #%* \textcolor{ForestGreen}{Let's find out the position of the most expensive movie}*) which(movies$Budget == max(movies$Budget, na.rm=TRUE)) [1] 1211 #%* \textcolor{ForestGreen}{What's the name of this movie?}*) movies$Title[1211] [1] "Enthiran" #%* \textcolor{ForestGreen}{Never heard of it!}*) #%* \textcolor{ForestGreen}{Which movies were made for free?!}*) movies$Title[which(movies$Budget==0)] [1] "Homefront" "Just Friends" "Knock Knock" "Now You See Me" [5] "Safe House" "Snitch" \end{lstlisting} What about some graphical summaries? By far the most commonly-used graphics are histograms and boxplots. Of course, there are so many more graphical options available in R; you should consult your MAS1608 notes for a more comprehensive coverage. \begin{lstlisting} #%* \textcolor{ForestGreen}{Partitions the plotting space: 1 row, 3 columns}*) par(mfrow = c(1,3)) #%* \textcolor{ForestGreen}{Now we'll try three different methods for producing histograms, based on class width/number of bars}*) #%* \textcolor{ForestGreen}{1. Sturges' rule; the default, tries to find an optimal number of bars}*) hist(movies$Budget, col="wheat",breaks="sturges", xlab="Budget (USD)", main="Sturges' rule") #%* \textcolor{ForestGreen}{2. Freedman-Diaconis's rule; tries to optimise the class width}*) hist(movies$Budget, col="gold",breaks="fd", xlab="Budget (USD)", main="FD rule") #%* \textcolor{ForestGreen}{3. Scott's rule; also tries to optimise the class width}*) hist(movies$Budget ,breaks="scott", xlab="Budget (USD)", main="Scott's rule") \end{lstlisting} \begin{lstlisting} par(mfrow=c(1,2)) #%* \textcolor{ForestGreen}{Let's try a boxplot now}*) boxplot(movies$Budget, col="gold", main="Boxplot of movie budgets", horizontal=TRUE, xlab="Budget (USD)") #%* \textcolor{ForestGreen}{Outliers are spoiling things a bit. Let's remove them}*) boxplot(movies$Budget, col="gold", main="Boxplot of movie budgets", horizontal=TRUE, xlab="Budget (USD)", outline=FALSE) \end{lstlisting} The histograms shown in Figure \ref{fig:hist} reveal just how positively skewed the movie budget data are. They also show how care should be taken with the method chosen to produce the histogram. Of course, this is a subjective choice, and some experimentation might be necessary -- \textit{especially} with such a heavily skewed dataset. The boxplots also indicate just how skewed our dataset is; to get a handle on the distribution of most of of our data, it might be preferable to produce the boxplot with the outliers removed. Table \ref{table:summaries} gives a summary of some R commands for commonly-used tables and graphs. \begin{table}[!b] \centering \begin{tabular}{l l l} \hline Command &Comment &Example \\ \hline \texttt{table} &Contingency table &\texttt{table(x)} \\ \texttt{barplot} &Generate a bar chart &\texttt{barplot(table(x))} \\ \texttt{hist} &Histogram &\texttt{hist(x)} \\ \texttt{plot} &Scatterplot &\texttt{plot(x, y)} \\ \texttt{points} &Add points to a plot &\texttt{points(x, y)} \\ \texttt{lines} &Add lines to a plot &\texttt{lines(x, y)} \\ \texttt{boxplot} &Box and whisker plot &\texttt{boxplot} \\ \hline \end{tabular} \caption{Summary of R commands for commonly-used tables/graphs/graphical features} \label{table:summaries} \end{table} \begin{center} \begin{figure}[!p] \centering \includegraphics[]{hist_2905.pdf} \includegraphics[]{boxplot_2905.pdf} \caption{Histograms for the movie budgets, constructed using three different methods (top row); boxplots (bottom row), with outliers removed (right).} \label{fig:hist} \end{figure} \end{center} \newpage \begin{example}\label{example:music2} Load the \texttt{music} dataset into RStudio, just as you did in Example \ref{example:music}. \begin{itemize} \item [(a)] Find the mean and standard deviation for the duration of tracks in this dataset. \item [(b)] Produce histograms for track duration and danceability; put these histograms on the same panel of plots, use different colours for the histograms (you can choose whichever colours you like!) and make sure you provide suitable titles and labels on the axes. \item [(c)] You are asked to investigate whether or not the duration of tracks has changed over time. \begin{itemize} \item [(i)] Use the \texttt{plot} command in R to produce a scatterplot of the year each track was released against track duration. \item [(ii)] Try out the following code, and comment on what you see in the plots in part (c). \begin{lstlisting} boxplot(music$duration~music$decade) \end{lstlisting} \vspace{-0.7cm} \end{itemize} \item [(d)] On average, how do energy scores for ``Hip Hop" tracks compare to ``Electronic" tracks? \end{itemize} \end{example} \vspace{1cm} \begin{example} The code in the following listing will download COVID-19 data from the website of the \textit{European Centre for Disease Prevention and Control}. Note: this dataset contains information for 2019-2020 only. \begin{lstlisting} covid = read.csv("https://opendata.ecdc.europa.eu/covid19/casedistribution/ csv", na.strings = "", fileEncoding = "UTF-8-BOM") \end{lstlisting} \vspace{-0.7cm} \begin{itemize} \item [(a)] Use the \texttt{head} command to take a look at the first 6 rows of the dataframe, to get a feel for the information contained in the dataset. Use the \texttt{dim} command to find out the size of the dataframe. \item [(b)] We can use the \texttt{attach} command to directly access the columns in the dataframe by their headings. This might make it a bit easier to isolate data for certain countries. For example: \begin{lstlisting} #%* \textcolor{ForestGreen}{Enable access to columns via their headings}*) attach(covid) #%* \textcolor{ForestGreen}{Find the total number of cases, and deaths}*) sum(cases); sum(deaths) [1] 71503614 [1] 1612833 #%* \textcolor{ForestGreen}{What about Great Britain?}*) sum(cases[countryterritoryCode=="GBR"], na.rm=TRUE) [1] 1849403 sum(deaths[countryterritoryCode=="GBR"], na.rm=TRUE) [1] 64170 \end{lstlisting} \newpage \begin{lstlisting} #%* \textcolor{ForestGreen}{What about a plot of cases/deaths through time in Great Britain? See Figure \ref{fig:covid}}*) par(mfrow=c(1,2)) plot(rev(cases[countryterritoryCode=="GBR"]), type="l", xlab="Day", ylab="Confirmed cases", main="GBR COVID-19 cases") plot(rev(deaths[countryterritoryCode=="GBR"]), type="l", xlab="Day", ylab="Confirmed deaths", main="GBR COVID-19 deaths") \end{lstlisting} \vspace{-0.7cm} \begin{itemize} \item [(i)] Find the total number of confirmed cases, and deaths, due to COVID-19, for three other countries of your choosing. \item [(ii)] Why is it necessary to use the \texttt{rev} command in lines 3 and 4 in the above listing? \item [(iii)] For the three other countries you chose in part (i), produce plots showing confirmed cases and deaths through time. \end{itemize} \end{itemize} \begin{center} \begin{figure}[!h] \centering \includegraphics[]{covid19v2.pdf} \caption{Daily cases and deaths due to COVID-19 in Great Britain.} \label{fig:covid} \end{figure} \end{center} \end{example} \subsection{Functions, for loops and if statements} The final part of this R Refresher will focus on what was probably the biggest part of the stage 1 R course: functions, for loops and if statements. In this section of the notes, we will focus on a single example that will pull all of these aspects of R programming together: How to fill a sticker album! The 2018 \textit{Panini} sticker album for the Fifa world cup football finals had a whopping 682 stickers to collect (including 40 awesome glitter stickers!). Each packet of stickers cost 80p and contained 5 stickers. The obvious question to ask is: \begin{center} \textit{``How much does it cost to fill a Panini World Cup `18 sticker album?"} \end{center} Now look at the code in the listing below. \begin{lstlisting} #%* \textcolor{ForestGreen}{Create a vector for the sticker sequence - suppose for now that we only have 50 stickers}*) (stickers = seq(1, 50)) #%* \textcolor{ForestGreen}{Can actually be written more concisely as \texttt{stickers = 1:50}}*) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 [25] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 [49] 49 50 #%* \textcolor{ForestGreen}{Now let's buy a packet!}*) (packet = sample(1:50, 5, replace=TRUE)) #%* \textcolor{ForestGreen}{Random sample of size 5 from $\{1, 2, \ldots, 50\}$}*) [1] 10 23 20 48 11 #%* \textcolor{ForestGreen}{We got stickers 10, 11, 20, 23 and 48!}*) #%* \textcolor{ForestGreen}{As a way of keeping track of the stickers we have collected, we could set the corresponding sticker}*) #%* \textcolor{ForestGreen}{sequence entry to zero}*) stickers[packet[1]] = 0 #%* \textcolor{ForestGreen}{Will set \texttt{stickers[10]} to 0}*) stickers[packet[2]] = 0 #%* \textcolor{ForestGreen}{Will set \texttt{stickers[23]} to 0}*) stickers[packet[3]] = 0 #%* \textcolor{ForestGreen}{Will set \texttt{stickers[20]} to 0}*) stickers[packet[4]] = 0 #%* \textcolor{ForestGreen}{Will set \texttt{stickers[48]} to 0}*) stickers[packet[5]] = 0 #%* \textcolor{ForestGreen}{Will set \texttt{stickers[11]} to 0}*) #%* \textcolor{ForestGreen}{That's a bit repetitive! We could put this inside a loop}*) for(j in 1:5){ stickers[packet[j]] = 0} stickers [1] 1 2 3 4 5 6 7 8 9 0 0 12 13 14 15 16 17 18 19 0 21 22 0 24 [25] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 0 [49] 49 50 \end{lstlisting} After running the code in the above listing a few times, you should start to get a feel for what's going on. The aim is to keep buying packets of stickers and setting the corresponding entry in the \texttt{sticker} sequence to zero as a way of keeping track of the stickers we have collected. The code above adjusts this sequence after buying just a single packet of stickers. We need to do this repeatedly, an unknown (and probably very large!) number of times. So, another \texttt{for} loop is required. \begin{lstlisting} #%* \textcolor{ForestGreen}{Re-set the sticker sequence}*) (stickers = seq(1, 50)) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 [25] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 [49] 49 50 #%* \textcolor{ForestGreen}{Now let's start a big for loop, in which we repeatedly buy packets of stickers}*) for(i in 1:10000){ #%* \textcolor{ForestGreen}{10,000 is a big number - surely this will be enough packets?}*) packet = sample(1:50, 5, replace=TRUE) for(j in 1:5){ stickers[packet[j]] = 0}#%* \textcolor{ForestGreen}{Close the j loop}*) }#%* \textcolor{ForestGreen}{Close the i loop}*) stickers [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [37] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 \end{lstlisting} Looking at the listing above, it looks like we've filled the sticker album after buying 10,000 packets of stickers -- notice that all entries in the sticker sequence vector are zero. But, we don't know \underline{exactly} when this happened. Did it happen after we bought the 10,000th packet? Or after we bough the 200th packet? Or the 1524th packet?! We need to `break' the code as soon as all entries in the sticker sequence are zero. \begin{lstlisting} (stickers = seq(1, 50)) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 [25] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 [49] 49 50 for(i in 1:10000){ packet = sample(1:50, 5, replace=TRUE) for(j in 1:5){ stickers[packet[j]] = 0} #%* \textcolor{ForestGreen}{Now inside the for loop use an if statement to break the code when all entries in the sticker sequence are zero}*) if(sum(stickers)==0) { break} } #%* \textcolor{ForestGreen}{When did the code break?}*) i [1] 37 \end{lstlisting} So this time it only took 37 packets of stickers to completely fill the album! Of course, in this toy example we only have 50 stickers to collect. For the \textit{Panini} football sticker album, there are 682! Really, what we'd like to do is run the code in the above listing repeatedly, to see what happens for different `people' trying to fill their sticker albums. It might be a good idea to put all of this into a \texttt{function} now. \begin{lstlisting} FillAlbum = function(N, k){ #%* \textcolor{ForestGreen}{N = something large, k = number of stickers to collect}*) stickers = seq(1, k) for(i in 1:N){ packet = sample(1:k, 5, replace=TRUE) for(j in 1:5){ stickers[packet[j]] = 0} if(sum(stickers)==0) { break} } return(i)} \end{lstlisting} Notice, from the code listing above, the general structure of a function in R: \begin{itemize} \item The word \texttt{function} (line 1), followed by round brackets and some ``arguments'' to be passed through the function (arguments that we must specify when we execute the function) \item An appropriate name for the function (\texttt{FillAlbum} here, see line 1) \item The use of braces surrounding the operations within the function (lines 1 and 13; all other braces relate to the for loops/if statement) \item A return statement (line 13), declaring what we would like the function to spit out \end{itemize} We can now execute our \texttt{FillAlbum} function to see how many packets of stickers different ``people'' might need to buy in order to fill their albums. \begin{lstlisting} #%* \textcolor{ForestGreen}{Let's try for an album containing just 50 stickers again}*) FillAlbum(10000, 50) [1] 52 FillAlbum(10000, 50) [1] 32 FillAlbum(10000, 50) [1] 40 #%* \textcolor{ForestGreen}{Now let's try for the full album of 682 stickers!}*) FillAlbum(10000, 682) [1] 1290 FillAlbum(10000, 682) [1] 1611 FillAlbum(10000, 682) [1] 1252 \end{lstlisting} Referring to the code in the listing above, we have simulated three people trying to fill the full \textit{Panini} sticker album. Respectively, they needed to buy 1,290, 1,611 and 1,252 packets of stickers to do this. We could now execute this function many times, within another \texttt{for} loop. \begin{lstlisting} #%* \textcolor{ForestGreen}{Let's simulate 1000 people trying to fill their albums}*) packets_needed = vector("numeric", 1000) for(k in 1:1000){ packets_needed[k] = FillAlbum(10000, 682)} #%* \textcolor{ForestGreen}{Each packet costs 80p}*) cost = packets_needed*0.8 mean(packets_needed); mean(cost) [1] 974.609 [1] 779.6872 sd(packets_needed); sd(cost) [1] 181.351 [1] 145.0808 par(mfrow=c(1,2)) hist(packets_needed, col="misty rose", xlab="Packets", main="Packets needed") hist(cost, col="chocolate", xlab="Pounds",main="Cost") \end{lstlisting} \newpage \begin{center} \begin{figure}[!h] \centering \includegraphics[]{sticker.pdf} \caption{Histograms showing the number of packets of stickers needed to fill the \textit{Panini} football sticker album (left), and the associated cost (right).} \label{fig:sticker} \end{figure} \end{center} This analysis is really insightful! We can see that, on average, it takes around 975 packets of stickers to fill the album (standard deviation 181), with an associated cost of around \textrm{\pounds}780 (standard deviation \textrm{\pounds}145). The histograms in Figure \ref{fig:sticker} reveal the full distribution for the number of packets needed and the associated cost. The larger the number of people simulated (we used 1,000 here), the more accurate this distribution will be! \newpage \lhead{MAS2906} \rhead{Chapter 3: Simulating Random Variables} \section{Simulating Random Variables} \subsection{Discrete distributions: Binomial, geometric and Poisson} \noindent\textbf{The binomial distribution:} If $X\sim Bin(n,p)$ then $X$ has PMF\ (probability mass function) given by: \scalebox{1.1}{\includegraphics[trim={1cm 25cm 2cm 2cm},clip,page=9]{white}} \vspace{1cm} There is no closed formula for the CDF (cumulative distribution function; i.e. $\pr{X\leq x}$). The R commands \texttt{dbinom}, \texttt{pbinom}, and \texttt{rbinom} are used respectively to compute the PMF, compute the CDF, and simulate from a binomial distribution. Here are examples: \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} Here's an example of how to create a bar plot of a sample from a binomial distribution: %<>= %pdf(file="fig1-1.pdf",width=3.5,height=3) %@ %\begin{Schunk} %\begin{Sinput} \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} %\end{Sinput} %\end{Schunk} %<>= %dev.off() %@ \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{fig1-1} \end{center} \newpage \noindent\textbf{The geometric distribution:} If $Y\sim Geom(p)$ then $Y$ has PMF and CDF given by: \scalebox{1.2}{\includegraphics[trim={1cm 17.5cm 2cm 6.8cm},clip,page=9]{white}} Note that $Y$ takes values in $1,2,3,\ldots$. R uses a slightly different definition of the geometric distribution. We use the definition that $Y$ is the number of Bernoulli trials with parameter $p$ up to and including the first trial which is a success. The definition in R is the number of trials up to, but \textbf{not} including the first success, so in R geometric random variables take values $0,1,2,\ldots$. Like the binomial distribution, the R commands \texttt{dgeom}, \texttt{pgeom}, and \texttt{rgeom} are used respectively to compute the PDF, compute the CDF, and simulate from the geometric distribution. However, they need to be adjusted to match our 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} \noindent\textbf{The Poisson distribution:} If $Z\sim Po(\lambda)$ then $Z$ has PMF given by: \scalebox{0.9}{\includegraphics[trim={1cm 15cm 2cm 12.5cm},clip,page=9]{white}} There is no closed formula for the CDF. Examples of R commands: \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} \noindent\textbf{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} \subsection{Continuous random variables} R has functions for the uniform, exponential and normal distributions, where we now refer to the probability \textit{density} function, rather than probability \textit{mass} function: \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} Examples of R commands: \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} \vspace{-0.75cm} 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} \vspace{-0.75cm} Quantile functions such as \texttt{qexp} and \texttt{qnorm} solve equations of the form $F_X(\alpha) = p$, for $\alpha$, given a probability $p$. For example: \begin{lstlisting} qnorm(0.9750021) [1] 1.96 \end{lstlisting} \vspace{-0.4cm} \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. \textbf{Answer:} \underline{\textcolor{white}{Answer here...}} \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? \textbf{Answer:} \underline{\textcolor{white}{Answer here...Answer here... Answer here...Answer}} \end{itemize} \end{example} \vspace{0.15cm} \begin{example} Use R to generate a sample of size 10,000 from an exponential distribution with rate 2.5, and store this sample in the vector x. \begin{itemize} \item [(a)] Use your vector x to estimate $\pr{X \leq 0.5}$, where $X \sim \textit{Exp}(2.5)$. \\ \textbf{Answer:} \underline{\textcolor{white}{Answer to go here}} \item [(b)] How does your estimate in part (a) compare to the true probability? Use R here -- do not perform any calculations by hand. \textbf{Answer:} \underline{\textcolor{white}{Answer here......Answer here please..... Answer...}} \end{itemize} \end{example} \begin{example} Write down the R code, and corresponding result, to answer each of the questions below. The first one has been done for you. \begin{center} \begin{tabular}{|c l |l |l|} \hline Model &Question &R code &Answer \\ \hline & & & \\ &Find $\pr{X=6}$. &\texttt{dpois(6,3)} &0.0504 \\ & & & \\ \cline{2-4} & & & \\ $X\sim Po(3)$ &Find $\pr{X<5}$. & & \textcolor{white}{Answers to be inserted} \\ & & & \\ \cline{2-4} & & & \\ &Generate 3 observations.& & \\ & & & \\ \hline & & & \\ &Find $\pr{Y\leq 4}$. & & \\ & & & \\ \cline{2-4} & & & \\ $Y \sim Bin(10, 0.2)$ &Find $\pr{Y\geq 4}$. & & \\ & & & \\ \cline{2-4} & & & \\ &Find $y$ such that $\pr{Y\leq y} = 0.677$. & & \\ & & & \\ \hline & & & \\ &Find $\pr{Z\leq 0.2}$. & & \\ & & & \\ \cline{2-4} & & & \\ $Z \sim N(0,1)$ &Find $\pr{Z > -2.5}$. & & \\ & & & \\ \cline{2-4} & & & \\ &Generate 3 observations. & & \\ & & & \\ \hline \end{tabular} \end{center} \end{example} \subsection{Some more advanced probability models} The main reason to use computers for Mathematics and Statistics is to solve problems which cannot be done by hand or which require very lengthy calculations. This section contains a simple motivating example to illustrate this. The number of arrivals per day at a university IT help-desk can be modelled using a Poisson distribution. However, it is thought that the mean of the Poisson distribution might vary from day-to-day. For example, the desk is often busier nearer the end of term due to more students having deadlines. The following model is proposed. Suppose the number of arrivals $X|\Lambda\sim Po(\Lambda)$ where the rate $\Lambda$ is itself a random variable, with $\Lambda\sim Exp(c)$ for a constant $c=0.05$. What can we say about the distribution of $X$? The random variable $X$ does not have a distribution which we can immediately recognize. In fact, much can be done analytically, but we will use simulation to answer the following questions: \begin{enumerate} \item What are the expectation and variance of $X$? \item What is $\pr{X>30}$? \end{enumerate} Here is code to simulate a number $n$ of samples from the \hi{distribution} of $X$: \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} \vspace{-0.5cm} Here's how to create a bar plot of the sample and compute the \hi{sample} mean and variance: %<>= %pdf(file="fig1-2.pdf",width=3.5,height=3) %@ \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} %<>= %dev.off() %@ \vspace{-1.5cm} \begin{center} \includegraphics[width=8cm]{fig1-2} \end{center} \vspace{-1cm} \begin{lstlisting} mean(x); var(x) [1] 19.618 [1] 382.1764 \end{lstlisting} \vspace{-0.5cm} The mean and variance of the distribution of $X$ will be close to the sample mean and variance that we see in the output above. To calculate $\pr{X> 30}$ we calculate the proportion of times this occurs in the sample: \begin{lstlisting} sum(x>30)/500 [1] 0.202 \end{lstlisting} The bigger the value we use for $n$, the closer the sample distribution will be to the true distribution for $X$; and the better our estimates of the mean, variance and probabilities will be. \newpage \subsection{Mixtures of normals} Suppose $\mu_1,\mu_2,\sigma_1>0,\sigma_2>0$ are fixed constants, and $w_1,w_2$ are positive constants with $w_1+w_2=1$. 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: \scalebox{1.2}{\includegraphics[trim={1cm 2cm 2cm 22cm},clip,page=9, height=3cm]{white}} 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*} The right hand side is actually another normal distribution (as seen in MAS1608), but $X$ \hi{does not have} the PDF of a normal random variable. For example: suppose $\mu_1=3, \sigma_1=1$ and $\mu_2=6, \sigma_2=2$ with $w_1=w_2=1/2$. Here is a plot of $f(x)$: \vspace{-1cm} \begin{center} \includegraphics[width=10cm]{fig1-3} \end{center} \newpage Consider how to sample $X$ from this 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$. 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} Here is an R function to perform this sampling $n$ times: \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} \hi{Write down any notes about this function.} Normal mixture distributions can contain more than two components -- the generalization is very straightforward. They offer a very flexible way to model data with different shaped distributions. \vspace{1cm} \begin{example} Investigations into infants suffering from \textit{severe idiopathic respiratory distress syndrome} have shown that whether the infant survives may be related to their weight at birth. Suppose that you are interested in developing a distribution for the birth weight $X$ of such infants. You have a normal $N(2.3, 0.52^2)$ distribution for the birth weight (in kg) of infants who survive and a normal $N(1.7, 0.66^2)$ distribution for infants who die. \begin{itemize} \item [(a)] If you believe that the proportion of infants that survive is 0.6, write down a suitable Normal mixture distribution for $X$. \item [(b)] Produce some R code to sample from the normal mixture distribution you propose in part (a). \item [(c)] Use your code in part (b) to simulate a sample of size 1000 from your normal mixture distribution, and produce a histogram of your sample. \end{itemize} \end{example} \newpage \lhead{MAS2906} \rhead{Chapter 4: The Inverse CDF Method} \section{The Inverse CDF Method} Let $U\sim U(0,1)$. Suppose $F(x)$ is a well-defined CDF which is invertible. Then the random variable $X=F^{-1}(U)$ has CDF $F(x)$. \textbf{Proof}: \scalebox{1}{\includegraphics[trim={0.5cm 16.9cm 1cm 1.2cm},clip,page=6]{white}} \vspace{1.75cm} \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} \scalebox{1.1}{\includegraphics[trim={0.5cm 5.2cm 1cm 15.2cm},clip,page=6]{white}} \vspace{5cm} \begin{example} Write some R code to sample $n$ observations from an \textit{Exp}(2) distribution. Do not use the built-in \texttt{rexp} function. Check your code by comparing a histogram of some simulated realisations against the PDF for an \textit{Exp}(2) distribution. Also check your code against results from the built-in \texttt{rexp} function. \end{example} \newpage \begin{example} The \textit{Pareto distribution}, named after the Italian civil engineer, economist, and sociologist Vilfredo Pareto, is a power-law probability distribution that is used in the description of social, scientific, geophysical, actuarial, and many other types of observable phenomena. It was originally used to describe the distribution of wealth in a society, fitting the trend that a large portion of wealth is held by a small fraction of the population. If $X \sim Pareto(\sigma)$,then it has cumulative distribution function \[ F_{X}(x; \sigma) = 1 - \left(\sigma/x\right)^{\log_{4}5}, \quad x \geq \sigma, \sigma>0. \] Write a function in R to produce realisations from this distribution, and test your function out for $\sigma = \{1, 2, 5\}$. \end{example} \vspace{7cm} \begin{example} The \textit{Weibull distribution} is often used in survival data analysis to model lifetimes (e.g. the lifetimes of patients receiving a particular drug). If $X \sim Weibull(\lambda, \kappa)$, then it has cumulative distribution function \[ F_{X}(x; \lambda, \kappa) = 1 - e^{-(x/\lambda)^{\kappa}}, \quad x\geq 0, \lambda, \kappa>0. \] Write a function in R to produce realisations from this distribution, and test your function out (a recent study suggested that, for a particular cancer drug, $\lambda= 2$ and $\kappa = 1.5$ for survival times in months after starting this drug). \end{example} \newpage \begin{example} Patients arriving at a doctor's surgery have appointment times $15$ minutes apart. The random variable $L$ measures how late a patient is for his/her appointment, with negative values denoting early arrivals. 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*} \begin{enumerate}[a.] \item[(a)] Sketch $f_L(\ell)$. \item[(b)] 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 [(c)] 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{example} \scalebox{0.9}{\includegraphics[trim={0.5cm 7cm 2cm 0.1cm},clip,page=7]{white}} \scalebox{1.2}{\includegraphics[trim={0.5cm 2cm 0.5cm 0.1cm},clip,page=8]{white}} \newpage \newpage \lhead{MAS2906} \rhead{Chapter 5: Transformation of Random Variables} \section{Transformation of Random Variables} \subsection{Transforming a normal random variable} The normal distribution is very widely used to model data. However, normal random variables take values on the entirety of $\R$ and they are symmetric around the mean. This is unrealistic in many situations -- for example, if we want to model people's heights or weights. One way to obtain a positive random variable is to define \begin{equation*} Y = e^{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$. What does the distribution of $Y$ look like? We can use sampling to find out, at least approximately: \begin{lstlisting} mu = 1.5 sigma = 0.3 n = 250 y = exp(rnorm(n, mean = mu, sd = sigma)) hist(y, xlab = 'y', ylab = 'frequency', main = 'Histogram of Y') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=9cm]{fig3-1new} \end{center} The random variable $Y$ can only adopt positive values, and the distribution is positively skewed (long tail on the right). This kind of distribution is useful to model quantities like height -- but what is the PDF of $Y$? A first guess is $f_Y(y)=f_X(\log(y))$ since $\log$ maps $Y$ back to $X$, but this is \hi{not} the correct answer. \newpage \subsection{Computing PDFs for transformed random variables} \begin{example} If $X\sim U(0,1)$, what is the PDF of $Y=X^3$? \end{example} \scalebox{1.2}{\includegraphics[trim={1cm 17cm 3cm 2cm},clip,page=1]{white}} To confirm this result, we can simulate $Y$. Does the histogram look right? \begin{lstlisting} y = (runif(500))^3 hist(y, xlab = 'Y', ylab = 'density', main="", prob=TRUE) curve((1/3)*(x^(-2/3)), from=0, to=1, col="purple", lty=2, add=TRUE) \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{fig3-2new} \end{center} \vspace{-0.5cm} \textbf{Notes:} \newpage \begin{example} If $X\sim N(0,1)$, what is the PDF of $Y=X^2$? \end{example} \scalebox{1.2}{\includegraphics[trim={1cm 12cm 3cm 14cm},clip,page=1]{white}} \scalebox{1.2}{\includegraphics[trim={1cm 5.5cm 3cm 17.5cm},clip,page=1]{white}} \begin{example} If $X\sim N(\mu,\sigma^2)$, what is the PDF of $Y=e^{X}$? This is the example used at the start of this section. We say $Y$ has the \hi{log-normal} distribution and write $Y\sim LN(\mu,\sigma^2)$. \end{example} \scalebox{1.2}{\includegraphics[trim={1cm 19cm 2cm 2cm},clip,page=2]{white}} \newpage \subsection{The transformation method} These examples are all special cases of the following result. If $X$ has PDF $f_X(x)$ and $h(x)$ is a monotonically increasing or monotonically decreasing function then $h$ is invertible and 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. This result allows direct calculation of $f_Y(y)$ and is called the \hi{univariate transformation method}. We will see this result again in MAS2903. \textbf{Proof:} \scalebox{0.8}{\includegraphics[trim={0.5cm 0.5cm 2cm 1.8cm},clip,page=3]{white}} \newpage \begin{example} When $X\sim U(0,1)$, use the transformation method to find the PDF of $Y=X^3$. \end{example} \scalebox{1.2}{\includegraphics[trim={1cm 23cm 2cm 1.5cm},clip,page=4]{white}} \vspace{1cm} \begin{example} When $X\sim U(-1,1)$ and $Y=e^{X}$, use the transformation method to find the PDF of $Y$. \end{example} \scalebox{1.2}{\includegraphics[trim={1cm 15cm 2cm 8.7cm},clip,page=4]{white}} \vspace{1cm} \begin{example} Use the transformation method to find the PDF of the log-normal distribution. \end{example} \scalebox{1.2}{\includegraphics[trim={1cm 8cm 2cm 16cm},clip,page=4]{white}} \vspace{1cm} \newpage \begin{example}\label{ex:cauchy} 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} \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]{fig3-3} \end{center} This distribution looks odd! Calculate the PDF:- \scalebox{1.2}{\includegraphics[trim={1cm 20cm 2cm 1.5cm},clip,page=5]{white}} \newpage \newpage \rhead{Chapter 6: Bivariate distributions} \section{Bivariate distributions} In this section we consider general simulation from a bivariate distribution, as well as simulation from the bivariate Normal distribution. \subsection{Sampling from bivariate distributions}\label{sec:bivart_sim} Suppose we have a pair of continuous random variables $(X,Y)$ with joint PDF $f_{X,Y}(x,y)$. Recall that \begin{equation*} f_{X,Y}(x,y) = f_X(x)f_{Y|X}(y\given X=x) \end{equation*} where $f_X(x)$ is the \hi{marginal} PDF of $X$ and $f_{Y|X}(y\given X=x)$ is the \hi{conditional} PDF of $Y$ given $X=x$. One way to simulate realizations of the pair $(X,Y)$ is: \begin{enumerate} \item Simulate a realization $x$ of $X$ from the marginal distribution of $X$. \item Simulate a realization $y$ of $Y$ from the conditional distribution of $Y$ given $X=x$. \end{enumerate} In the \hi{trivariate case} we can add a third step: simulate $Z$ given $Y=y$ and $X=x$. \vspace{0.5cm} \begin{example} Suppose $X\sim U(0.1,0.5)$ and \begin{equation*} f_{Y|X}(Y\given X=x) = \begin{cases} xe^{-xy}, &\text{when\ }y\geq 0,\\ 0&\text{otherwise}. \end{cases} \end{equation*} Write an R function to generate a $2\times n$ \hi{matrix} of samples from this distribution, and produce a scatter plot. \end{example} Note that, given that $X=x$, $Y$ has an exponential distribution with parameter $\lambda=X$: \begin{equation*} Y\given X=x \sim Exp(x). \end{equation*} \begin{lstlisting} example6.1 = function(n) { output = matrix(0, 2, n) x = runif(n, 0.1, 0.5) y = rexp(n, rate = x) output[1, ] = x output[2, ] = y return(output) } \end{lstlisting} \textbf{Make some notes about the R function.} \newpage To obtain a scatterplot: \begin{lstlisting} test = example6.1(200) plot(test[1,], test[2,], xlab = 'x', ylab = 'y') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{fig4-1} \end{center} \begin{example} (Trivariate case.) Suppose that $X\sim U(-1,1)$ and $Y\sim Exp(\lambda)$ independently of $X$ and that \begin{equation*} f_{Z|X,Y}(Z\given X=x, Y=y) = \frac{1}{\sqrt{2\pi}y}e^{-\frac{1}{2}\left( \frac{z-x}{y} \right)^2}. \end{equation*} Write an R function to sample realizations $(X,Y,Z)$ in the case $\lambda=5$ and produce a scatterplot of $(X,Z)$. What would the plot look like if $\lambda=1$? \end{example} Note that, given that $X=x$, $Y=y$ then $Z$ has a normal distribution: \begin{equation*} Z\given X=x, Y=y \sim N(x,y^2). \end{equation*} \begin{lstlisting} example6.2 = function(n) { output = matrix(0, 3, n) x = runif(n, -1, 1) y = rexp(n, rate = 5) z = rnorm(n, mean = x, sd = y) output[1, ] = x output[2, ] = y output[3, ] = z return(output) } \end{lstlisting} \textbf{Make some notes about the R function.} \newpage To obtain a scatterplot: \begin{lstlisting} output = example6.2(2000) plot(output[1,], output[3,], xlab = 'x', ylab = 'z') \end{lstlisting} \vspace{-1.4cm} \begin{center} \includegraphics[width=8cm]{fig4-2} \end{center} \vspace{-1cm} \subsection{The multivariate normal distribution} We met the bivariate normal distribution in MAS1608. The bivariate normal distribution is a generalization of the normal distribution to pairs of random variables, or equivalently, to a distribution on vectors in $\R^2$. The multivariate normal distribution is the analogous distribution on vectors in $\R^n$. \subsubsection{The bivariate normal distribution} \fbox{% \begin{minipage}{\textwidth} Suppose we have constants $\mu_x$, $\mu_y$, $\sigma_x$, $\sigma_y$ and $\rho$ such that $\sigma_x\geq 0$, $\sigma_y\geq 0$ and $-1\leq \rho\leq 1$. Define the $2\times 2$ matrix $\Sigma$ by \begin{equation*} \Sigma = \begin{pmatrix} \sigma_x^2 & \rho\sigma_x\sigma_y \\ \rho\sigma_x\sigma_y & \sigma_y^2 \end{pmatrix}. \end{equation*} Then define a joint probability density function by \begin{equation}\label{equ:MVN_pdf} f_{X,Y}(x,y) = \frac{1}{2\pi\sqrt{\det\Sigma}}\exp\left(-\frac{1}{2} Q(x,y) \right) \end{equation} where \begin{equation}\label{equ:quad_form} Q(x,y) = (\vec{x}-\vec{\mu})^T \Sigma^{-1} (\vec{x}-\vec{\mu}) \end{equation} and the column vectors $\vec{x},\vec{\mu}$ are defined by $\vec{x}=(x,y)^T$, $\vec{\mu}=(\mu_x,\mu_y)^T$. If random variables $(X,Y)$ have joint probability density given by $f_{X,Y}$ above, then we say that $(X,Y)$ have a \hi{bivariate normal distribution} and write \begin{equation*} (X,Y)^T\sim N_2(\vec{\mu},\Sigma). \end{equation*} \end{minipage}}\\ Some of the terms in the definition can be expanded: \begin{equation*} \sqrt{\det\Sigma} = \sigma_x\sigma_y\sqrt{1-\rho^2} \end{equation*} and \begin{align*} Q(x,y) &= \frac{1}{\sigma_x^2\sigma_y^2(1-\rho^2)} \begin{pmatrix}x-\mu_x & y-\mu_y\end{pmatrix} \begin{pmatrix} \sigma_y^2 & -\rho\sigma_x\sigma_y \\ -\rho\sigma_x\sigma_y & \sigma_x^2 \end{pmatrix} \begin{pmatrix}x-\mu_x \\ y-\mu_y\end{pmatrix}\\ &= \frac{1}{1-\rho^2}\left[ \left(\frac{x-\mu_x}{\sigma_x}\right)^2 +\left(\frac{y-\mu_y}{\sigma_y}\right)^2 -\frac{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x\sigma_y} \right]. \end{align*} However, equations~$\eqref{equ:MVN_pdf}$ and~$\eqref{equ:quad_form}$ are \hi{much} more important than the expansions just given. It can be proved that the function $f_{X,Y}(x,y)$ in equation~$\eqref{equ:MVN_pdf}$ integrates to $1$ and therefore defines a valid joint PDF. The proof is omitted. Remarks: \begin{enumerate} \item The vector $\vec{\mu}=(\mu_x,\mu_y)^T$ is called the \hi{mean vector} and the matrix $\Sigma$ is called the \hi{covariance} matrix (or sometimes \hi{variance-covariance} matrix). \item Functions of the form $F(\vec{x})=\vec{x}^T\Sigma^{-1} \vec{x}$ are called \hi{quadratic forms}. Quadratic forms are functions $\R^n\rightarrow \R$ which satisfy certain properties. They crop up in several areas of mathematics and statistics. \item The matrix $\Sigma$ and its inverse $\Sigma^{-1}$ are \hi{positive definite}. A matrix $A$ is positive definite if \begin{equation*} \vec{x}^TA\vec{x}\geq 0 \end{equation*} for all non-zero vectors $\vec{x}$. \item It follows that when $\mu_x=\mu_y=0$, $Q(x,y)$ is of positive definite quadratic form. \end{enumerate} Over the next two pages we will visualize the bivariate Normal distribution using 3-D plots and contour plots. \newpage \begin{center} \scalebox{1.2}{\includegraphics{bivariateNormalA}} \end{center} \newpage \begin{center} \scalebox{1.2}{\includegraphics{bivariateNormalB}} \end{center} \newpage The previous two pages show contour and surface plots for the density function on the $xy$-plane for different parameter values. Comments: \begin{enumerate} \item $Q(x,y)\geq 0$ with equality only when $\vec{x}=\vec{\mu}$. It follows that the density function has its mode at $\vec{x}=\vec{\mu}$. The above plots all have $\mu_x=\mu_y=0$. \item Changing the values of $\mu_x,\mu_y$ does not change the shape of the plots, but corresponds to a translation of the $xy$-plane i.e. changing $\mu_x,\mu_y$ just shifts the contours / surface to a new mode position. \item The contours of equal density are \hi{circular} when $\sigma_x=\sigma_y$ and $\rho=0$ and \hi{elliptical} when $\sigma_x\neq\sigma_y$ or $\rho\neq 0$. \item $\sigma_x$ and $\sigma_y$ control the extent to which the distribution is \hi{dispersed}. \item The parameter $\rho$ is the \hi{correlation} of $X,Y$ i.e.~$\cor{X,Y}=\rho$. (The proof is omitted.) Thus for non-zero $\rho$, the contours are at an angle to the axes. \end{enumerate} When $(X,Y)$ has a bivariate normal distribution then given a general set $A\subseteq\R^2$, calculation of the probability $\pr{(X,Y)\in A}$ requires a computer. However, the marginal and conditional distributions of the bivariate normal distribution are particularly simple -- one of the reasons the distribution is so useful in statistics. \vspace{1ex} \fbox{% \begin{minipage}{\textwidth} Suppose $(X,Y)^T\sim N_2(\vec{\mu},\Sigma)$. The following results hold. \vspace{2mm} \begin{enumerate} \item The marginal distributions are normal: \begin{equation*} X\sim N(\mu_x,\sigma_x^2)\quad\text{and} \quad Y\sim N(\mu_y,\sigma_y^2). \end{equation*} \item The conditional distributions are normal: \begin{align*} X|Y=y\;&\sim\; N(\mu_x+\rho\frac{\sigma_x}{\sigma_y}(y-\mu_y),\sigma_x^2(1-\rho^2))\quad\text{and}\\ Y|X=x\;&\sim\; N(\mu_y+\rho\frac{\sigma_y}{\sigma_x}(x-\mu_x),\sigma_y^2(1-\rho^2)). \end{align*} \item When $\rho=0$, $X$ and $Y$ are independent. \item Linear combinations of $X$ and $Y$ are also normally distributed: \begin{equation*} aX+bY \sim N(a\mu_x+b\mu_y,a^2\sigma_x^2+b^2\sigma_y^2+2ab\rho\sigma_x\sigma_y) \end{equation*} where $a,b$ are constants. Thus a \hi{linear combination} of normal random variables is also normally distributed itself. \end{enumerate} \end{minipage}}\\ Proofs are omitted. \subsubsection{Simulating bivariate normal random variables} R has built-in functions to simulate bivariate and multivariate normal random variables. However, we can simulate using the strategy introduced in section~\ref{sec:bivart_sim}: first simulate from the marginal for $X$ (or marginal for $Y$ if preferred), and then simulate $Y$ from the conditional $Y|X=x$. \begin{example} Suppose $(X,Y)^T\sim N_2(\vec{\mu},\Sigma)$ where $\mu_x=2$, $\mu_y=3$, $\sigma_x=1$, $\sigma_y=1$ and $\rho=0.5$. Simulate a sample of size $500$ from this distribution and draw a scatter plot. Use simulation to find $\pr{X^2+Y^2<9}$. \end{example} The marginal distribution of $X$ is $X\sim N(2,1^2)$. Using the formula for the conditional \begin{align*} Y|X=x\;&\sim\; N(\mu_y+\rho\frac{\sigma_y}{\sigma_x}(x-\mu_x),\sigma_y^2(1-\rho^2))\\ &\sim\; N(3+0.5(x-2), 0.75). \end{align*} To sample: \begin{lstlisting} npts = 500 x = rnorm(npts, mean = 2, sd = 1) y = rnorm(npts, mean = 3+0.5*(x-2), sd = sqrt(0.75)) \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{fig5-1} \end{center} To find $\pr{X^2+Y^2<9}$ approximately, count the number of points in the region: \begin{lstlisting} npts = 10000 x = rnorm(npts, mean = 2, sd = 1) y = rnorm(npts, mean = 3+0.5*(x-2), sd = sqrt(0.75)) f = x^2+y^2 sum(f<9)/npts [1] 0.2877 \end{lstlisting} \vspace{-0.5cm} \subsubsection{The multivariate normal distribution} The multivariate normal distribution is a distribution on vectors in $\R^n$. Suppose that $\vec{X}$ is a random vector with $n$ entries, i.e.\ $\vec{X}=(X_1,\ldots,X_n)^T$. Then \begin{equation*} \vec{X}\sim N_n(\vec{\mu},\Sigma) \end{equation*} if $X_1,\ldots,X_n$ have joint PDF given by \begin{equation*} f_{\vec{X}}(\vec{x}) = \frac{1}{2\pi\sqrt{\det\Sigma}}\exp\left(-\frac{1}{2} Q(\vec{x}) \right) \end{equation*} where \begin{equation*} Q(\vec{x}) = (\vec{x}-\vec{\mu})^T \Sigma^{-1} (\vec{x}-\vec{\mu}). \end{equation*} This definition makes sense for any column vector $\vec{\mu}\in\R^n$ and any \hi{positive definite} $n\times n$ matrix $\Sigma$. Some remarks: \begin{enumerate} \item The vector $\vec{\mu}$ is the \hi{mean} of the distribution and $\Sigma$ is called the \hi{covariance} matrix. \item All the marginal distributions of $\vec{X}$ are normal. (We do not specify their parameters here, however). \item Similarly, all the conditional distributions of $\vec{X}$ are normal. (Again, we do not specify the parameters of these distributions here). \end{enumerate} \begin{example} Suppose $\vec{X}\sim N_3(0,\Sigma)$ where \begin{equation*} \Sigma = \begin{bmatrix} 3 & 2 & 1 \\ 2 & 2 & 1 \\ 1 & 1 & 3 \end{bmatrix}. \end{equation*} Use R to sample from this distribution. Create a scatterplot of $(X, Y)$ samples and compute the correlation. \end{example} In R, you use the \texttt{mvrnorm} function to sample: \begin{lstlisting} library(MASS) npts = 500 mu = c(0, 0, 0) sigma = matrix(data = c(3,2,1,2,2,1,1,1,3), nrow = 3, byrow = TRUE) x = mvrnorm(npts, mu, sigma) \end{lstlisting} \vspace{-0.5cm} Now \texttt{x} has dimension \texttt{npts}$\times 3$, so column 1 is $X$ and column $2$ is $Y$: \begin{lstlisting} plot(x[,1], x[,2], xlab = 'x', ylab = 'y') \end{lstlisting} \vspace{-1.3cm} \begin{center} \includegraphics[width=6.5cm]{fig5-2} \end{center} \vspace{-0.75cm} The sample correlation is: \begin{lstlisting} cor(x[,1],x[,2]) [1] 0.812629 \end{lstlisting} \end{document} \newpage \section{Samples and populations} In this section we consider the fundamental statistical concepts of \hi{populations} and \hi{samples}. We investigate the properties of samples, particularly as the sample size gets large, via computer simulation. This material is intended to reinforce the ideas studied in MAS2901. \begin{example} A factory has been depositing mercuric waste into a river near its estuary. It is feared that this waste may have caused the mercury content in the prawns there to have risen to a level which may be toxic for human consumption. The maximum permissible level is one part per million and the factory claims the average mercury content in the prawns does not exceed this level. This claim is to be tested. It is obviously impractical to measure the mercury content of all the prawns. A random sample of $20$ prawns is collected. The concentration of mercury (parts per million) in each of the $20$ prawns is measured, giving the following results: the observed sample mean is $\bar{x}=1.2$, and the observed sample standard deviation $=s=0.5$. Given these measurements do we believe the factory's claim? \end{example} \subsection{Populations} Suppose we measure some random quantity $X$: $X$ can adopt a range of possible values, and some values are more likely than others. This is, of course, the familiar idea of a \hi{distribution} for $X$. Usually we do not know this distribution exactly, but we want to learn something about it. The unknown distribution is called the \hi{population distribution}. Referring to the prawn example: \begin{itemize} \item the population consists of the prawns in the estuary; \item the random quantity $X$ is the mercury concentration in a randomly selected prawn; and \item the population distribution is the distribution of $X$. \end{itemize} We are usually interested in key properties of the population distribution such as: \begin{itemize} \item the expectation of $X$ -- usually called the \hi{population mean}; \item the variance of $X$ -- usually called the \hi{population variance}; or \item the $95$th percentile of $X$ (for example). \end{itemize} Often we make some simplifying assumptions about the population distribution. For example, we might assume: \begin{enumerate}[(a)] \item $X$ is normally distributed with unknown mean and variance; \item $X$ is exponentially distributed with rate parameter $\lambda$, where $\lambda$ is unknown but lies on the interval $(0,1)$; \item $X$ is normally distributed with unknown mean and variance $\sigma^2=5$. \end{enumerate} A set of assumptions like this is referred to as a \hi{model}. Typically a model involves some unknown parameters. In the examples above we have the unknown rate $\lambda$ in (b), and unknown mean $\mu$ in (c). In some situations -- usually rather artificial ones -- we know the population distribution exactly. For example: \begin{itemize} \item let $X$ be the score obtained from rolling a fair die; or \item let $X$ be the number on a card drawn at random from a full deck. (Assume Jack, Queen, King numbered 11, 12, 13 respectively.) \end{itemize} \subsection{Samples} Typically we do not know everything about the population distribution. We learn about the population distribution by drawing a \hi{sample}. Referring to the prawn example: we cannot measure the mercury concentration in every prawn in the estuary, so we sample 20 prawns at random and analyse those. A sample of size $n$ corresponds to taking $n$ independent measurements from the distribution. Each measurement is itself a random variable with the same distribution as $X$: the sample measurements denoted $X_1,X_2,\ldots,X_n$. So in our example, $X_1$ represents the mercury concentration in the first prawn collected -- and this is a random quantity. The actual measurements obtained are denoted $x_1,x_2,\ldots,x_n$. The distinction between the population distribution and how we learn about the population from limited samples is probably the most \hi{important concept} in statistics. \subsection{Estimators} Suppose we wish to learn about some aspect of the population distribution: the population mean or the population variance for example. To do this we construct an \hi{estimator} for the quantity of interest. For example, if we wish to learn about the population mean, a good estimator is the sample mean \begin{equation*} \bar{X}=\frac{1}{n}\sum_{i=1}^n X_i. \end{equation*} (We won't address what a `good' estimator might be in this module.) Formally, an estimator is defined to be some \hi{function} of the sample, i.e.\ an estimator $S$ is a random variable for the form \begin{equation*} S = g(X_1,X_2,\ldots,X_n) \end{equation*} for some function $g$. It follows that every estimator is a \hi{random variable}. When we observe some measurements $X_1=x_1, \ldots, X_n=x_n$ then we can compute the corresponding value of the estimator: \begin{equation*} s = g(x_1,x_2,\ldots,x_n). \end{equation*} This value is called an \hi{estimate}. To summarize: an estimator $S$ is a random quantity which is a function of the sample; whereas an estimate is an observed value of that random quantity. Since any estimator $S$ is a random variable it makes sense to talk about its \hi{distribution}. The rest of this section is concerned with studying the distribution of some simple estimators. \begin{example} Suppose the population distribution is normal, and we wish to estimate the population mean. Suppose the sample size is $n=4$ and our estimator is $\bar{X}=(X_1+X_2+X_3+X_4)/4$. What is the distribution of $\bar{X}$ when the population distribution is $N(170,20^2)$? \end{example} Simulation can be used to answer this question approximately: we can simulate lots of samples and look at the resulting distribution. \begin{lstlisting} simulate.xbar = function(n) { xbar = vector(mode = 'numeric', length = n) for (i in 1:n) { x = rnorm(4, 170, 20) # Generate a sample of size 4 xbar[i] = 0.25*sum(x) } xbar } \end{lstlisting} Histogram of $\bar{X}$: \begin{lstlisting} test = simulate.xbar(500) hist(test, xlab = 'sample mean', ylab = 'frequency') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{fig6-1} \end{center} \newpage \begin{example} Suppose the population distribution is normal, and we wish to estimate the 90th percentile using a sample of size 10. A sensible estimator is to define $S$ to be the second largest value in the sample (i.e.\ the 9th value when the samples are ordered from smallest to largest). What is the distribution of $S$ when the population distribution is $N(0,1)$? \end{example} Code to simulate the distribution of $S$: \begin{lstlisting} simulate.percentile = function(n) { s = vector(mode = 'numeric', length = n) for (i in 1:n) { x = rnorm(10, 0, 1) # Generate a sample of size 10 x = sort(x) s[i] = x[9] # Get 9th value on sorted list } s } \end{lstlisting} Histogram of $S$: \begin{lstlisting} S = simulate.percentile(500) hist(S, xlab = 's', ylab = 'frequency', main = '') \end{lstlisting} \vspace{-1cm} \begin{center} \includegraphics[width=8cm]{fig6-2} \end{center} \subsection{Simulation experiment} Consider the following two examples for the density of the population distribution. For each example, decide which histogram on the slides (A, B, C or D) is most likely to represent the distribution of the sample mean $\bar{X}$ when the sample size is 10. \newpage \begin{example} Suppose $X$ (the population quantity) has the following PDF: \begin{center} \includegraphics[width=7cm]{fig6-3} \end{center} \end{example} \hi{Most likely histogram:} \vspace{0.5cm} \begin{example}\label{ex:bimodal} Suppose $X$ (the population quantity) has the following PDF: \begin{center} \includegraphics[width=8cm]{fig6-4} \end{center} \end{example} \hi{Most likely histogram:} Conclusion: \begin{itemize} \item The sample mean is distributed around the population mean. \item The distribution of sample mean values `forgets' the underlying shape of the population distribution. \item As $n$ increases we expect the distribution of $\bar{X}$ to become more clustered around the true value. \end{itemize} \subsection{The Central Limit Theorem} The Central Limit Theorem is the most important theorem in statistics. It describes the distribution of the sample mean $\bar{X}$ as the sample size increases. \fbox{% \begin{minipage}{\textwidth} Suppose $X_1,X_2,\ldots,X_n$ are independent and \hi{identically distributed} random variables with common mean $\mu$ and variance $\sigma^2$ which are both finite. Define \begin{equation*} Z = \frac{\bar{X}-\mu}{\sigma/\sqrt{n}}. \end{equation*} Then as $n\rightarrow\infty$ the distribution of $Z$ tends to $N(0,1)$. \end{minipage}}\\ While this results describes the approximate distribution of the sample mean as $n$ gets large, it can also be adapted and extended to understand the asymptotic behaviour of other estimators. Remarks: \begin{enumerate} \item Equivalently, the Central Limit Theorem shows that for large $n$, $\bar{X}\sim N(\mu,\sigma^2/n)$ approximately. \item It follows that $\bar{X}$ is distributed somewhere around the population mean $\mu$. \item As $n$ increases, the variance of the distribution of $\bar{X}$ decreases. In other words, $\bar{X}$ is more likely to lie closer to $\mu$ -- ``larger samples give better estimates''. \item The underlying population distribution is forgotten as $n$ increases, as $\bar{X}$ has a normal distribution. \end{enumerate} \begin{example} Suppose $X$ has the bimodal distribution considered in example~\ref{ex:bimodal}: \begin{center} \includegraphics[width=8cm]{fig6-4} \end{center} This is a normal mixture distribution, and the following function can be used to sample values $\bar{X}$ from the distribution when the sample size is $n$. \newpage \begin{lstlisting} simulate.bimod = function(k, n) { # Generate k samples of size n s = vector(mode = 'numeric', length = k) for (i in 1:k) { u = rnorm(n, 3, 0.6) v = rnorm(n, 7, 0.6) r = runif(n) x = c(u[r>0.5], v[r<=0.5]) s[i] = mean(x) } s } \end{lstlisting} Use this function to study the distribution of $\bar{X}$ for sample sizes of $n=2$, $n=5$ and $n=10$. The population mean is $\mu=5$ and $\sigma^2=4.3$. \end{example} Generate samples: \begin{lstlisting} samples.size.2 = simulate.bimod(1000, 2) samples.size.5 = simulate.bimod(1000, 5) samples.size.10 = simulate.bimod(1000, 10) \end{lstlisting} Histograms: \begin{center} \includegraphics[width=5cm]{fig6-5} \includegraphics[width=5cm]{fig6-6} \includegraphics[width=5cm]{fig6-7} \end{center} We can compare the observed means and variances to the values predicted from the Central Limit Theorem: \begin{center} \begin{tabular}{|c|cc|cc|} \hline & & & & \\ Sample size $n$ & $\mu$ & $\sigma^2/n$ & Observed mean of $\bar{X}$& Observed variance of $\bar{X}$ \\ & & & & \\ \hline & & & & \\ 2 & 5.0 & 2.15 & $4.94$ & $2.27$\\ 5 & 5.0 & 0.86 & $4.98$ & $0.862$\\ 10 & 5.0 & 0.43 & $4.96$ & $0.443$\\ \hline \end{tabular} \end{center} \newpage The observed mean and variance of $\bar{X}$ for $n=2$ are calculated by: \begin{lstlisting} mean(samples.size.2) [1] 4.939386 var(samples.size.2) [1] 2.270115 \end{lstlisting} Further examples will be studied in the practicals. \newpage %%%Practicals \section*{Aims} This introductory practical is split into two parts: \begin{itemize} \item In the first hour or so, we will focus on the basics of R, revising some of the key concepts from last year's MAS1608 material \item In the second hour we will focus on material from this week's lecture -- using R to simulate from standard probability distributions, as well as more complex models including mixture distributions \end{itemize} \vspace{0.25cm} \noindent Please read through all of this handout \textit{very carefully}, and follow the instructions \textit{meticulously}. Whenever you are typing commands into R, make sure you do so \textit{extremely carefully}. If you get an error message, try typing the code in again, just to make sure you've not made a spelling mistake or other simple typo. A common error is to forget that R is case sensitive, so use upper/lower cases exactly when you're told to do so. \vspace{0.25cm} \noindent For all R practicals, make sure you've got your lecture notes handy. \section*{Getting started} \noindent Log-in to a PC in the usual way. Now open RStudio. If you can't remember how to do this, click on the ``Start'' button in the bottom left-hand-corner and search for ``RStudio''. Once you've opened RStudio you should see four panels; if there is no text editor (usually the top-left panel), click on the green button with a white plus symbol and open a new R script. Remember -- to execute a line of code in the text editor, hit \texttt{Ctrl}-\texttt{Enter}. \vspace{0.25cm} \noindent It is good practice to always use the text editor, rather than type into the console window (bottom- left panel) directly. This is especially true in future practicals when you may be building large functions; it is much easier to de-bug and edit code in the text editor than in the console window. \vspace{0.25cm} \noindent \textbf{Hint: To commit several lines of code in one go, highlight your block of code in the text editor and hit \texttt{Ctrl}-\texttt{Enter} just once!} \vspace{0.25cm} \noindent Now you've got RStudio open, you're ready to go! Try out \textit{all} of the following exercises. \newpage \section*{1 R Refresher} \subsection*{Querying databases; plots and summaries} \noindent You might remember the movies database from MAS1608. This database is available from the MAS2602 webpage -- you can load it into R by typing the following code: \begin{lstlisting}[frame=tb, label=2, caption = Loading the movies database] movies = read.table("http://www.mas.ncl.ac.uk/~nlf8/movies.txt", header=T) attach(movies) \end{lstlisting} \begin{itemize} \item [(a)] Use the \texttt{head} command to look at the first six rows of the database. How many variables are there in the dataset? \item [(b)] Use the \texttt{dim} command to find out the dimensions of the database. How many movies are there? \item [(c)] How many people submitted votes for the films \textit{"AVP: Alien Vs. Predator"}, \textit{"Midnight in the Garden of Good and Evil"} and \textit{"Zoolander"}? \textcolor{gray}{(Hint: the \texttt{which} and \texttt{tail} commands might come in useful here)} \item [(d)] What is the mean number of votes per film across all films in the database? What is the standard deviation? \item [(e)] How many films are classified as "Action"? How many are classified as "Romance"? How many "Short" films are there in the database? \textcolor{gray}{(Hint: use the \texttt{sum} command)} \item [(f)] Type the following code: \begin{lstlisting}[frame = tb, label=3, caption = Querying the movies database] Budgets_known = Budget[Budget!=(-1)] par(mfrow=c(2,2)) hist(Budgets_known) abline(v = mean(Budgets_known), col = "red") boxplot(Rating, main = "Boxplot for movie ratings") barplot(table(mpaa), xlab = "MPAA rating", col = "misty rose") apply(movies[ ,3:6], 2, mean) Length2 = Length[Year>=1970 & Year<=2000] \end{lstlisting} Using the code in Listing 2: \begin{itemize} \item [(i)] Can you figure out what is going on in line 1? \item [(ii)] What is the purpose of line 2? \item [(iii)] Re-produce the plots, but this time colour the histogram "forest green", change its title, and draw a vertical line at the median budget; colour the boxplot "wheat"; and add the word "frequency" to the y-axis of the bar chart. \item [(iv)] Produce a scatterplot of Rating against Votes, and explain what you see. \item [(v)] What is going on in line 7? \item [(vi)] Find the median and interquartile range of movie lengths, for movies made between the years 1970 and 2000 (inclusive) \textcolor{gray}{(Hint: the command \texttt{quantile} might be useful here)} \item [(vii)] Write a line of code that will extract the movie names, and the years in which the movies were made, into a new dataframe - and store this new dataframe in \texttt{movies2}. \end{itemize} \end{itemize} \subsection*{Functions, for loops and if statements} Suppose the random variable $Z$ follows a standard Normal distribution, that is, $Z \sim N(0,1)$. Then $Z$ has probability density function given by \[ f_{Z}(z) = \frac{1}{\sqrt{2 \pi}} e^{-\frac{1}{2}z^{2}}, \quad z \in \mathbb{R}. \] Type the following code into R: \begin{lstlisting}[frame = tb, label=4, caption = Drawing the standard Normal distribution with simulation grid] z = seq(-4, 4, 0.001) d = 1/(sqrt(2*pi))*exp(-(z^2)/2) plot(z, d, xlab = "z", ylab = "density", type = "l") abline(v = -2, lty = 2, col = "red") abline(v = 1, lty = 2, col = "red") abline(h = 0.4, lty = 2, col = "red") \end{lstlisting} We will now use Monte Carlo integration to estimate $P(-2 < Z < 1)$ by following the steps below: \begin{itemize} \item [1.] Throw randomly generated points at a "sampling grid" that covers the region of interest (see red dashed lines on your plot) \item [2.] Calculate the proportion of points $P$ that lie under the density curve for $Z$ \item [3.] Find \[ P(-20. \] Use the following code to superimpose your fitted PDF on the histogram: \begin{lstlisting}[frame = tb, label=7, caption = Overlaying fitted exponential PDF] x = seq(0, 80, 0.001) d = lambda*exp(-lambda*x) lines(x, d, type="l", col="red") \end{lstlisting} \item [(d)] Does your fitted exponential distribution provide a good fit to the wind speed data? \item [(e)] During a hurricane, observed wind speeds exceed 70 miles per hour. Use R to find the probability, according to your fitted model, of wind speeds exceeding 70 miles per hour. \item [(f)] How does your calculation in part (e) change if a new data point of 72 miles per hour is included in the analysis? \end{itemize} \item [3.] Use the \texttt{pnorm} command in R to confirm your Monte Carlo integration approximation to $P(-22,\\ 0& \text{otherwise}. \end{cases} \end{equation*} \begin{itemize} \item [(a)] Note that your are given the PDF $f_{X}$, not the CDF $F_{X}$. Thus, you should obtain the CDF by integrating, i.e. find \[ \int_{2}^{x} \left(\frac{2}{t}\right)^3 dt. \] \item [(b)] Find the inverse of the CDF, that is, $F_{X}^{-1}$. \item [(c)] Write an R function called \texttt{inv.cdf} that takes a single input variable: \texttt{N}, the number of realisations to generate. Within this function, you should implement the inverse CDF method to generate values from $X$ by applying your answer for $F_{X}^{-1}$ in part (b) to observations generated randomly from a $U(0,1)$ distribution. \item [(d)] Execute the following code to produce a histogram of 1000 realisations from $X$, on a density scale, and to superimpose the PDF of $X$: \begin{lstlisting} X = inv.cdf(1000) hist(X, probability = TRUE) curve((2/x)^3, from = 0, to = 100, col = "red", add = TRUE) \end{lstlisting} \end{itemize} \vspace{1.5cm} \item [2.] It is suggested that the survival times of patients ($X$ years) receiving a new cancer drug will follow a Weibull distribution with CDF \[ F_{X}(x; \lambda, \kappa) = 1-e^{-(x/\lambda)^{\kappa}}, \quad x \geq 0, \quad \lambda, \kappa>0, \] where $\lambda$ is referred to as the scale parameter and $\kappa$ is the shape parameter. \begin{itemize} \item [(a)] Write an R function that uses the inverse CDF method to generate a sample of size \texttt{N} from the Weibull distribution. Your function should have \textit{three} input variables: \texttt{N}, \texttt{lambda} and \texttt{kappa}. \item [(b)] Past studies from a similar drug indicate that $\lambda = 2$ and $\kappa = 1.5$. Use these values, and your function in part (a), to simulate 10,000 realisations from $X$, and plot a histogram of these values. \item [(c)] Use your simulated values in part (b) to estimate \begin{itemize} \item [(i)] the mean survival time of patients receiving this new cancer drug; \item [(ii)] the probability that a randomly chosen patient receiving this new drug will survive beyond five years. \end{itemize} \item [(d)] Find out if there is a built-in function in R to generate samples from a Weibull distribution. If there is, use this to generate a sample from a Weibull distribution with $\lambda = 2$ and $\kappa=1.5$, and re-produce your answers to parts (b) and (c). \end{itemize} \newpage \item [3.] In a hurricane-prone region, wind speed excesses ($X$ metres per second) over some high threshold are modelled with a \textit{Generalised Pareto Distribution}, with PDF \[ f_{X}(x; \sigma) = \sigma^{-1}e^{-x/\sigma} \quad x\geq 0, \] and CDF \[ F_{X}(x; \sigma) = 1-e^{-x/\sigma}, \quad x \geq 0. \] Hourly maximum wind speed excesses during a storm, for a 24 hour period, resulted in a sample mean of $11.6\textrm{ ms}^{-1}$, with lower and upper quartiles of $3.4\textrm{ ms}^{-1}$ and $15.9\textrm{ ms}^{-1}$ respectively. From these wind speed excesses, the maximum likelihood estimate for $\sigma$ can be shown to be $\hat{\sigma} = 12$. \begin{itemize} \item [(a)] Write an R function to return realisations from the Generalised Pareto Distribution. \item [(b)] Use your function in part (a) to obtain a sample of size 1000 from the Generalised Pareto Distribution, and store this sample in the vector \texttt{X}. \item [(c)] Produce a histogram of your sample in \texttt{X}, using the argument \texttt{probability = TRUE}, and overlay the PDF as provided at the start of this question. Make sure you give your histogram a suitable title, axis labels etc., and experiment with different colours and/or line types, as demonstrated in lectures. \item [(d)] Produce summary statistics for your sample of realisations from this model. How does your simulated sample compare to the real life data? \item [(e)] Using your simulated series, how often might we expect hourly wind speeds to exceed $50\textrm{ ms}^{-1}$? \item [(f)] Use R to directly generate 1000 realisations from an exponential distribution with rate $1/\hat{\sigma}$, and explore this simulated series. How does this compare to your series stored in \texttt{X}? \end{itemize} \end{itemize} \newpage \section*{Transformation of random variables} \begin{enumerate} \item Suppose $X\sim U(1,5)$, and that we consider two transformations of $X$: $Y_{1} = \sqrt{X}$ and $Y_{2} = \log{X}$. \begin{itemize} \item [(a)] Use R to simulate 1000 values on each of $Y_{1}$ and $Y_{2}$, and store the results in the vectors \texttt{Y1} and \texttt{Y2} respectively. \item [(b)] Derive the probability density function (PDF) for each of $Y_{1}$ and $Y_{2}$. \item [(c)] On a $1\times 2$ panel of plots in R (i.e. 1 row, 2 columns), produce a histogram of your samples from $Y_1$ and $Y_2$. On each histogram, overlay the PDF you derived in part (b) to confirm your simulations. Make sure you give your histograms suitable titles, axis labels etc; produce the histograms on the density scale by using the argument \texttt{probability=TRUE}; and experiment with different colours and/or line types, as demonstrated in lectures. [Hint: You will need to use the \texttt{par(mfrow = c(1, 2))} command to partition the plotting space, and you may need to experiment with the \texttt{curve} command for overlaying your derived densities.] \end{itemize} \vspace{1cm} \item Suppose $Z \sim N(0, 1)$, and that we consider the transformation $Y = Z^{2}$. \begin{itemize} \item [(a)] Use R to simulate 1000 values on $Y$, and store the results in the vector \texttt{Y}. \item [(b)] Produce a histogram of your sample from $Y$, and comment on the shape of the distribution. Make sure you use the argument \texttt{probability = TRUE} when constructing your histogram. \item [(c)] The \textit{chi-squared distribution on 1 degree of freedom} has PDF \[ f_{X}(x) = \frac{1}{2^{1/2}\Gamma(1/2)}x^{-1/2}e^{-x/2}, \quad x \geq 0. \] Superimpose the PDF of the chi-squared distribution onto your histogram of realised values of $Y$, and comment. [Hint: the Gamma function $\Gamma(\cdot)$ can be evaluated using the built-in R command \texttt{gamma}] \item [(d)] Use R to find the probability $\pr{Y>3.84}$ by considering the proportion of your simulated values in \texttt{Y} that exceed 3.84. \item [(e)] R has built-in functions for random simulation, and computing cumulative probabilities, from the chi-squared distribution. Find out what these are so you can simulate 1000 realisations from the chi-squared distribution directly, and so you can find a more accurate answer to the probability calculation in part (d). [Hint: recall the prefixes \texttt{r} and \texttt{p} for random generation, and obtaining cumulative probabilities, from standard models such as the Normal distribution; i.e. \texttt{rnorm}, \texttt{pnorm}. What do you \textit{think} you might use for the chi-squared distribution?] \end{itemize} \vspace{1cm} \item 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. Some pointers: \begin{itemize} \item [--] Produce a histogram of a sample of realisations from $Y$, and find some numerical summaries \item [--] $Y$ actually follows \textit{Cauchy} distribution. Is there a built-in command in R to sample directly from this distribution? \end{itemize} \end{enumerate} \newpage \section*{Bivariate distributions} \begin{enumerate} \item Suppose that $X\sim N(100,25)$ and that \[ f_{Y|X}(Y|X=x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(y-x)^2}. \] Write an R function to sample realizations from $(X,Y)$, and explore the dependence between $X$ and $Y$. \textcolor{gray}{[Hint: Do you recognize the distribution of $Y|X$? See Examples 4.1 and 4.2 in the lecture notes for more help]} \vspace{1cm} \item Suppose \[ \left(\begin{array}{c} X \\ Y \\ \end{array}\right) \sim N_{2}\left(\left(\begin{array}{c} 1 \\ 3 \\ \end{array}\right), \left(\begin{array}{c c} 2 &-0.8\sqrt{2} \\ -0.8\sqrt{2} & 1 \end{array}\right)\right). \] \begin{itemize} \item [(a)] Write down the values of $\mu_{x}$, $\mu_{y}$, $\sigma_{x}$ and $\sigma_{y}$. \item [(b)] Find the correlation $\rho$. \item [(c)] Find the conditional distribution for $Y|X=x$, and use this to simulate a sample of 1000 from the bivariate distibution for $(X,Y)^{T}$. Produce a plot of $X$ versus $Y$, and use this to confirm your answer to part (b). \textcolor{gray}{[Hint: See result 2 on page 25 of the lecture notes, and Example 4.3]} \item [(d)] We will now repeat the simulation exercise in part (c) using the \texttt{MASS} package. Set up your mean vector and variance matrix with the following code: \begin{lstlisting} mu=c(1,3) sigma=matrix(data=c(2,-0.8*sqrt(2),-0.8*sqrt(2),1),nrow=2,byrow=TRUE) \end{lstlisting} \vspace{-0.6cm} Now install the \texttt{MASS} package, which can be used to simulate observations from a multivariate normal distribution; then explore the \texttt{mvrnorm} function: \begin{lstlisting} library(MASS) ?mvrnorm \end{lstlisting} \vspace{-0.6cm} Simulate a sample of 1000 pairs from the bivariate distribution for $(X,Y)^{T}$ using \texttt{mvrnorm}, and produce a plot of $X$ versus $Y$. \textcolor{gray}{[Hint: See Example 4.4 in the lecture notes]} \end{itemize} \end{enumerate} \newpage \section*{Assignment} This assignment is worth 10\% of the MAS2602 module marks. Solutions should be handed in to the General Office n the usual way by \textcolor{red}{\bf 3pm on Thursday 21st November}. For information about the late work policy, see Blackboard or the preliminary module information given in the module notes. Your work should be typed up and any necessary graphs or code included in your typed document. Marks will be awarded for style and presentation. \subsection*{Question 1 [Total 25 marks]} In this question, you are required to use Monte Carlo Integration to approximate $\displaystyle \int_{0}^{1}f(x) dx$, where \[ f(x) = a b x^{a-1}(1-x^a)^{b-1} \] for $0 \le x \le 1$ and where $a>0$ and $b>0$. \begin{itemize} \item [(a)] Each student will need to generate their own unique values for $a$ and $b$. To generate \textit{your} values for $a$ and $b$, type the following code into R: \begin{lstlisting} set.seed(LOGIN_ID) a = sample(3:14, 1); b = sample(3:14, 1) \end{lstlisting} \vspace{-0.75cm} where \lstinline{LOGIN_ID} is your computer ID without the `a' or `b' at the start. In your solutions, write down your values for $a$ and $b$ and the resulting integrand. For example, if your generated values for $a$ and $b$ were 10 and 12 respectively, you would write down: \[ a = 10; \quad b = 12; \quad f(x) = 120 x^{9}(1-x^{10})^{11}. \] \item [(b)] Use R to produce a plot of $f(x)$. Within the \texttt{plot} command, use the argument \texttt{xlim = c(-1, 2)}. Use the \texttt{abline} command to indicate a suitable \textit{simulation grid} on your plot. Include this plot in your solutions, with appropriately labelled axes and title. \textbf{[7 marks]} \item [(c)] Using a \texttt{for} loop over $N=10^6$ replications, estimate your integral using Monte Carlo Integration. Include your R code in your solutions, which should be nicely formatted/indented. Where necessary, include comments in your code so it is clear what you are trying to do. In your solutions, make it obvious what your estimated integral is. \textbf{[18 marks]} \end{itemize} \vspace{2cm} \begin{center} \textit{\bf Turn over for questions 2 and 3} \end{center} \newpage \subsection*{Question 2 [Total 25 marks]} Suppose that $X$ is a random variable with PDF \begin{equation*} f_X(x) = \frac{x}{\sigma^2}e^{-x^2/(2\sigma^2)}\quad\text{for\ }x>0 \end{equation*} where $\sigma$ is a fixed constant. \begin{enumerate}[(a)] %\item By considering the derivative of $g(x)=e^{-x^2/(2\sigma^2)}$, find the CDF $F_X(x)$ of $X$. \item Find the CDF $F_X(x)$ of $X$. Make sure you include some working, in your solutions. \textbf{[4 marks]} \item Write a function to generate simulated values of $X$ using the inverse CDF method. Include the R code for this function in your solutions. \textbf{[8 marks]} \item Let $d$ be the last digit of your student ID. Let $\sigma = 1+(d/2.5)$. Write down your value of $d$, and plot a density histogram of a sample of $X$ of size 10,000 using your value of $\sigma$. Include this histogram in your solutions; again, make sure axes are appropriately labelled. \textbf{[5 marks]} \item Check that your function has worked by superimposing a curve on your histogram showing the PDF. \textbf{[2 marks]} \item The distribution you have simulated from is often used in the physical sciences to model wind speeds and wave heights; it is also used to model noise variation in magnetic resonance imaging (MRI). Can you find the name of this distribution? \textbf{[2 marks]} \item Use your simulated valued of $X$ to verify the expectation and variance of the distribution you identified in part (e). \textbf{[4 marks]} \end{enumerate} \vspace{1.25cm} \subsection*{Question 3 [Total 25 marks]} \begin{enumerate}[(a)] \item Write an R function to sample from the log-normal distribution $LN(\mu,\sigma^2)$ (see section 3.1 of the notes and example 3.3). Make sure you function takes three input arguments: the number of realisations to generate (\texttt{N}), the mean (\texttt{mu}) and the standard deviation (\texttt{sigma}). Include your R code in your solutions. \textbf{[5 marks]} \item For $\sigma=2$ generate samples from the log-normal distribution for a number of different values for $\mu$ in the range $\mu\in(0,5)$. For each value of $\mu$ compute the log of the sample mean. Plot a graph of the log of the sample mean versus $\mu$, and include this plot in your solutions. \textbf{[8 marks]} \item Repeat part (b), plotting the log of the sample median versus $\mu$. \textbf{[8 marks]} \item What can you conclude about the mean and median of the $LN(\mu,\sigma^2)$ distribution? \textbf{[4 marks]} \end{enumerate} \newpage \section*{Test-style questions} \begin{itemize} \item [\textbf{1.}] Suppose $X \sim N(0,1)$ and that $Y = X^2$. Which option, A--E, correctly specifies the probability density function of $Y$? \begin{enumerate}[A.] \item $f_Y(y)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}y}$ for $0\leq y<\infty$ \item $f_Y(y)=\frac{1}{2\pi}e^{-\frac{1}{4}y^2}$ for $0\leq y<\infty$ \item $f_Y(y)=\frac{1}{\sqrt{8\pi y}}e^{-\frac{1}{2}y}$ for $0<\leq y<\infty$ \item $f_Y(y)=\frac{1}{2^{1/2}y^{1/2}\Gamma(1/2)}e^{-\frac{1}{2}y}$ for $0\leq y <\infty$ \item $f_Y(y)=\lambda e^{(y-\lambda e^y)}$ for $-\infty8)$. \item [(d)] Use your sample in part (b) to estimate the correlation between $X$ and $Y$. How does this compare to the true correlation? \end{itemize} \item [\textbf{4.}]It is suggested that the survival times of patients ($X$ years) receiving a new cancer drug will follow a Weibull distribution with CDF \[ F_{X}(x; \lambda, \kappa) = 1-e^{-(x/\lambda)^{\kappa}}, \quad x \geq 0, \quad \lambda, \kappa>0, \] where $\lambda$ is referred to as the scale parameter and $\kappa$ is the shape parameter. \begin{itemize} \item [(a)] Write an R function that uses the inverse CDF method to generate a sample of size \texttt{N} from the Weibull distribution. Your function should have \textit{three} input variables: \texttt{N}, \texttt{lambda} and \texttt{kappa}. \item [(b)] Past studies from a similar drug indicate that $\lambda = 2$ and $\kappa = 1.5$. Use these values, and your function in part (a), to simulate 10,000 realisations from $X$, and plot a histogram of these values. \end{itemize} \end{itemize} \newpage \section*{Solutions: Practical 1} \section*{R Refresher} \subsection*{Querying databases; plots and summaries} \begin{itemize} \item [(a)] 24 \item [(b)] 4847 \item [(c)] 14651; 8517; 18277 \item [(d)] 4013.53; 10620.21 \item [(e)] 927; 739; 16 \item [(f)] \begin{itemize} \item [(i)] Removing all the missing budgets \item [(ii)] Partition the plotting space into a 2 by 2 space (giving four plots) \item [(iii)] Re-do plots... \item [(iv)] Use \texttt{plot(Rating, Votes)} \item [(v)] Working out the mean for columns 3, 4, 5 and 6 \item [(vi)] Median = 97.0; IQR = 17.0 \item [(vii)] \texttt{movies2 = movies[ ,1:2]} \end{itemize} \end{itemize} \subsection*{Functions, for loops and if statements} \begin{itemize} \item [(a)] Answer should be around 0.819 \item [(b)] Need to have \texttt{MCintegrate = function(N, l, u)}, and then replace \texttt{x = runif(1, -2, 1)} with \texttt{x = runif(1, l, u)}. \end{itemize} \section*{Simulation} \begin{itemize} \item [1.] \begin{itemize} \item [(a)] \texttt{W = rpois(1000, 2.5); barplot(table(W))} \item [(b)] \begin{itemize} \item [(i)] \texttt{dnorm(3, 12, 0.01)} \item [(ii)] \texttt{pnorm(2, 12, 0.01)} \item [(iii)] \texttt{X = rbinon(1000, 12, 0.01)}; \texttt{mean(x)} (should be close to $n \times p$) ; \texttt{var(x)} (should be close to $n \times p \times (1-p)$ \end{itemize} \item [(c)] \begin{itemize} \item [(i)] No - why have such definite cutoffs, and equal probability in between? \item [(ii)] \texttt{Y = runif(1000, -15, 45)}; \texttt{mean(Y)}; \texttt{var(Y)} \item [(iii)] \texttt{hist(Y, xlab = "Minutes late")} \end{itemize} \end{itemize} \item [2.] \begin{lstlisting} wind = c(26,10,10,3,5,33,13,23,5,8,14,24,3,2,7,10,46,15,2,35) hist(wind, freq = FALSE, xlab="Wind speed (mph)", main="Histogram of wind speeds", xlim=c(0,80)) lambda = 1/mean(wind) x = seq(0, 80, 0.001) d = lambda*exp(-lambda*x) lines(x, d, type="l", col="red") length(wind[wind>70])/length(wind) \end{lstlisting} \item [3.] Use \texttt{pnorm(1)-pnorm(-2)}. \item [4.] See page 6 of the lecture notes - amend code! \end{itemize} \newpage \section*{Solutions: Practical 2} \begin{itemize} \item [1.] \begin{itemize} \item [(a)] We find \begin{eqnarray*} F_{X}(x) &=& \int_{2}^{x} \left(\frac{2}{t}\right)dt \\ & & \\ &=&2^{3}\int_{2}^{x} t^{-3} dt \\ & & \\ &=& 8\left[-\frac{1}{2t^{2}}\right]_{2}^{x} \\ & & \\ &=&1 - \frac{4}{x^{2}}, \end{eqnarray*} giving \[ F_{X}(x) = \left\{ \begin{array}{c l} 0 & x\leq 2 \\ 1-\frac{4}{x^{2}} &x > 2 \end{array} \right. \] \item [(b)] Finding the inverse - rearranging $F_{X}$ to make $x$ the subject gives \[ F_{X}^{-1}(x) = \sqrt{\frac{4}{1-x}}, \quad x>2. \] \item [(c)] R function: \begin{lstlisting} inv.cdf = function(N){ X = vector("numeric", length = N) for(i in 1:N){ U = runif(1) X[i] = sqrt(4/(1 - U))} return(X)} \end{lstlisting} \item [(d)] \textcolor{white}{Space} \begin{center} \includegraphics[width=9cm, height=7cm]{sol1_2019.pdf} \end{center} \end{itemize} \newpage \item [2.] \begin{itemize} \item [(a)] Following the same steps as in question 1 (without having to integrate, as we are given the CDF here!), gives: \begin{lstlisting} inv.cdf.weibull = function(N, lambda, kappa){ X = vector("numeric", length = N) for(i in 1:N){ U = runif(1) X[i] = lambda*((-log(1-U))^(1/kappa))} return(X)} \end{lstlisting} \item [(b)] Implementing the function in (a), with \texttt{N=10000}, \texttt{lambda=2} and \texttt{kapp=1.5} and storing the values in \texttt{X}, gives: \begin{center} \includegraphics[width=9cm, height=7cm]{sol2_2019.pdf} \end{center} \item[(c)] \begin{lstlisting} mean(X) [1] 1.800228 length(X[X > 5])/length(X) [1] 0.018 \end{lstlisting} \item [(d)] There is! Students should explore the intrinsic function \texttt{rweibull}, noting that the first parameter is the shape and the second is the scale, so: \begin{lstlisting} X2 = rweibull(10000, 1.5, 2) \end{lstlisting} \end{itemize} \newpage \item [3.] \begin{itemize} \item [(a)] Let $y = 1-e^{-x/\sigma}$. Then \begin{eqnarray*} -x/\sigma &=&\log(1-y) \\ x&=&-\sigma\log(1-y), \end{eqnarray*} and so \[ F_{X}^{-1} = -\sigma \log(1-x). \] Code: \begin{lstlisting} inv.cdf.gpd = function(N, sigma){ X = vector("numeric", length = N) for(i in 1:N){ U = runif(1) X[i] = -sigma*log(1-U)} return(X)} \end{lstlisting} \item [(b)] Code: \begin{lstlisting} X = inv.cdf.gpd(1000, 12) \end{lstlisting} \item [(c)] Code: \begin{lstlisting} hist(X, probability=TRUE, col="ForestGreen", xlab="Whatever", ylab="Whatever", main="Insert canny title") curve((1/12)*exp(-x/12), from=0, to=80, col="red", lty=2, add=TRUE) \end{lstlisting} \begin{center} \includegraphics[width=5cm]{gpd.pdf} \end{center} \item [(d)] Code: \begin{lstlisting} summary(X) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00201 3.29576 7.96130 11.63675 16.07250 98.41827 \end{lstlisting} The mean and quartiles match up very closely to those obtained from the wind speed data. \newpage \item [(e)] Code: \begin{lstlisting} length(X[X>50])/length(X) [1] 0.01 \end{lstlisting} So we might expect this wind speed to be exceeded once every 100 observations -- or once every 4 days. \item [(f)] Code: \begin{lstlisting} X1 = rexp(1000, 1/12) summary(X1) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00696 3.20706 7.66012 11.89842 16.94143 76.07372 \end{lstlisting} Summaries for \texttt{X} and \texttt{X1} very similar - this form of the Generalised Pareto Distribution is equivalent to an exponential distribution with rate $1/\sigma$. \end{itemize} \end{itemize} \newpage \section*{Solutions: Practical 3} \begin{itemize} \item [1.]\begin{itemize} \item [(a)] Code: \begin{lstlisting} Y1 = sqrt(runif(1000, 1, 5)) Y2 = log(runif(1000, 1, 5)) \end{lstlisting} \item [(b)] For $Y1$: \[ F_{Y1}(y) = \pr{Y_{1}\leq y} = \pr{\sqrt{X} \leq y} = \pr{X \leq y^{2}} = F_{X}(y^{2}). \] Now \begin{eqnarray*} f_{Y1}(y) &=& \frac{d}{dy}F_{X}(y^{2}) \\ &=& f_{X}(y^{2}) \times 2y \quad \textrm{by chain rule} \\ &=& 0.25 \times 2y = 0.5y, \quad \sqrt{1} \leq y \leq \sqrt{5}, \end{eqnarray*} since $f_{X}(x) = 0.25$ for all $x$ (Uniform distribution with `base length' 4, hence `height' must be 0.25). For $Y2$: \[ F_{Y2}(y) = \pr{Y_{2} \leq y} = \pr{\log{X} \leq y} = \pr{X \leq e^{y}} = F_{X}(e^{y}). \] Now \begin{eqnarray*} f_{Y2}(y) &=& \frac{d}{dy}F_{X}(e^{y}) \\ &=& f_{X}(e^{y}) \times e^{y} \quad \textrm{by chain rule} \\ &=& \frac{e^{y}}{4}, \quad \log{1} \leq y \leq \log{5}. \end{eqnarray*} \item [(c)] Code: \begin{lstlisting} par(mfrow=c(1, 2)) hist(Y1, probability = TRUE) curve(0.5*x, from=1, to=2.3, add=TRUE, col="red", lwd=2, lty=2) hist(Y2, probability=TRUE) curve(0.25*exp(x), from=0, to=1.6, add=TRUE, col="blue", lwd=2, lty=2) \end{lstlisting} \begin{center} \includegraphics[width=7.5cm]{transformQ.pdf} \end{center} \end{itemize} \item [2.] \begin{itemize} \item [(a)] Code: \begin{lstlisting} Y = (rnorm(1000))^2 \end{lstlisting} \item [(b)] Code: \begin{lstlisting} hist(Y, probability=TRUE, main="", xlab="whatever", ylab="whatever") \end{lstlisting} \item [(c)] Code: \begin{lstlisting} curve((1/(2^(0.5)*gamma(0.5)))*x^(-0.5)*exp(-x/2), from=0, to=12, col="red", lty=2, add=TRUE) \end{lstlisting} \begin{center} \includegraphics[width=5cm]{chisq.pdf} \end{center} \item[(d)] Code: \begin{lstlisting} length(Y[Y>3.84])/length(Y) [1] 0.041 \end{lstlisting} \item[(e)] Code: \begin{lstlisting} Y1 = rchisq(1000,1) 1-pchisq(3.84,1) [1] 0.05004352 \end{lstlisting} \end{itemize} \item [3.] Code: \begin{lstlisting} Y = tan(runif(1000, -p1/2, pi/2)) hist(Y, probability=TRUE, col="misty rose") \end{lstlisting} \begin{center} \includegraphics[width=8cm]{cauchy999.pdf} \end{center} Looks very odd - with \textit{very} heavy tails! What is it?? We can use the transformation method, in exactly the same way we have in previous questions, to find that \[ f_{Y}(y) = \frac{1}{\pi(1+y^{2})}, \quad -\infty < y <\infty. \] This is the PDF of the Cauchy distribution: \begin{lstlisting} Y2 = rcauchy(1000) \end{lstlisting} \end{itemize} \newpage \section*{Solutions: Practical 4} \begin{itemize} \item [1.] Recognise that $Y|X=x\sim N(\mu=x, \sigma^2 = 1)$. Then: \begin{lstlisting} Q1 = function(N){ output = matrix(0, 2, N) x = rnorm(N, 100, 5) y = rnorm(N, x, 1) output[1, ] = x output[2, ] = y return(output)} \end{lstlisting} \begin{lstlisting} test = Q1(1000) plot(test[1,], test[2,], xlab="X", ylab="Y") cor(test[1,], test[2,]) [1] 0.980935 \end{lstlisting} \begin{center} \includegraphics[width=8cm]{ex4a.pdf} \end{center} \item [2.] \begin{enumerate}[(a)] \item We have $\mu_x = 1$; $\mu_y = 3$; $\sigma_x = \sqrt{2}$; $\sigma_y = \sqrt{1}$ \item $\rho = \textsf{cov}(X,Y)/\sigma_x \sigma_y = -0.8\sqrt{2}/\sqrt{2}\sqrt{1} = -0.8$ (or just recognize from form of variance matrix) \item From property (2) on page 25 of the lecture notes, we have \[ Y|X = x \sim N\left(3+-0.8\frac{\sqrt{1}}{\sqrt{2}}(x-1), 1^{2}(1-0.64)\right) \equiv N(3-0.5656854(x-1), 0.36). \] Then: \begin{lstlisting} Q2 = function(N){ output = matrix(0, 2, N) x = rnorm(N, 1, sqrt(2)) y = rnorm(N, 3-0.5656854*(x-1), sqrt(0.36)) output[1, ] = x output[2, ] = y return(output)} test = Q2(1000) plot(test[1,], test[2,], xlab="X", ylab="Y") cor(test[1,], test[2,]) [1] -0.8164408 \end{lstlisting} \begin{center} \includegraphics[width=8cm]{ex4b.pdf} \end{center} \item \begin{lstlisting} library(MASS) npts = 1000 mu = c(1, 3) sigma = matrix(data = c(2, -0.8*sqrt(2), -0.8*sqrt(2), 1), nrow=2, byrow=TRUE) test2 = mvrnorm(npts, mu, sigma) plot(test2[,1], test2[,2]) cor(test2[,1], test2[,2]) [1] -0.7970329 \end{lstlisting} \end{enumerate} \end{itemize} \newpage \section*{Solutions: Assignment} \subsection*{Question 1} \begin{enumerate}[(a)] \item $a=10$; $b=12$; $f(x) = 120x^9(1-x^{10})^{11}, 0\leq x \leq 1.$ \textcolor{magenta}{[No marks]} \item \textcolor{magenta}{[7 marks] -- no marks for code; [4] for plot; [2] for simulation grid; [1] for labels} \begin{lstlisting} x = seq(0, 1, 0.0001) y = 120*(x^9)*(1-(x^10))^11 plot(x, y, type="l", main="Plot of my function", xlim=c(-1,2)) abline(h=5, lty=2, col="red") abline(v=0, lty=2, col="red") abline(v=1, lty=2, col="red") \end{lstlisting} \begin{center} \includegraphics[width=6cm]{ass1.pdf} \end{center} \item \textcolor{magenta}{[18 marks] -- [2] for correct generation of points; [2] for evaluation of function; [4] for \texttt{if} statement and counter; [2] for calculation of \texttt{P}; [2] for adjustment by area of simulation grid (1.5 in this example, but look at their plot; [3] for correct overall structure of the main function; [1] for inclusion of comments; [2] for correct implementation and answer (everyone's answer should be close to 1)} \begin{lstlisting} MCintegrate = function(N){ no_of_hits = 0 for(i in 1:N){ x = runif(1, 0, 1) #Generate the x co-ordinate y = runif(1, 0, 5) #Generate the y co-ordinate f = 120*(x^9)*(1-(x^10))^11 #Evaluate the function at x if(y