\documentclass[article,nojss]{jss}
\DeclareGraphicsExtensions{.pdf,.eps}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Add-on packages and fonts
\usepackage{graphicx}
\usepackage{amsmath}

\newcommand{\noun}[1]{\textsc{#1}}
%% Bold symbol macro for standard LaTeX users
\providecommand{\boldsymbol}[1]{\mbox{\boldmath $#1$}}

%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}
\usepackage{array} % tabel commands
\setlength{\extrarowheight}{0.1cm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\newcommand{\ls}{\textbf{\textsf{limSolve }}}
\newcommand{\rs}{\textbf{\textsf{rootSolve }}}
\newcommand{\R}{\proglang{R }}
\newcommand{\ds}{\textbf{\textsf{deSolve }}}
\newcommand{\dd}{\textbf{\textsf{ddesolve }}}
\newcommand{\rb}[1]{\raisebox{1.5ex}{#1}}

\title{Using \R for scientific computing}
\Plaintitle{Using R for scientific computing}

\Keywords{Variables, Functions, Figures, Interpolation, Fitting, Roots,
  Ordinary differential equations, \proglang{R}}

\Plainkeywords{Variables, Functions, Figures, Interpolation, Fitting,
  Roots, Ordinary differential equations, R}


\author{Karline Soetaert\\
Centre for Estuarine and Marine Ecology\\
Netherlands Institute of Ecology\\
The Netherlands\\
June 2009
}

\Plainauthor{Karline Soetaert}

\Abstract{
  \R \citep{R2008} is the open-source (read: free-of-charge) version of
  the language S.
  It is best known as a package that performs statistical analysis and graphics.
  However, \R is so much more: it is a high-level language
  in which one can perform complex calculations, implement new methods, and
  make high-quality figures.

  \R has high-level functions to operate on matrices, perform numerical
  integration, advanced statistics,... which are easily triggered and which
  make it ideally suited for data-visualization, statistical analysis and
  mathematical modeling.

  It is the aim of these lecture notes to make you acquainted with the
  \R language.
  The lecture notes are based on a book \citep{Soetaert08} about ecological
  modelling in which \R is extensively used for developing, applying and
  visualizing simulation models.

  There are many excellent sources for learning the \R (or S) language.
  \R comes with several manuals that can be consulted from
  the main \R program (Help/Manuals). R-intro.pdf is a good start.
  Many other good introductions to \R are available, some freely on the web,
  and accessible via the \R web site (www.r-project.org).
  My favorite is the \R introduction by Petra Kuhnert and Bill Venables
  \citep{Kuhnert05}, but beware: this "introduction" comprises more than
  300 pages!

}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Karline Soetaert\\
  Centre for Estuarine and Marine Ecology (CEME)\\
  Netherlands Institute of Ecology (NIOO)\\
  4401 NT Yerseke, Netherlands
  E-mail: \email{k.soetaert@nioo.knaw.nl}\\
  URL: \url{http://www.nioo.knaw.nl/users/ksoetaert}\\
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands.
%% need no \usepackage{Sweave}
%\VignetteIndexEntry{Using R for Scientific Computing - Lecture notes: questions}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document
\begin{document}
\SweaveOpts{engine=R,eps=FALSE}
\SweaveOpts{keep.source=TRUE}
<<preliminaries,echo=FALSE,results=hide>>=
library(marelac)
library(marelacTeaching)
options(prompt = ">")
options(width=70)
@

\maketitle

\section{introduction}
\subsection{The R-software}

\subsubsection{Installing R}
  R is downloadable from the following web site: http://www.r-project.org/
  Choose the precompiled binary distribution.
  On this website, you will also find useful documentation.
  To use \R for the examples in this course, several packages need to be
  downloaded.
  \begin{itemize}
    \item	\pkg{deSolve}. Performs integration. \citep{deSolve}
    \item	\pkg{rootSolve}. Finds the root of equations \citep{rootSolve}.
    \item	\pkg{scatterplot3d}. For 3-D graphics. \citep{scatterplot3D}
    \item	\pkg{seacarb}. Aquatic chemistry. \citep{seacarb}
    \item	\pkg{marelac}. Functions and constants from the marine
      and lacustrine sciences \citep{marelac}.
    \item	\pkg{marelacTeaching}. Data sets, used in this tutorial.
      \citep{marelacTeaching}.
  \end{itemize}

  If you run \R within windows, downloading specific packages can best be
  done within the \R program itself.
  Select menu item \code{packages / install packages},
  choose a nearby site (e.g. France (Paris)) and select the package you
  need. If you install package \pkg{marelacTeaching}
  then all other packages will be automatically installed as well.
\subsubsection	{Other useful software}
  I  prefer to run \R from within the Tinn-R editor, which can be
  downloaded from URL http://sourceforge.net/projects/tinn-r and
  http://www.sciviews.org/Tinn-R. This editor provides R-sensitive syntax
  and help. Download the latest Tinn-R setup file and install it.
  From within the Tinn-R program, you launch \R via the menu (\code{R/start
  preferred Rgui}).
 
\subsection{Quick overview of R}
  R-code is highly readable, once you realise that:
  \begin{itemize}
    \item	\code{<-} is the assignment operator.
    \item	everything starting with \code{\#} is considered a comment.
    \item	R is case-sensitive: \code{a} and \code{A} are two different objects.
  \end{itemize}
\subsection{Console versus scripts}
  There are two ways in which to work with R.
  \begin{enumerate}
    \item We can type commands into the \R console window at the command
      prompt (>) and use \R as a powerful scientific calculator.

      For instance, enter in the console window:
<<>>=
 pi*0.795^2 ; 25*6/sqrt(67) ; log(25)
@
     Here \code{sqrt} and \code{log} are built-in functions in R; \code{pi}
     is a built-in constant; the semi-colon (;) is used to separate R-commands.

     In the console window, the <UP> and <DOWN> arrow keys are used to
     navigate through previously typed sentences.


    \item Alternatively, we can create R-scripts in an editor (e.g. Tinn-R)
      and save them in a file ("filename.R") for later re-use.
      R-scripts are sequences of R-commands and expressions. These scripts
      should be submitted to \R before they are executed.
      This can be done in several ways:
      \begin{itemize}
        \item	by typing, in the R-console window:
\begin{verbatim}
> source ("filename.R") 
\end{verbatim}
         \item	by opening the file, copying the R-script to the clipboard
           (ctrl-C) and pasting it (ctrl-V) into the R-console window
         \item	If you do not use the tinn-R editor, the file is opened
           as an R-script from within the \R console.
           After selecting the script, and pressing the "send" button
           the statements are executed and the cursor moved to the next line.
         \item	If you do use the Tinn-R editor,
           \includegraphics[width=0.4\textwidth]{tinn1}
           you can either submit the entire file (buttons 1,2), selected
           parts of the text (buttons 3,4), submit marked blocks (buttons 5,6)
           or line-by-line (last button).
      \end{itemize}
  \end{enumerate}
  Throughout these notes, the following convention is used:
\begin{verbatim}
> 3/2
\end{verbatim}
  denotes input to the console window (\code{>} is the prompt)
\begin{verbatim}
[1] 1
\end{verbatim}
  is \R output, as written in the console window
\begin{verbatim}
 getwd()
\end{verbatim}
  is an \R statement in a script file (it lacks the prompt).

  A screen capture of a typical Tinn-R session, with the Tinn-R editor
  (upper window) and the R-console (lower window) is given below.
  A script file is opened in the Tinn-R editor. Note the context-sensitive
  syntax (green=comments, blue= reserved words, rose = R-parameters).
  Several lines of R-code have been selected (blue area) and sent to the
  R-console, which has produced the graphics window that floats
  independently from the other windows.

\includegraphics[width=0.8\textwidth]{capture}


\subsubsection{Getting help, examples, demonstrations}
  R has an extensive help facility. Apart from the Help window launched
  from the Help menu, or the HTML help facility, it is also available from
  the command line prompt.
  For instance, typing
\begin{verbatim}

> ?log
> ?sin
> ?sqrt
> ?round
> ?Special
\end{verbatim}

  will explain about logarithms and exponential functions, trigonometric
  functions, and other functions.
\begin{verbatim}
> ?Arithmetic
\end{verbatim}
  lists the arithmetic operations in R.
\begin{verbatim}
> help.search("factor")
\end{verbatim}
  will list occurrences of the word <factor> in R-commands.

  Sometimes the best help is provided by the very active mailing list.
  If you have a specific problem, just type \code{R: <problem>}
  on your search engine. Chances are that someone encountered the problem
  and it was already solved.

  Most of the help files also include examples. You can run all of them by
  using R-statement \code{example}.

  For instance, typing into the console window:
\begin{verbatim}
> example(matrix)
\end{verbatim}
  will run all the examples from the matrix help file.

\begin{verbatim}
> example(pairs)
\end{verbatim}
  will run all the examples from the pairs help file.
  (! try this ! \code{pairs} is a very powerful way of visualizing
  pair-wise relationships).

  Alternatively, you may select one example, copy it to the clipboard
  (ctrl-C for  windows users) and then paste it (ctrl-V) in the console window.
  In addition, the \R main software and many R-packages come with
  demonstration material. Typing
\begin{verbatim}
> demo()
\end{verbatim}
  will give a list of available demonstrations in the main software.
\begin{verbatim}
> demo(graphics)
\end{verbatim}
  will demonstrate some simple graphical capabilities.
\subsubsection{	Small things to remember}
  \begin{itemize}

    \item 	Pathnames in \R are written with forward slashes ("/") ,
      although in windows, backslashes, ($\backslash$), are used.
      Thus, to set a working directory in R:
\begin{verbatim}
setwd("C:/R code/") 
\end{verbatim}
    \item	If a sentence on one line is syntaxically correct, \R will
      execute it, even if it is the intention that it proceeds on the next
      line. For instance if we write:
<<>>=
A <- 3 + cos(pi)
    - sqrt(5)
@
      then A will get the value ($3+cos(\pi)$) and \R will print the value
      of ($–\sqrt(5)$).
      In contrast, in the next lines:
\begin{verbatim}
> 3 + cos(pi) -
>     sqrt(5)
[1] -0.236068
\end{verbatim}

      \R will print the value of $3+cos(\pi)-\sqrt(5)$: as the sentence
      on the first line was not syntaxically finished,
      \R has (correctly) assumed that it continued on the next line.

      Be careful if you want to split a complex statement over several lines !
      These errors are very difficult to trace, so it is best to avoid them.
  \end{itemize}
\subsubsection{Exercises. Using R as a calculator}
  It is very convenient to use \R as a powerful calculator. This can best
  be done from within the R-console.

  \begin{enumerate}
    \item Use the console to calculate the value of:
    \begin{itemize}
      \item	$(4/6*8-1)^{2/3}$
      \item	ln(20)
      \item	$log_2(4096)$
      \item	$2*\pi*3$
      \item $\sqrt {2.3^2  + 5.4^2  - 2*2.3*5.4*\cos (\pi /8)}$
    \end{itemize}
    tip: you may need to look at the help files for some of these
    functions. typing \code{?"+"} will open a help file with the common
    arithmetic operators.

    \item Now write the R-statements in a script file, using the Tinn-R
    editor. Try the various ways in which to submit the statements to \R.
  \end{enumerate}
\clearpage

\section{R-variables}
  \R calculates as easily with vectors, matrices and arrays as with single
  numbers.

  \R also includes more complex structures such as data frames and lists,
  which allow to combine several types of data.

  Learning how to create these variables, how to address them and modify
  them is essential if you want the make good use of the \R software.

\subsection{Numbers, vectors, matrices and arrays}
\subsubsection{Assignment}
  When variables are used, they need to be initialised with numbers.
<<>>=
A <- 1
B <- 2
A+B
@
  R can take as arguments for its functions single numbers, vectors,
  matrices, or arrays.

<<>>=
V<-factorial(10)
@
  Calculates $10! (=1\cdot 2\cdot 3\cdot 4\cdot 5\cdot \ldots \cdot 10)$.
  The operator \code{<-} assigns the result of this calculation to variable V.
  V can then be used in subsequent calculations:
<<>>=
V/10
@
  Note that the assignment of a value to V does not display it on the window.
  To display V we simply write:
<<>>=
 V
@
  Alternatively, we may assign the result of calculations to a variable AND
  view the results, by embracing the statement between brackets:
<<>>=
 (X <- sin(3/2*pi) )
@
  Apart from integers, real and complex numbers, \R also recognizes infinity
  (Inf) and Not a Number (NaN). Try:
\begin{verbatim}
> 1/0
> 0/0
> 1e-8 * 1000
\end{verbatim}
  (where the \code{e-8} notation denotes $10^{-8}$).

\subsubsection{Vectors}
  Vectors can be created in many ways:
  \begin{itemize}
    \item	Using R-function \code{vector}
    \item	The function \code{c()} combines numbers into a vector
      \footnote{This is perhaps THE most important function in R}
    \item	The operator "\code{:}" creates a sequence of values, each
      1 larger (or smaller) than the previous one
    \item	A more general sequence can be generated by R-function \code{seq}
    \item	The same quantity is repeated using R-function \code{rep}
  \end{itemize}
  For instance, the commands:
<<>>=
c(0, pi/2, pi, 3*pi/2, 2*pi)
seq(from=0,to=2*pi, by=pi/2 )
seq(0, 2*pi, pi/2 )
@
  will all create a vector, consisting of: 0, $\pi,\ldots 2*\pi$.

  Note that R-function \code{seq} takes as input (amongst others)
  parameters \code{from}, \code{to} and \code{by} ($2^{nd}$ example).
  If the order is kept, they not be specified by name ($3^{rd}$ example).

  The next command calculates the sine of this vector and outputs the result:
<<>>=
sin( seq(0, 2*pi, pi/2 ))
@
  Function \code{rep} is used to repeat elements:
<<>>=
rep(1,times=5)
rep(c(1,2),times=5)
c(rep(1,5),rep(2,5))
@
  The next statements:
<<>>=
 V <- 1:20
 sqrt(V)
@
  create a sequence of integers between 1 and 20 and take the square root
  of all of them, displaying the result to the screen.
  The operator \code{<-} assigns the sequence to V.

  Some other examples of the \code{:} operator are:
<<>>=
(V <- 0.5:10.5)
6:1
@
  Finally, the statements:
<<>>=
V <- vector(length=10)
FF<- vector()
@
  generate a vector V comprising 10 elements, and a vector FF of unknown length.

  Note: a peculiar feature of \R is that the elements of a vector can also
  be given \code{names}:
<<>>=
(fruit <- c(banana=1, apple=2, orange =3))
names(fruit)
@
\subsubsection{Matrices}
  Matrices can also be created in several ways:
  \begin{itemize}
    \item	By means of R-function \code{matrix}
    \item	By means of R-function \code{diag} which constructs a diagonal matrix
    \item	The functions \code{cbind} and \code{rbind} add columns and
      rows to an existing matrix, or to another vector
  \end{itemize}
  The statement:
<<>>=
A <-matrix(nrow=2,data=c(1,2,3,4))
@
  creates a matrix A, with two rows, and, as there are four elements,
  two columns.
  Note that the data are inputted as a vector (using the \code{c()} function).

  The next two statements display the matrix followed by the square root
  of its elements:
<<>>=
A
sqrt(A)
@
  By default, \R fills a matrix column-wise (see the example above). However,
  this can easily be overruled, using parameter \code{byrow}:
<<>>=
(M <-matrix(nrow=4, ncol=3, byrow=TRUE, data=1:12))
@
  The unity matrix (I) is created using R-function diag:
<<>>=
diag(1,nrow=2)
@

  The \emph{names} of columns and rows are set as follows:
<<>>=
rownames(A) <- c("x","y")
colnames(A) <- c("c","b")
A
@
  note that we also use the \code{c()} function here ! Row names and column
  names are in fact vectors containing strings.

  Matrices can also be created by combining (binding) vectors, e.g. rowwise:
<<>>=
V <- 0.5:5.5
rbind(V,sqrt(V))
@
  t(A) will transpose matrix A (interchange rows and columns).
<<>>=
t(A)
@
\subsubsection{Arrays}
  Arrays are multidimensional generalizations of matrices; matrices and
  arrays in \R are actually vectors with a dimension attribute.
  A multi-dimensional array is created as follows:
<<>>=
AR <-array(dim=c(2,3,2),data=1)
@
  In this case AR is a 2*3*2 array, and its elements are all 1.
\subsection{Dimensions}
  The commands
\begin{verbatim}
> length(V)
> dim(A)
> ncol(M)
> nrow(M)
\end{verbatim}
  Will return the length (total number of elements) of (vector or matrix) V,
  the dimension of matrix or array A, and the number of columns and rows of
   matrix M respectively.

\subsection{Selecting and extracting elements}
  To select subsets of vectors or matrices, we can either
  \begin{itemize}
    \item	specify the numbers of the elements that we want (simple indexing)
    \item	specify a vector of logical values (TRUE/FALSE) to indicate which
    elements to include (TRUE) and which not to include (FALSE). This uses
    logical expressions
  \end{itemize}
\subsubsection{Simple indexing}
  The elements of vectors, matrices and arrays are indexed using the
  \code{[]} operator:
\begin{verbatim}
M[1 , 1]
M[1 , 1:2]
M[1:3 , c(2,4)]
\end{verbatim}
  Takes the element on the first row, first column of a matrix M
  ($1^{st}$ line),
  then selects the entries in the first row and first two columns
  ($2^{nd}$ line) and then the elements on the first three rows, and
  $2^{nd}$ and $4^{th}$ column of matrix M ($3^{rd}$ line).

  If an index is omitted, then all the rows ($1^{st}$ index omitted) or
  columns ($2^{nd}$ index omitted) are selected. In the following:
\begin{verbatim}
M[   ,2] <-0
M[1:3, ] <- M[1:3, ] * 2
\end{verbatim}
  first all the elements on the $2^{nd}$ column ($1^{st}$ line) of M are
  zeroed and then the elements on the first three rows of M multiplied with
  2 ($2^{nd}$ line).
  Similar selection methods apply to vectors:
\begin{verbatim}
V[1:10]
V[seq(from=1,to=5,by=2)]
\end{verbatim}
  The statement on the $1^{st}$ line takes the first 10 elements of vector V,
  whilst on the $2^{nd}$ line, the $1^{st}$, $3^{rd}$ and $5^{th}$ element of
  vector V are selected.
\subsubsection{Logical expressions}
  Logical expressions are often used to select elements from vectors and
  matrices that obey certain criteria.

  R distinguishes logical variables TRUE and FALSE, represented by the
  integers 1 and 0.
\begin{verbatim}
> ?Comparison
> ?Logic
\end{verbatim}
  will list the relational and logic operators available in R.
  The following will return TRUE for values of sequence V that are positive:
<<>>=
(V <- seq(-2,2,0.5))
 V>0
@
  while
<<>>=
V [V > 0]
@
  will select the positive values from V,
<<>>=
V [V > 0] <- 0
@
  will zero all positive elements in V,
<<>>=
 sum(V < 0)
@
  will return the number of negative elements: it sums the TRUE (=1) values, and
<<>>=
V [V != 0]
@
  will display all nonzero elements from V ( "!" is the "not" operator).
  Logical tests can also be combined, using $\mid$  (the "or" operator),
  and  $\&$ ("and").
<<>>=
V [V<(-1) | V>1]
@
  will display all values from V that are < -1 and > 1. Note that we have
  enclosed "-1" between brackets (can you see why this is necessary?)
  Finally,
<<>>=
which (V == 0)
which.min (V)
@
  will return the element index of the 0-value, and of the minimum.
\subsection{Removing elements}
  When the index is preceded by a "-", the element is removed.
\begin{verbatim}
> M[ ,-1]
\end{verbatim}
  Will show the contents of matrix M, except the first column.

\begin{verbatim}
> x <- x[-1]
> M <- M[-1, ]
> V <- V[-V>=0]
\end{verbatim}
  will delete the $1^{st}$ element of x, ($1^{st}$ line),  the $1^{st}$
  row of M ($2^{nd}$ line), and the positive elements of V ($3^{rd}$ line).
  For more information, type
\begin{verbatim}
> ?Extract
\end{verbatim}
\subsection{More complex data structures}
  \R also allows creating more complex structures such as data frames and lists.
\subsubsection{lists}
  A list is a combination of several objects; each object can be of
  different length:
\begin{verbatim}
> list(Array = AR, Matrix = M)
\end{verbatim}
  will combine the previously defined array AR and matrix M.
\subsubsection{data.frames}
  These are combinations of different data types (e.g. characters,
  integers, logicals, reals), arranged in tabular format:
<<>>=
genus <- c("Sabatieria","Molgolaimus")
dens  <- c(1,2)
Nematode <-data.frame(genus=genus,density=dens)
Nematode
@
  In the example above, the data.frame \code{Nematode} contains two columns,
  one with strings (the genus name), one with values (the densities).
  Data.frames are in fact special cases of lists, consisting of vectors
  with equal length.
  Many matrix-operations work on data.frames with a single data type, but
  there exist also special operations on data.frames.
\subsubsection{Selecting data from data.frames and lists}
  Data.frames and lists can be accessed by their names, or by the \code{[ ]}
  or \code{[[ ]]} operator.
  The object resulting from a selection using single brackets \code{[ ]},
  will be a data.frame respectively a list itself;
  with double brackets \code{[[ ]]}, one obtains a vector (data.frames)
  or a variable data-type (lists).

  For instance:
<<>>=
Nematode$density/sum(Nematode$density)
@
  will divide all density values (1,2) by the summed density.
<<>>=
mean(Nematode[,2])
@
  will calculate the mean of Nematode density (the $2^{nd}$ column).

  Try:
\begin{verbatim}
> Nematode$genus
> Nematode[1]
> Nematode[[1]]
\end{verbatim}

  These statements will all output the two genus names, but in a
  different format.
\begin{verbatim}
> ?Extract
\end{verbatim}
  also explains the various ways in which to extract elements from lists
  and data frames.
\subsection{Data Conversion}
  Conversion from one type of data structures to another can easily be done,
  e.g. by:
\begin{verbatim}
> as.data.frame(M)
> as.vector(A)
\end{verbatim}
  If unsure about the type, you can write:
\begin{verbatim}
> is.data.frame(M)
> is.vector(A)
\end{verbatim}
  Or you can display the data type by:
\begin{verbatim}
> class(M)
\end{verbatim}
\subsection{Data import from external sources}
  In all previous examples, data was entered from the console (or from a
  script file).
  There are ample facilities to import data from external sources. Most
  often, we will use functions \code{read.table}, \code{read.csv}, or
  \code{read.delim} to read matrices or data frames written in tabular
  form as text files.
  R also has plenty of built-in data sets. They are listed by:
\begin{verbatim}
> data()
\end{verbatim}

\subsection{Exercises}
  Creating and manipulating matrices and vectors is essential if we want to
  use \R as a mathematical tool.
  Although this has been implemented in a consistent way in \R, it is
  not simple for novice users!
  Practice is the best teacher, so you will get plenty of exercise.

  Most of the exercises can be answered with one single \R-statement.
  However, as these statement smay be quite complicated, it is often simpler
  to first break them up into smaller parts, after which they are
  merged into one.
\subsubsection{Vectors, sequences.}
  \begin{itemize}
    \item	Use R-function \code{mean} to estimate the mean of two numbers,
      9 and 17. (you may notice that this is not as simple as you might think!).
    \item Vector V
    \begin{itemize}
      \item	Create a vector, called V, with even numbers, between 16 and 56.
        Do not use loops. (tip: use R-function  \code{seq} )
      \item	Display this vector
      \item	What is the sum of all elements of V? Do not use loops; there
        exists an R-function that does this; the name of this function
        is trivial.
      \item	Display the first 4 elements of V
      \item	Calculate the product of the first 4 elements of V
      \item	Display the $4^{th}$, $9^{th}$ and $11^{th}$ element of V .
        (tip: use the \code{c()} function).
    \end{itemize}
    \item Vector W
    \begin{itemize}
      \item	Create a new vector, W, which equals vector V, multiplied with 3;
        display its content.
      \item	How many elements of W are smaller than 100?

        First create a new vector that contains only the elements from W<100
          (call it W100), then calculate the length of this new vector.
      \item	Now perform the same calculation, in one \R statement.
    \end{itemize}
    \item	Create a sequence that contains the values (1,$1/2$,$1/3$,$1/4$,
      $\ldots$, $1/10$)
    \item	Compute the square root of each element
    \item	Compute the square ($^2$) of each element
    \item	Create a sequence with values ($0/1$,$1/2$,$2/3$,$3/4$,
      $\ldots$, $9/10$)
    \item Vector U
    \begin{itemize}
      \item	Create a vector, U, with 100 random numbers, uniformly
        distributed between -1 and 1.

      tip: R-statement \code{runif} generates uniformly distributed
        random numbers; use \code{?runif} to see how it works.
      \item	Check the range of U; all values should be within -1 and +1.

      tip: there exists an R-function to do that-its name is trivial.
      \item	Calculate the sum and the product of the elements of U
      \item	How many elements of U are positive?
      \item	Zero all negative values of U.
      \item	Sort U
    \end{itemize}
    \item Vectors x, y
    \begin{itemize}
      \item	Create two vectors: vector x, with the elements: 2,9,0,2,7,4,0
        and vector y with the elements 3,5,0,2,5,4,6 (in that order). (tip:
        use the c() function).
      \item	Divide all the elements of y by the elements of x.
      \item	Type in the following commands; try to understand:
      \begin{itemize}
        \item	x>y
        \item	x==0
      \end{itemize}
      \item	Select all values of y that are larger than the corresponding
        values of x
      \item	Select all values of y for which the corresponding values of
        x are 0.
      \item	Remove all values of y for which the corresponding values of
        x equal 0.
      \item	Zero all elements of x that are larger or equal than 7. Show x.
    \end{itemize}
  \end{itemize}

\subsubsection{Matrices}
  \begin{itemize}
    \item	Use \R-function \code{matrix} to create a matrix with the
      following contents:
      \[
        \left[ {\begin{array}{*{20}c}
           3 & 9  \\
           7 & 4  \\
        \end{array}} \right]
      \]

    \item	display it to the screen
    \item	Use \R-function \code{matrix} to create a matrix called "A":
      \[
        \left[ {\begin{array}{*{20}c}
         1 & {1/2} & {1/3}  \\
         {1/4} & {1/5} & {1/6}  \\
         {1/7} & {1/8} & {1/9}  \\
        \end{array}} \right]
      \]

      \begin{itemize}
        \item	Take the transpose of A.
        \item	Create a new matrix, B, by extracting the first two rows
          and first two columns of A. Display it to the screen.
      \end{itemize}
    \item	Use \code{diag} to create the following matrix, called "D":
      \[
      \left[ {\begin{array}{*{20}c}
         1 & 0 & 0  \\
         0 & 2 & 0  \\
         0 & 0 & 3  \\
      \end{array}} \right]
      \]

    \item	Use \code{cbind} and \code{rbind} to augment this matrix,
      such that you obtain:
      \[
        \left[ {\begin{array}{*{20}c}
         1 & 0 & 0 & 4  \\
         0 & 2 & 0 & 4  \\
         0 & 0 & 3 & 4  \\
         5 & 5 & 5 & 5  \\
        \end{array}} \right]
      \]

      It is simplest to do this in two statements (but it can be done in one!)
    \item	Remove the second row and second column of the previous matrix
  \end{itemize}

\subsubsection{Diversity of deep-sea nematodes}
  We will now work on a data set consisting of nematode species densities,
  found in Mediterranean deep-sea sediments, at depths ranging from 160 m
  to 1220 m. The densities are expressed in number of individuals per 10 cm2.

  Nematodes are small (<1mm) worms, and they are generally very abundant in
  all marine sediments.

  The data (from \citep{Soetaert91}, originally present in an ACCESS database,
  have been exported as a table in a comma-separated-values (so-called csv)
  format. This format has the advantage that it can easily be read by text
  editors (such as the TINN-R editor) as well as in spreadsheet programs.

  Open the file "nemaspec.csv" in Tinn-R\footnote{: if you do not have this
  file, it can be found in the \pkg{marelacTeaching} package directory.
  Within R, type:

  browseURL(paste(system.file(package="marelacTeaching"), "/lecture/", sep=""))}.

  Check its structure.

  You may also open the file in EXCEL, but do not forget to close it
  before proceeding.
  EXCEL is very territorial, and will not allow another program, such as
  \R, to access a file that is open in EXCEL.

  On the first line is the heading (the names of the stations), the first
  column contains the species names.
  Before importing the file in \R, check the working directory:
\begin{verbatim}
> getwd()
\end{verbatim}
  If the file called "nemaspec.csv" is not in this directory, you may
  need to change the working directory:
\begin{verbatim}
> setwd("directory name")
\end{verbatim}
  (do not forget that \R requires "/" where windows uses " $\backslash$") .

  Make a script file in which you write the next steps; submit each line
  to \R to check its correctness.
  Read the comma-delimited file, using R-command \code{read.csv}. Type
  \code{?read.csv} if you need help.

  Specify that the first row is the heading (\code{header=TRUE}) and the
  first column contains the rownames (\code{row.names=1}).

  Put the data in data.frame \code{Nemaspec}.
\begin{verbatim}
Nemaspec <- read.csv("nemaspec.csv", header=TRUE, row.names=1)
\end{verbatim}
  Check the contents of \code{Nemaspec}. As the dataset is quite substantial,
  it is best to output only the first part of the data:
\begin{verbatim}
head(Nemaspec)
\end{verbatim}

  The rest is up to you:
  \begin{itemize}
    \item	Select the data from station M160b (the $2^{nd}$ column of Nemaspec);
      put these data in a vector called \code{dens}.

      (remember: to select a complete column, you select all rows by leaving
      the first index blanc).
    \item	Remove from vector dens, the densities that are 0. Display this
      vector on the screen.  		(Answer: [1]  6.580261  5.919719   etc$\ldots$)
    \item	Calculate N, the total nematode density of this station. The
      total density is simply the sum of all species densities (i.e. the sum
      of values in vector dens). What is the value of N ? 				(Answer :699).
    \item	Divide the values in vector dens by the total nematode density N.
      Put the results in vector p, which now contains the relative
      proportions for all species. The sum of all values in p should now
      equal 1. Check that.
    \item	Calculate S, the number of species: this is simply the length of p;
      call this value S. 				    	(Answer: S=126)
    \item	Estimate the values of diversity indices N1 and N2 and Ni, given
      by the following formulae:
      \[
      \begin{array}{l}
       N1 = e^{\sum { - p_i  \cdot \log _e (p_i )} }  \\
       N2 = 1/\left( {\sum {p_i ^2 } } \right) \\
       Ni = 1/\max (p_i ) \\
       \end{array}
      \]
      You can calculate each of these values using only one \R statement !
      (A: 90.15358, 66.77841, 22.56157)
    \item	The 126 nematode species per 10 cm2 were obtained by looking at
      all 699 individuals. Of course, the fewer individuals are determined
      to species, the fewer species will be encountered. Some researchers
      determine 100 individuals, other 200 individuals. To standardize
      their results, the expected number of species in a sample can be
      recalculated based on a common number of individuals.
      The expected number of species in a sample with size n, drawn from a
      population which size N, which has S species is given by:
      \[
        ES(n) = \sum\limits_{i = 1}^S {\left[ {1 - \frac{{\left( {_n^{N - Ni} } \right)}}{{\left( {_n^N } \right)}}} \right]}
      \]
      where Ni is the number of individuals in the ith species in the full
      sample and   is the so-called "binomial coefficient", the number of
      different sets with size n that can be chosen from a set with total
      size N.

      In R, binomial coefficients are estimated with statement \code{choose(N,n)}.

      What is the expected number of species per 100 individuals ?
      (n=100,N=699).				(A: ES(100) = 60.68971).
    \item	Print all diversity indices to the screen, which should look like:
\begin{verbatim}
        N        N0        N1        N2        Ni       ESS
699.00000 126.00000  90.15358  66.77841  22.56157  60.68971
\end{verbatim}
\end{itemize}
\clearpage

\section{R functions}
  One of the strengths of \R is that one can make user-defined functions
  that add to \R-s built-in functions.

\subsection{Function definition}
  Typically, complex functions are written in \R-script files, which are
  then submitted to \R. For instance,
<<echo=FALSE,results=hide>>=
options(prompt = " ")
@

<<>>=
Circlesurface <- function (radius) pi*radius^2
@
  defines a function (called \code{Circlesurface})  which takes as input
  argument a variable called \code{radius}  and which returns the value
  $\pi*radius^2$ (which is the surface of a circle).

  After submitting this function to \R, we can use it to calculate the
  surfaces of circles with given radius:
<<echo=FALSE,results=hide>>=
options(prompt = ">")
@

<<>>=
Circlesurface(10)
Circlesurface(1:20)
@
  the latter statement will calculate the surface of circles with
  radiuses 1, 2, ... ,20.

  More complicated functions may return more than one element:
\begin{verbatim}
Sphere <- function(radius)
{
 volume  <- 4/3*pi*radius^3
 surface <- 4 *pi*radius^2
 return(list(volume=volume,surface=surface))
}
\end{verbatim}
<<echo=FALSE,results=hide>>=
Sphere <- function(radius)
{
 volume  <- 4/3*pi*radius^3
 surface <- 4 *pi*radius^2
 return(list(volume=volume,surface=surface))
}
@
  Here we recognize
  \begin{itemize}
    \item	the function heading ($1^{st}$ line), specifying the name of
      the function (Sphere) and the input parameter (radius)
    \item	the function specification. As the function comprises multiple
      statements, the function specification is embraced by curly braces
      $\{ \ldots \}$.
    \item	The return values (last line). \code{Sphere} will return the
      volume and surface of a sphere, as a list.
  \end{itemize}
  The earth has approximate radius 6371 km, so its volume (km3) and
  surface (km2) are:
<<>>=
Sphere(6371)
@
  The next statement will only display the volume of spheres with
  radius 1, 2, $\ldots$ 5
<<>>=
Sphere(1:5)$volume
@
  Sometimes it is convenient to provide default values for the input
  parameters.

  For instance, the next function estimates the density of "standard mean
  ocean water" (in $kg~m^{-3}$),
  as a function of temperature, T, (and for salinity=0, pressure = 1 atm)
  \citep{Millero81};
  the input parameter T is by default equal to 20 $^{\circ}$C:
\begin{verbatim}
Rho_W <- function(T=20)
{
999.842594 + 0.06793952 * T - 0.00909529 * T^2 +
0.0001001685 * T^3 - 1.120083e-06 * T^4 + 6.536332e-09 * T^5
}
\end{verbatim}
<<echo=FALSE,results=hide>>=
Rho_W <- function(T=20)
{
999.842594 + 0.06793952 * T - 0.00909529 * T^2 +
0.0001001685 * T^3 - 1.120083e-06 * T^4 + 6.536332e-09 * T^5
}
@
  By ending the first sentence with a \code{+} we made clear that the
  statement is not finished and continues on the next line. It would have
  been wrong to put the \code{+} on the next line.

  Calling the function without specifying temperature, uses the default value:
<<>>=
Rho_W()
Rho_W(20)
Rho_W(0)
@
\subsection{Programming}
  \R has all the features of a high-level programming language:
\subsubsection{If, else, ifelse constructs}
  Try to understand the following:
\begin{verbatim}
Dummy <- function (x)
{
if ( x<0  ) string <- "x<0"     else
if ( x<2  ) string <- "0>=x<2"  else
            string <- "x>=2"
print(string)
}
\end{verbatim}
<<echo=FALSE,results=hide>>=
Dummy <- function (x)
{
if ( x<0  ) string <- "x<0"     else
if ( x<2  ) string <- "0>=x<2"  else
            string <- "x>=2"
print(string)
}
@
<<>>=
Dummy(-1)
Dummy(1)
Dummy(2)
@
  Note that we have specified the \code{else} clause on the same line as
  the \code{if} part so that \R knows that the statement is continued on
  the next line!

  If and else constructs involving only one statement can be combined:
<<>>=
x<-2
ifelse (x>0, "positive", "negative,0")
@
\subsubsection{Loops}
  Loops allow a set of statements to be executed multiple times:

  The \code{for} loop iterates over a specified set of values. In the
  example below, the variable "i" takes on the values (1,2,3):
<<>>=
for (i in 1:3) print(c(i,2*i,3*i))
@
  \code{while} and \code{repeat} will execute until a specified
  condition is met.
<<>>=
i<-1 ; while(i<3) {print(i); i<-i+1}
@
  \code{break} exits the loop

  \code{next} stops the current iteration and advances to the next iteration.
<<>>=
i<-1
repeat
{
  print(i)
  i <-i+1
  if(i>2) break
}
@
  The curly braces $\{ \ldots \}$ embrace multiple statements that are
  executed in each iteration.

  Note: loops are implemented very inefficiently in \R and should be
  avoided as often as possible.
  Fortunately, \R offers many high-level commands that operate on vectors
  and matrices. These should be used as much as possible!

  For more information about if constructs and loops, type
\begin{verbatim}
> ?Control
\end{verbatim}
\subsection{R- packages}
  A package in \R is a file containing many functions that perform certain
  related tasks.
  Packages can be downloaded from the \R website.

  Once installed, we generate a list of all available packages, we load
  a package and we obtain a list with its contents by the following commands:
\begin{verbatim}
>library()
>library(deSolve)
>library(help=deSolve)
>help(package=deSolve)
\end{verbatim}

\subsection{Exercises}
\subsubsection{R-function sphere}
  Extend the \code{Sphere} function with the circumference of the sphere
  at the place of maximal radius.
  The formula for estimating the circumference of a circle with radius r
  is: $2\cdot \pi \cdot r$.
  What is the circumference of the earth near the equator?

\subsubsection{An R-function to estimate saturated oxygen concentrations}
  The saturated oxygen concentration in water (µmol kg-1), as function of
  temperature (T), and salinity (S) can be calculated by:
  $SatOx = e^A $      where :

  A= -173.9894 + 25559.07/T + 146.4813* loge(T/100) -22.204*T/100 + S *

  (-0.037362+0.016504*T/100-0.0020564 *T/100*T/100)

  and T is temperature in Kelvin (Tkelvin = Tcelsius+273.15).

  Tasks:
  \begin{itemize}
    \item	Make a function that implements this formula; the default
      values for temperature and salinity are 20$^{\circ}$C and 35
      respectively.
    \item	What is the saturated oxygen concentration at the default
      conditions? (A: 225.2346)
    \item	Estimate the saturated oxygen concentration for a range of
      temperatures from 0 to 30°C, and salinity 35. (Tip: no need to use loops).
  \end{itemize}

\subsubsection{Loops}
  The Fibonacci numbers are calculated by the following relation:
  $F_n=F_{n-1}+F_{n-2}$

  with $F_1= F_2 =1$

  Tasks:
  \begin{itemize}
    \item	Compute the first 50 Fibonacci numbers; store the results in
      a vector (use R-command \code{vector} to create it). You have to
       use a loop here
    \item	For large n, the ratio $F_n/F_{n-1}$ approaches the
      "golden mean": $(1 + \sqrt 5 )/2$
    \item	What is the value of $F_{50}/F_{49}$; is it equal to the
      golden mean?
    \item	When is n large enough? (i.e. sufficiently close (<$1e^{-6}$)
      to the golden mean)
  \end{itemize}

\subsubsection{Diversity of deep-sea nematodes for all stations}
  \begin{itemize}
    \item	Starting from your code to estimate diversity indices for
      deepsea station M160b, now write a loop that does so for all
        the stations in Nemaspec.
    \item	First create a matrix called div, with the number of rows
      equal to the number of deepsea stations, and with 6 columns, one
        for each diversity index. This matrix will contain the diversity values.
    \item	The column names of div are: "N", "N0", "N1", "N2", "Ninf", "ESS"

      The row names of matrix div are the station names (= the column names
      of Nemaspec). Tip: Use R-command \code{colnames()}, \code{rownames()}
    \item	Now loop over all columns of data frame Nemaspec, estimate the
      diversity indices and put the results in the correct row of matrix div:
\begin{verbatim}
for (i in 1:ncol(Nemaspec))
 {
   # you have to write this part of the code
 }
\end{verbatim}
    \item	Show matrix div to the screen
  \end{itemize}

\subsubsection{Diversity indices as a function}
  Based on the results obtained in previous section, make a function that
  will calculate the diversity indices for any data matrix.

\subsubsection	{Rarefaction diversity}
  \footnote{This question requires significant thought and imagination;
  there are several ways to do this.}
  If you still have time and the courage: try an alternative way of
  estimating the number of species per 100 individuals
  by taking random "subsamples" of 100 individuals and estimating the
  number of species from this subsample.

  If the procedure is repeated often enough, the mean value should converge
  to the expected number of species, ESS(100);
  this is the rarefaction method of \citep{Sanders68}.

  You may need the following \R-functions:
  \begin{itemize}
    \item	\code{round} (converting reals to integers),
    \item	\code{cumsum} (take a cumulative sum),
    \item	\code{sample} (take random selection of elements),
    \item	\code{table} (to make a table of counts),
  \end{itemize}
  as well as \code{length}, \code{mean}.

  \citep{Hurlbert71} showed that rarefaction generally overestimates
  the true estimated number of species ; can you corroborate this finding?
\clearpage

\section{Statistics}
  \R originated as a statistical package, and it is still predominantly used
  for this purpose.

  You can do virtually any statistical analysis in \R.

  As there exist many documents that may help you with statistical analyses in R,
  we will not deal with the subject here.

  Statistics is used just to show you how to use efficiently use \R,
  in cases where you have no clue where to begin!

\subsection{Using R in four steps}
  Suppose you want to perform a hierarchic clustering and plot the
  dendrogram of a multivariate data set.

  If you have never done that in \R here are the steps:
  \begin{enumerate}
    \item	Find a function that performs the requested task.

      for instance, use \code{help.search ("cluster")} to help you.

      Depending on the number of packages that you have installed, \R will
      list a number of possible functions whose help file contains the
      word "cluster".
      Use the function from the package \pkg{stats} (part of the core of R).
    \item	Open the help file (?<name-of-the-function>) and look up the
      syntax for this function.
      If you have no time to read it completely, at least read (part of)
      the section "Description", "Usage", and "Examples".
    \item	Try the examples in the help file. You may:
    \begin{itemize}
      \item	Try them all at once (example(<name-of-the-function>).
      \item	Alternatively, you may select the statements in the Examples
        section that look applicable to your problem, copy-paste them into
        your script file (Ctrl-C / Ctrl-V) and execute them; e.g. line by
        line. Chances are real that, if they are suited for your case, you
        will transform them anyway.
    \end{itemize}
    \item Transform a promising example so that it suits your problem.
  \end{enumerate}

\subsection{Exercise: multivariate statistics on the nematode species data.}
  Use \R to perform a multivariate statistical analysis of the nematode data.

  Beware: the nematode data have stations as columns and species as rows.

  \begin{itemize}
    \item	Perform a hierarchic clustering and plot the dendrogram
    \item	Perform a principal component analysis (PCA) and plot the results;
      you may also repeat the PCA analysis, with the first two stations removed!
  \end{itemize}
\clearpage

\section{Graphics}
  R has extensive graphical capabilities, and allows making simple (1-D, x-y),
  image-like (2-D) and perspective (3-D) figures.

  Try:
\begin{verbatim}
> demo(graphics)
> demo(image)
> demo(persp)
\end{verbatim}
  to obtain a display of R’s simple (1-D, x-y), image-like (2-D) and
  perspective (3-D) capabilities.

  Graphics are plotted in the figure window which floats independently from
  the other windows. If not already present, it is launched by writing
  (in windows):
\begin{verbatim}
> windows()
\end{verbatim}
  or
\begin{verbatim}
> x11()
\end{verbatim}

  A figure consists of a plot region surrounded by 4 margins, which are
  numbered clockwise, from 1 to 4, starting from the bottom.
  R distinguishes between:
  \begin{enumerate}
    \item high-level commands. By default, these create a new figure, e.g.
    \begin{itemize}
      \item	\code{hist}, \code{barplot}, \code{pie}, \code{boxplot}, ...
        (1-D plot)
      \item	\code{plot}, \code{curve}, \code{matplot}, \code{pairs},...
        ((x-y)plots)
      \item	\code{image}, \code{contour}, \code{filled.contour},...
        (2-D surface plots)
      \item	\code{persp}, \code{scatterplot3d},...			(3-D plots)
        \footnote{scatterplot3d is in \R-package \pkg{scatterplot3d} which
        has to be loaded first}.
    \end{itemize}

    \item	low-level commands that add new objects to an existing figure, e.g.
    \begin{itemize}
      \item	\code{lines}, \code{points}, \code{segments}, \code{polygon},
      \code{rect}, \code{text}, \code{arrows}, \code{legend}, \code{abline},
      \code{locator}, \code{rug}, ...
      These add objects within the plot region
      \item	\code{box}, \code{axis}, \code{mtext} (text in margin),
      \code{title}, ... which add objects in the plot margin
    \end{itemize}
    \item	graphical parameters that control the appearance of.
    \begin{itemize}
      \item plotting objects:

        \code{cex} (size of text and symbols), \code{col} (colors),
        \code{font}, \code{las} (axis label orientation), \code{lty}
        (line type), \code{lwd} (line width), \code{pch} (the type of points),…
      \item	graphic window:

        \code{mar} (the size of the margins), \code{mfrow} (the number of
        figures on a row), \code{mfcol} (number figures on a column),…
    \end{itemize}

  \end{enumerate}

\begin{verbatim}
> ?plot.default
> ?par
> ?plot.window
> ?points
\end{verbatim}
will open the help files, while
\begin{verbatim}
> example(plot.default)
> example(points)
\end{verbatim}
  will run the examples, displaying each new graph, after you have pressed
  "<ENTER>" (try it!)

\subsection{X-Y plots}
  A circle can be plotted by (x,y) points, where $x= r \cdot \cos(a)$  and
  $y= r \cdot \sin(a)$, with $a$ the angle, from 0 to $2\pi$, and $r$ the radius.
  In the following script, we first generate a sequence of angle values, a,
  from 0 to $2 \pi$, comprising 100 values (\code{length.out}) and then plot
  a circle with unit radius:
<<label=f1,include=FALSE>>=
a <- seq(0,2*pi, length.out=100)
plot(x=cos(a),y=sin(a))
@
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=f1fig,fig=TRUE,echo=FALSE>>=
<<f1>>
@
\end{center}
\caption{Simple figure with \code{plot} - see text for R-code}
\label{fig:f1}
\end{figure}

  As \code{plot} is a high-level command, it starts a new figure.

  By default, \R adds axes, and labels, and represents the (x,y) data as
  small dots (points). Note that the graph is not symmetrical.

  We will now make a more complex figure that resembles a "target face",
  e.g. for practicing archery or to throw darts.
  We first use the same command (\code{plot}) as above, but we add a
  number of graphical parameters that specify that:
  \begin{itemize}
    \item	Rather than dots, the points should be connected by
      lines (\code{type}).
    \item	The line should be twice as wide as the default (\code{lwd})
    \item	The x- and y-axes labels (\code{xlab},\code{ylab}) have to
      be empty
    \item	The axes and axes annotations (\code{axes}) are removed
    \item	The graph has to be symmetrical, i.e. the x/y aspect ratio = 1
      (\code{asp}).
  \end{itemize}
\begin{verbatim}
plot(cos(a),sin(a),type="l",lwd=2,xlab="",ylab="",axes=FALSE,
     asp=1)
\end{verbatim}

  To this figure, we can now add several low-level objects:
  \begin{itemize}
    \item	a series of lines, representing smaller and smaller circles
    (\code{lines}).
\begin{verbatim}
for (i in seq( 0.1,0.9,by=0.1)) lines(i*sin(a), i*cos(a))
\end{verbatim}
    \item	an innermost red polygon (\code{polygon}).
\begin{verbatim}
polygon(sin(a)*0.1,cos(a)*0.1,col="red")
\end{verbatim}

    \item	point marks as text labels, ranging from from 10 to 1
      (\code{text}). The closer to the centre, the higher the score
\begin{verbatim}
for (i in 1:10) text(x=0,y=i/10-0.025,labels=11-i,font=2)
\end{verbatim}

      Now two archers take 10 shots at the target face.

      We mimic their arrows by generating normally distributed (x,y)
      numbers, with mean=0 (the centre!) and where the experience of the
      archer is mimicked by the standard deviation.
      The more experienced, the closer the arrows will be to the centre,
      i.e. the lower the standard deviation.
    \item	R-statement \code{rnorm} generates normally distributed numbers;
      we need 20 of them, arranged as a matrix with 2 columns.
\begin{verbatim}
shots1 <- matrix(ncol=2, data=rnorm(n=20,sd=0.2))
shots2 <- matrix(ncol=2, data=rnorm(n=20,sd=0.5))
\end{verbatim}
     \item	The shots are added to the plot as points, colored darkblue
       (experienced archer) and darkgreen (beginners level).
       Note that we choose a 50\% enlarged point size (cex), and we choose
       a circular shaped point (\code{pch=16})
\begin{verbatim}
points(shots1,col="darkblue",pch=16,cex=1.5)
points(shots2,col="darkgreen",pch=16,cex=1.5)
\end{verbatim}

     \item	Finally, we add a legend, explaining who has done the shooting:
\begin{verbatim}
legend("topright",legend=c("A","B"),pch=16,
        col=c("darkblue","darkgreen"),pt.cex=1.5)
\end{verbatim}

<<label=f2g,include=FALSE,echo=FALSE>>=
plot(cos(a),sin(a),type="l",lwd=2,xlab="",ylab="",axes=FALSE,
 asp=1)
for (i in seq( 0.1,0.9,by=0.1)) lines(i*sin(a), i*cos(a))
polygon(sin(a)*0.1,cos(a)*0.1,col="red")
for (i in 1:10) text(x=0,y=i/10-0.025,labels=11-i,font=2)
shots1 <- matrix(ncol=2, data=rnorm(n=20,sd=0.2))
shots2 <- matrix(ncol=2, data=rnorm(n=20,sd=0.5))
points(shots1,col="darkblue",pch=16,cex=1.5)
points(shots2,col="darkgreen",pch=16,cex=1.5)
legend("topright",legend=c("A","B"),pch=16,
 col=c("darkblue","darkgreen"),pt.cex=1.5)
@

      Note that the legend text and the colors are inputted as a vector
      of strings, using the \code{c()} function (e.g. \code{c("A", "B")}).
  \end{itemize}

\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=f2fig,fig=TRUE,echo=FALSE>>=
<<f2g>>
@
\end{center}
\caption{Figure with several low-level objects - see text for R-code}
\label{fig:f2}
\end{figure}

\subsection{X-Y plots; conditional plotting}
  As a more sophisticated demonstration of the use of symbols in R-graphs,
  we work on a biological example, from the \R data set called "Orange".
  This data set contains the circumference (in mm, at breast height) measured
  at different ages for five orange trees.
  We start by looking at the data (only part is displayed).
<<>>=
head(Orange)
tail(Orange)
@

  and make a rough plot of circumference versus age:
<<label=f3,include=FALSE>>=
plot(Orange$age, Orange$circumference,xlab="age, days",
     ylab="circumference, mm", main= "Orange tree growth")
@
  (as Orange is a dataframe, columns can be addressed by their names,
  \code{Orange$age} and \code{Orange$circumference}).

  The output (figure) shows that there is a lot of scatter,
  which is due to the fact that the five trees did not grow at the same rate.
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=f3fig,fig=TRUE,echo=FALSE>>=
<<f3>>
@
\end{center}
\caption{Simple plot of the \code{orange} dataset - see text for R-code}
\label{fig:f3}
\end{figure}

  It is instructive to plot the relationship between circumference and
  age differently for each tree. In R, this is simple:
  we can make some graphical parameters (symbol types, colors, size,...)
  conditional to certain "factors".

  Factors play a very important part in the statistical applications of R;
  for our application, it suffices to know that the factors are integers,
  starting from 1.

  In the R-statement below, we simply use different symbols (\code{pch})
  and colors (\code{col}) for each tree: \code{pch=(15:20)[Orange$Tree]}
  means that, depending on the value of \code{Orange$Tree} (i.e. the tree
  number), the symbol (\code{pch})
  will take on the value 15 (tree=1), 16 (tree=2),... 20 (tree=5).
  \code{col=(1:5) [Orange$Tree]} does the same for the point color.
  The final statement adds a legend, positioned at the bottom, right.
<<label=f4,include=FALSE>>=
plot(Orange$age, Orange$circumference,xlab="age,
   days",ylab="circumference, mm", main= "Orange tree growth",
   pch=(15:20)[Orange$Tree],col=(1:5) [Orange$Tree],cex=1.3)
legend("bottomright",pch=15:20,col=1:5,legend=1:5)
@
  The output shows that tree number 5 grows fastest, tree number 1 is
  slowest growing. (note: it is also instructive to run the examples in
  the Orange help file. )
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=f4fig,fig=TRUE,echo=FALSE>>=
<<f4>>
@
\end{center}
\label{fig:f4}
\end{figure}


\subsection{Zooplankton growth rates}
  \code{Zoogrowth} from package \pkg{marelacTeaching} is a literature data set,
  compiled by Hansen et al. (1997)
  with measurements of maximal growth rates of zooplankton organisms as
  a function of body volume.

  Run the example for this data set (you will need to load package
  \pkg{marelacTeaching} first):
\begin{verbatim}
> require(marelacTeaching)
> example(Zoogrowth)
\end{verbatim}
\subsection{Images and contour plots}
  \R has some very powerful functions to create images and add contours.
  For example, the data set \code{Bathymetry} from the \pkg{marelac} package
  can be used to generate the bathymetry (and hypsometry) of the world oceans
  (and land):
<<label=bat,include=FALSE>>=
require(marelac)
image(Bathymetry$x,Bathymetry$y,Bathymetry$z,col=femmecol(100),
      asp=TRUE,xlab="",ylab="")
contour(Bathymetry$x,Bathymetry$y,Bathymetry$z,add=TRUE)
@
\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}
\begin{center}
<<label=batfig,fig=TRUE,echo=FALSE>>=
<<bat>>
@
\end{center}
\caption{Image plot of ocean bathymetry - see text for R-code}
\label{fig:f1}
\end{figure}

Note the use of "asp=TRUE", which maintains the aspect ratio.

\subsection{Plotting a mathematical function}
  Plot curves for mathematical functions are quickly generated with
  R-command "curve":
<<label=mat, include = FALSE>>=
curve(sin(3*pi*x))
curve(sin(3*pi*x),from=0,to=2,col="blue",
      xlab="x",ylab="f(x)",main="curve")
curve(cos(3*pi*x),add=TRUE,col="red",lty=2)
abline(h=0,lty=2)
legend("bottomleft",c("sin","cos"),text.col=c("blue","red"),lty=1:2)
@
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=matfig,fig=TRUE,echo=FALSE>>=
<<mat>>
@
\end{center}
\caption{Plotting a mathematical function - see text for R-code}
\label{fig:fmat}
\end{figure}

  The first command will plot the curve for $y=sin(3\cdot \pi \cdot x)$,
  using the default settings (Fig. left),
  while the next commands afirst draw a graph of $y=sin(3 \cdot \pi \cdot x)$,
  in blue (\code{col}), and for x values ranging between 0 and 2 (\code{from},
  \code{to}), adding a main title (\code{main}) and x- and y-axis labels
  (\code{xlab}, \code{ylab}) ($1^{st}$ sentence).

  The $2^{nd}$ R-sentence adds the function $y=cos(3\cdot \pi \cdot x)$,
  as a red (\code{col}) dashed line (\code{lty}).
  Note the use of parameter \code{add=TRUE}, as by default curve creates
  a new plot.

  The final statements adds the x-axis, i.e. a horizontal, dashed
  (\code{lty=2}), line (\code{abline}) at y=0 and a legend.

\subsection{Multiple figures}
  There are several ways in which to arrange multiple figures on a plot.
  \begin{enumerate}
    \item	The simplest is by specifying the number of figures on a row
      (mfrow) and on a column (mfcol):
\begin{verbatim}
> par(mfrow=c(3,2))
\end{verbatim}
       will arrange the next plots in 3 rows, 2 columns. Graphs will be
       plotted row-wise.
\begin{verbatim}
> par(mfcol=c(3,2))
\end{verbatim}
       will arrange the plots in 3 columns, 2 rows, in a columnwise sequence.
       Note that both mfrow and mfcol must be inputted as a vector. Try:
\begin{verbatim}
> par(mfrow=c(2,2))
> for ( i in 1:4) curve(sin(i*pi*x),0,1,main=i)
\end{verbatim}
    \item	R-function \code{layout} allows much more complex plot arrangements.
  \end{enumerate}
\subsection{Exercises}
\subsubsection{Simple curves}
  \begin{itemize}
    \item	Create a script file which draws a curve of the function
      $y=x^3 sin^2 (3\pi x)$ in the interval [-2, 2].
    \item	Make a curve of the function $y=1/cos(1+x^2)$ in the interval [-5,5].
  \end{itemize}

\subsubsection{Human population growth}
  The human population (N, millions of people) at a certain time t, can be
  described as a function of time (t),
  the initial population density at t=t0 (Nt0), the carrying capacity "K"
  and the rate of increase "a" by the following equation
  \footnote{This is the solution of a so-called logistic differential
  equation \citep{Verhulst38}}:
  \[
    N(t) = \frac{K}{{1 + [\frac{{K - N_{t0} }}{{N_{t0} }}]e^{ - a \cdot (t - t0)} }}
  \]

  For the US, the population density in 1900 (N0) was 76.1 million; the
  population growth can be described with parameter values:
  a=0.02 $yr^{-1}$, K =  500 million of people.

  Actual population values are:
\begin{verbatim}
1900  1910  1920  1930  1940  1950  1960  1970  1980
76.1	92.4	106.5	123.1	132.6	152.3	180.7	204.9	226.5
\end{verbatim}
  Tasks:
  \begin{enumerate}
    \item	Plot the population density curve as a thick line, using the
      US parameter values.
    \item	Add the measured population values as points. Finish the graph
      with titles, labels etc...
  \end{enumerate}

\subsubsection{Toxic ammonia}
  Ammonia nitrogen is present in two forms: the ammonium ion ($NH_4^+$)
  and unionized ammonia ($NH_3$).
  As ammonia can be toxic at sufficiently high levels, it is often desirable
  to know its concentration.

  The relative importance of ammonia, (the contribution of ammonia to
  total ammonia nitrogen,
  $NH_3/( NH_3+ NH_4^+)$) is a function of the proton concentration [$H^+$]
  and a parameter KN,
  the so-called stoichiometric equilibrium constant:
  \[
    p_{[NH_3 ]}  = \frac{{K_N }}{{K_N  + [H^ +  ]}}
  \]

  Tasks:
  \begin{itemize}
    \item	Plot the relative fraction of toxic ammonia to the total ammonia
      concentration as a function of pH,
      where $pH = -log_{10}([H^+])$ and for a temperature of 30$^{\circ}$C.
      Use a range of pH from 4 to 9.

      The value of KN is $8 10^{-10}$ at a temperature of 30$^{\circ}$C.
    \item	Add to this plot the relative fraction of ammonia at 0$^{\circ}$C;
      the value of KN at that temperature is 8 $10^{-11}$ mol $kg^{-1}$.
  \end{itemize}
\subsubsection{The iris data set}
  A famous data set that is part of \R is the "iris" data set (Fisher, 1936),
  which we will explore next.

  It gives measurements, in centimeters for sepal length and width and petal
  length and width, respectively,
  for 50 flowers of the species \emph{Iris setosa}, \emph{Iris versicolor}
  and \emph{Iris virginica}.

  Tasks:
  \begin{itemize}
    \item	Have a look at the data:
    \item	What is the class of the data set? why?
    \item	What are the dimensions of the data set? (number of rows, columns)
    \item	Produce a scatter plot of petal length against petal width;
      produce an informative title and labels of the two axes.
    \item	Repeat the same graph, using different symbol colours for
      the three species.
    \item	Add a legend to the graph. Copy-paste the result to a WORD
      document. If you do not have WORD, make a PDF file of the graph.
    \item	Create a box-and whisker plot for sepal length where the data
      values are split into species groups; use as template the first
      example in the "boxplot" help file.
    \item	Now produce a similar box-and whisker plot for all four
      morphological measurements, arranged in two rows and two columns.
      First specify the graphical parameter that arranges the plots two by two.
  \end{itemize}

\clearpage

\section{Matrix algebra}
  Matrix algebra is very simple in R. Practically everything is possible!
  Here are the most important R-functions that operate on matrices:
  \begin{itemize}
    \item	\%*\% 		Matrix multiplication
    \item	t(A)		transpose of A
    \item	diag(A)		diagonal of A
    \item	solve(A) 	inverse of A
    \item	solve(A,B)	solving Ax=B for x
    \item	eigen(A)	eigenvalues and eigenvectors for A
    \item	det(A)		determinant of A
  \end{itemize}
  For instance the following first inverts matrix A (\code{solve(A)}), and
  then multiplies the inverse with A , giving the unity matrix:

<<>>=
(A <-matrix(nrow=2,data=c(1,2,3,4)))
solve(A) %*% A
@

  whilst t(A) will transpose matrix A (interchange rows and columns).
<<>>=
t(A)
@
  The next set of statements will solve the linear system Ax=B for the
  unknown vector x:
<<>>=
B <- c(5,6)
solve(A,B)
@

  Finally, the eigenvalues and eigenvectors of A are estimated using
  R-function \code{eigen}.
  This function returns a list that contains both the eigenvalues (\code{$values})
  and the eigenvectors (\code{$vectors}), (the columns).
<<>>=
eigen(A)
@

\subsection{Exercises}
\subsubsection{Matrix algebra exercise 1}
  \begin{itemize}
    \item	Use R-function \code{matrix} to create the matrices called
      \code{A} and \code{B}:

      $ A = \left[ {\begin{array}{*{20}c}
         1 & 2 & 3  \\
         6 & 4 & 1  \\
         { - 2} & 1 & { - 1}  \\
      \end{array}} \right]$, $B=\left[ {\begin{array}{*{20}c}
         1 & 4 & 7  \\
         2 & 5 & 8  \\
         3 & 6 & 9  \\
      \end{array}} \right]
      $

    \item	Take the inverse of A and the transpose of A.
    \item	Multiply A with B.
    \item	Estimate the eigenvalues and eigenvectors of A.
    \item	For a matrix A, x is an eigenvector, and $\lambda$ the eigenvalue
    of a matrix A, if
    $A \cdot x =\lambda \cdot x$. Test it!
  \end{itemize}

\subsubsection{Matrix algebra exercise 2}
  \begin{itemize}
    \item	Create a matrix, called P:
      \[
      \left[ {\begin{array}{*{20}c}
         0 & {0.0043} & {0.1132} & 0  \\
         {0.9775} & {0.9111} & 0 & 0  \\
         0 & {0.0736} & {0.9534} & 0  \\
         0 & 0 & {0.0452} & {0.9804}  \\
      \end{array}} \right]
      \]

    \item	What is the value of the largest eigenvalue (the so-called
      dominant eigenvalue) and the corresponding eigenvector?.
    \item	Create a new matrix, T, which equals P, except for the first
      row, where the elements are 0.
    \item	Now estimate N= (I-T)-1, where I is the identity matrix
      \footnote{Note: this is a stage-model of a killer whale \citep{Caswell01}.
      The eigenvalue-eigenvectors estimate the rate of increase and stable
      age distribution, the matrix N contains the mean time spent in each stage.}.
  \end{itemize}

\subsubsection{System of linear equations}
  \begin{itemize}
    \item	Solve the following system of linear equations for the unknown xi:
\begin{verbatim}
3x1 + 4x2 + 5x3 = 0
6x1 + 2x2 + 7x3 = 5
7x1 +  x2       = 6
\end{verbatim}
    \item	You will first need to rewrite this problem in the form: Ax=B,
      where A contains the coefficients, x the unknowns, and B the right-hand
      side values. Then you solve the system, using R-function \code{solve}
    \item	Check the results (i.e. is Ax = B ?)
  \end{itemize}
\clearpage

\section{Roots of functions}
\subsection{Roots of a simple function}
  Suppose we want to solve the following problem: $\cos(x) = 2 \cdot x$ for x.

  Mathematically, we seek the root of the function $y= \cos(x) - 2 \cdot x$,
  this is the value of x for which $y = 0$.

  As the function is quite complex, it is not possible to find an exact solution
  (an explicit expression) for this root.

  It is always a good idea to plot the equation ($1^{st}$ line), and add the
  x-axis ($2^{nd}$ line).
\begin{verbatim}
curve(cos(x)-2*x,-10,10)
abline(h=0,lty=2)
\end{verbatim}

  This figure shows that there indeed exists a value x, for which y = 0.

  Now R-function \code{uniroot} can be used to locate this value.

  Functions that seek a root from a nonlinear equation generally work
  "iteratively", i.e. they move closer and closer to the root in
  successive steps (iterations).

  It is usually not feasible to find this root exactly, so it is approximated,
  i.e. up to a certain accuracy (tol, a very small number)
  \footnote{  More specifically: the root of $y=\cos (x)-2 \cdot x$ is the
  value x for which  |cos(x)-2*x | < tol or for which successive changes of
  x are < tol}.

  For the method to work, there should be at least one root in the interval.

  The statement below solves for the root; it returns several values, as a list.
<<>>=
(rr<-uniroot(f = function(x) cos(x)-2*x, interval=c(-10,10)))
@
  The most important value is the root itself (\code{$root}), which is
  0.45103686;

  the function value at the root was 3.66e-5, the function performed
  5 iterations.

  In this example, the function was simple enough to include it in the
  call to \code{uniroot}.

  The next chapter gives a more complex example from aquatic chemistry,
  where the equation to solve is significantly more complex.

  Finally, we add the root to the figure:
\begin{verbatim}
points(rr$root,0,pch=16,cex=2)
\end{verbatim}
<<root2, include=FALSE, echo=FALSE>>=
curve(cos(x)-2*x,-10,10)
abline(h=0,lty=2)
points(rr$root,0,pch=16,cex=2)
@
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=r1fig,fig=TRUE,echo=FALSE>>=
<<root2>>
@
\end{center}
\caption{Function drawn with \code{curve} and the root of the function
  plotted - see text for R-code}
\label{fig:f1}
\end{figure}
\subsection{Root of a complex function: solving for the pH}
  In aquatic systems, the buffering capacity of dissolved inorganic
  carbon (DIC) species:
  carbon dioxide ($CO_2$), bicarbonate ions ($HCO_3^-$) and carbonate
  ions ($CO_3^{2-}$) keep pH in a narrow range.
  The buffering capacity is measured by total alkalinity, TA, (defined below).

  If alkalinity and inorganic carbon concentration (DIC) are known, it is
  possible to calculate pH
  \citep{Zeebe03}, by solving the following equations for the unknown
  proton concentration $[H^+]$ \footnote{in practice, it is possible to
  merge these 3 equations such that only one equation is obtained,
  but this is neither didactically clearer nor computationally more efficient}.
  \[
  {\rm{ [HCO}}_{\rm{3}}^{\rm{ - }} ] = \frac{{K_{C1}  \cdot [H^ +  ]}}{{[H^ +  ] \cdot [H^ +  ] + K_{C1}  \cdot [H^ +  ] + K_{C1}  \cdot K_{C2} }} \cdot DIC
  \]
  \[
  {\rm{ [CO}}_{\rm{3}}^{{\rm{2 - }}} ] = \frac{{K_{C1}  \cdot K_{C2} }}{{[H^ +  ] \cdot [H^ +  ] + K_{C1}  \cdot [H^ +  ] + K_{C1}  \cdot K_{C2} }} \cdot DIC
  \]
  \[
  TA = 2{\rm{[CO}}_{\rm{3}}^{{\rm{2 - }}} ] + {\rm{[HCO}}_{\rm{3}}^{\rm{ - }} ] - {\rm{[H}}^ +  ]
  \]

  Here is how to solve for the proton concentration $[H^+]$ (or the pH value)
  in \R.

  The trick is to estimate alkalinity based on a guess of proton
  concentration, using equation (3) and compare that with the measured
  alkalinity value.
  If both are equal within the tolerance level, the proton concentration
  has been found.

  In the implementation below, the dissociation constants for carbonate
  (\code{kc1}, \code{kc2})
  and at salinity 0, temperature 20, and pressure 0 are calculated in
  \R s package \pkg{seacarb},
  which has to be loaded first (\code{require}).

  We then define a function whose root has to be solved (\code{pHfunction}).
  In this function we estimate total alkalinity, based on the guess of pH,
  the dissociation constants (\code{kc1},\code{kc2}) and the DIC concentration.
  The difference of this calculated alkalinity (\code{EstimatedAlk}) with
  the true alkalinity is then returned; if pH is correctly estimated, then
  true and estimated alkalinity will be equal, and the difference will be zero.
  So, to find the pH, we need to find the root of this function.

  Note that the conversion from pH to $[H^+]$ gives the proton concentration
  in $mol~kg^{-1}$.
  As the concentrations of the other substances are in $\mu mol~kg^{-1}$,
  we convert using a factor $10^6$.

  We restrict the region of the pH root in between 0 and 12 (which is more
  than large enough), and
  we set the tolerance (tol) to a very small number to increase precision.
\begin{verbatim}
require(seacarb)
kc1 <- K1(S=0,T=20,P=0)    # Carbonate k1
kc2 <- K2(S=0,T=20,P=0)    # Carbonate k2

pHfunction <- function(pH, kc1,kc2, DIC, Alkalinity )
{
   H    <- 10^(-pH)
   HCO3 <- H*kc1  /(H*kc1 + H*H + kc1*kc2)*DIC
   CO3  <- kc1*kc2 /(H*kc1 + H*H + kc1*kc2)*DIC

   EstimatedAlk  <- -H *1.e6  + HCO3 + 2*CO3
   return ( EstimatedAlk  - Alkalinity)
}
\end{verbatim}
<<echo=FALSE>>=
require(seacarb)
kc1 <- K1(S=0,T=20,P=0)    # Carbonate k1
kc2 <- K2(S=0,T=20,P=0)    # Carbonate k2

pHfunction <- function(pH, kc1,kc2, DIC, Alkalinity )
{
   H    <- 10^(-pH)
   HCO3 <- H*kc1  /(H*kc1 + H*H + kc1*kc2)*DIC
   CO3  <- kc1*kc2 /(H*kc1 + H*H + kc1*kc2)*DIC

   EstimatedAlk  <- -H *1.e6  + HCO3 + 2*CO3
   return ( EstimatedAlk  - Alkalinity)
}
@
<<>>=
uniroot(pHfunction, interval=c(0, 12), tol=1.e-20, kc1=kc1, kc2=kc2,
        DIC=2100, Alkalinity=2200)
@
\subsection{Exercises}
\subsubsection{Simple functions}
  \begin{itemize}
    \item	Find the root of the equation $e^x = 4 x^2$ in the interval [0,1].

      First draw the function curve.
    \item	Solve the equations $1000=y*(3+x)*(1+y)^4$ for y and with x
      varying over the range from 1 to 100. Plot the root as a function of x.

      Tip: first make a sequence of x-values, then loop over each "x"
      value, each time estimating the root and putting it in a vector.
  \end{itemize}
\subsubsection{Chemistry: pCO2 rises increase acidity}
  pH can also be estimated based on the measured alkalinity and $pCO_2$,
  the partial pressure of $CO_2$. To solve this equation, it is simplest to use
  another (equivalent) way to write the relationships between the DIC species:
  \[
  [HCO_3^ -  ] = K_{C1}  \cdot \frac{{[CO_2 ]}}{{[H^ +  ]}}
  \]
  \[
  [CO_3^{2 - } ] = K_{C2}  \cdot \frac{{[HCO_3^ -  ]}}{{[H^ +  ]}}
  \]

  $pCO_2$ relates to $[CO_2]$ through Henrys constant, Kh, which can also be
  estimated as a function of salinity, temperature and pressure, using
  R-package \pkg{seacarb}:
  \[
  pCO_2  = \frac{{[CO_2 ]}}{{Kh}}
  \]

  \begin{itemize}
    \item	Estimate the pH at equilibrium with alkalinity 2300 $\mu mol~ kg^{-1}$
      and the current $pCO_2$ of 360 ppm.

      Use package \code{seacarb} to estimate the dissociation constants and
      Henrys constants at temperature 20$^{\circ}$C, salinity 0, and
      pressure 0. 					(A: pH=8.19)
    \item	The Intergovernmental Panel on Climate Change predicts for 2100
      an atmospheric $CO_2$ concentration ranging between 490 and 1250 ppmv,
      depending on the socio-economic scenario (IPCC, 2007).
      These increases of $pCO_2$ make the water more acid. Make a plot of
      pH as a function of these increased atmospheric $pCO_2$ levels.
      (Assume that the $pCO_2$ of the ocean is at equilibrium with the
      atmospheric $pCO_2$).
      What is the maximal drop of pH ? 	(A: at $pCO_2$ of 1250 ppmv, pH=7.68).
  \end{itemize}

\clearpage

\section{Interpolation, smoothing}
  Interpolating and smoothing in \R can be done in several ways:
  \begin{itemize}
    \item	\code{approx} linearly interpolates through points
    \item	\code{spline} uses spline interpolation, which is smoother
    \item	\code{smooth.spline} smoothens data sets; this means that it
      does not connect the original points.
  \end{itemize}

  The use of these functions is exemplified in the following script and
  corresponding output:
<<>>=
x <- 1:10
y <- c(9,8,6,7,5,8,9,6,3,5)
@
<<smooth,include=FALSE>>=
plot(x,y,pch=16,cex=2,main="interpolation,smoothing")
lines (spline(x,y, n=100),lty=1)
points(approx(x,y, xout=seq(1,10,0.1)),pch=1)
lines (smooth.spline(x,y),lty=2)
legend("bottomleft",lty=c(1,NA,2),pch=c(NA,1,NA),
       legend=c("spline","approx","smooth.spline"))
@
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=smoothfig,fig=TRUE,echo=FALSE>>=
<<smooth>>
@
\end{center}
\caption{smoothing and interpolation - see text for R-code}
\label{fig:smooth}
\end{figure}

\subsection{Curve fitting}
  \R also has several curve fitting procedures. Depending on whether the
  function to be fitted is linear, or non-linear, you may use:
  \begin{itemize}
    \item	\code{lm} and \code{glm} for fitting linear models and
      generalised linear models
    \item	\code{nls}, \code{nlm}, \code{optim}, \code{constrOptim} for
      nonlinear models.
  \end{itemize}

  As an example, we now fit the US population density values, at
  10-year intervals, with the logistic growth model (see previous chapter).
  The model was:
  $N(t)=\frac{{K}}{{1+[\frac{K-N_{t0}}{N_{t0}}]\cdot e^{-a\cdot(t-t_0)}}}$ ,
  and the data:
\begin{verbatim}
1900	1910	1920	1930	1940	1950	1960	1970	1980
76.1	92.4	106.5	123.1	132.6	152.3	180.7	204.9	226.5
\end{verbatim}
  We start by inputting the data.
<<>>=
year<- seq(1900,1980,by=10)
pop <- c(76.1,92.4,106.5,123.1,132.6,152.3,180.7,204.9,226.5)
@

  The simplest method for non-linear curve fitting is by using R-function
  \code{nls}.

  This functions requires as input the formula (\code{y ~ f (x,parameters)})
  and starting values of the parameters.

  In the example, y are the population values, f is the logistic
  growth formulation.

  As starting conditions, we use: K=500,N0=76.1,a=0.02.
<<>>=
fit <- nls(pop~K/(1+(K-N0)/N0*exp(-a*(year-1900))),
           start=list(K=500,N0=76.1,a=0.02))
@
  We end the fitting by printing a summary of the fitting parameters, which
  shows the estimates of the parameters and their standard errors.
  Clearly, it is not possible to obtain reliable estimates of the value of
  K based on the data.
<<>>=
summary(fit)
@

  The values of the coefficients themselves are retrieved using
  R-function \code{coef}.
<<>>=
(cc<-coef(fit))
@
\subsection{Exercises}
\subsubsection{Smoothing}
  An anemometer measures wind-velocity at three hourly intervals.
  On a certain day, these velocities are: 5,6,7,9,4,6,3,7,9 at time 0, 3, ... 24
  o'clock respectively. In order to estimate air-sea exchange, we need
  hourly measures.

  Tasks:
  \begin{itemize}
    \item	Interpolate the three-hourly measurements to hourly measurements.
    \item	Make a plot of the interpolated values
  \end{itemize}
\subsubsection{Fitting}
  Primary production is measured by $^{14}C$ incubations from phytoplankton
  samples, at different light intensities.

The data are:
<<>>=
ll <- c(0.,1,10,20,40,80,120,160,300,480,700)
pp <- c(0.,1,3,4,6,8,10,11,10,9,8)
@
  Fit the resulting production estimates (pp), as a function of light
  intensity (ll) with the 3-parameter Eilers-Peeters equation.
  The primary production is calculated as:
  \[
    pp = p\max  \cdot \frac{{2 \cdot (1 + \beta ) \cdot {\raise0.7ex\hbox{$I$} \!\mathord{\left/
    {\vphantom {I {Iopt}}}\right.\kern-\nulldelimiterspace}
    \!\lower0.7ex\hbox{${Iopt}$}}}}{{({\raise0.7ex\hbox{$I$} \!\mathord{\left/
    {\vphantom {I {Iopt}}}\right.\kern-\nulldelimiterspace}
    \!\lower0.7ex\hbox{${Iopt}$}})^2  + 2 \cdot \beta  \cdot {\raise0.7ex\hbox{$I$} \!\mathord{\left/
     {\vphantom {I {Iopt}}}\right.\kern-\nulldelimiterspace}
    \!\lower0.7ex\hbox{${Iopt}$}} + 1}}
  \]

  where I is light and pmax, $\beta$ and Iopt are parameters.
  \begin{itemize}
    \item	First plot primary production (pp) versus light  (ll). Use
      large symbols.
    \item	Then use R-function \code{nls} to fit the model to the data
    \item	Add the best-fit line to the graph. (note: use \code{coef} to
      retrieve the best parameter values).
  \end{itemize}
\clearpage

\section{Differential equations}
  Differential equations express the rate of change of a constituent (C)
  along one or more dimensions, usually time and/or space.

  Consider the following set of two differential equations:
  \[
    \begin{array}{l}
     \frac{{dA}}{{dt}} = r \cdot (x - A) - k \cdot A \cdot B \\
     \frac{{dB}}{{dt}} = r \cdot (y - B) + k \cdot A \cdot B \\
    \end{array}
  \]

  \begin{itemize}
    \item	A and B are called differential variables (or state variables),
    \item $\frac{{dA}}{{dt}}$ is the derivative (or the rate of change),
    \item	r,x,y and k are parameters (constants).
  \end{itemize}

  To solve sets of differential equations in \R we define a function (here
  called \code{model}), and which has as input the time (\code{t}), the
  values of the state variables (\code{state}) and the values of the
  parameters (\code{pars}). This function simply calculates the rate of
  change of the state variables (\code{dA} and \code{dB}) and returns those
  as a list.

  The R-statement \code{with (as.list (c(state,pars)),{ }} ensures that
  the state variables and parameters can be addressed by their names.
\begin{verbatim}
model <- function(t,state,pars)
{
with (as.list(c(state,pars)),
{
 dA <- r*(x-A)-k*A*B
 dB <- r*(y-B)+k*A*B
 return (list(c(dA,dB)))
 }
     )
}
\end{verbatim}

  Before we can solve this model, we
  \begin{itemize}
    \item	generate a sequence of time values at which we want output
      (\code{times}),
    \item	assign initial conditions to the state variables (\code{state}) and
    \item	give values to the parameters (\code{parms}):
  \end{itemize}
\begin{verbatim}
times <- seq(0,300,1)
state <- c(A=1,B=1)
parms <- c(x =1, y = 0.1, k = 0.05, r = 0.05)
\end{verbatim}
  The model can now be solved. To do so, we use \R’s integration routine
  \code{ode}; which can be found in \R-package \pkg{deSolve}. This package
  is loaded first.
\begin{verbatim}
require(deSolve)
\end{verbatim}
  At each time \code{t}, \code{ode} will call function \code{model}, with
  the current values of the state variables and the parameter values.

  The output is stored in a data.frame, called \code{out}.
\begin{verbatim}
out   <- as.data.frame(ode(state,times,model,parms))
\end{verbatim}
<<echo=FALSE>>=
model <- function(t,state,pars)
{
with (as.list(c(state,pars)),
{
 dA <- r*(x-A)-k*A*B
 dB <- r*(y-B)+k*A*B
 return (list(c(dA,dB)))
 }
     )
}
times <- seq(0,300,1)
state <- c(A=1,B=1)
parms <- c(x =1, y = 0.1, k = 0.05, r = 0.05)
out   <- as.data.frame(ode(state,times,model,parms))
@
  All we need to do now is to plot the model output. Before we do so, we
  have a look at data.frame out:
<<>>=
head(out)
@
  The data are arranged in three columns: first the time, then the
  concentrations of A and B. As out is a data frame we can extract the
  data using their names (\code{out$time}, \code{out$A}, \code{out$B}).

  Before plotting the model output, the range of concentrations of
  substances A and B is estimated; this is used to set the limits of
  the y-axis (\code{ylim}).

  \R-function \code{plot} creates a new plot; \code{lines} adds a line
  to this plot; \code{lty} selects a line type; \code{lwd=2} makes the
  lines twice as thick as the default.  Finally a \code{legend} is added.
\begin{verbatim}
ylim   <- range(c(out$A,out$B))
plot(out$time,out$A,xlab="time",ylab="concentration",
      lwd=2,type="l",ylim=ylim,main="model")
lines(out$time,out$B,lwd=2,lty=2)
legend("topright",legend=c("A","B"),lwd=2,lty=c(1,2))
\end{verbatim}
<<ode,include=FALSE,echo=FALSE>>=
ylim   <- range(c(out$A,out$B))
plot(out$time,out$A,xlab="time",ylab="concentration",
      lwd=2,type="l",ylim=ylim,main="model")
lines(out$time,out$B,lwd=2,lty=2)
legend("topright",legend=c("A","B"),lwd=2,lty=c(1,2))
@
\setkeys{Gin}{width=0.4\textwidth}
\begin{figure}
\begin{center}
<<label=lvfig,fig=TRUE,echo=FALSE>>=
<<ode>>
@
\end{center}
\caption{ode model - see text for R-code}
\label{fig:smooth}
\end{figure}
\subsection{Exercises}
\subsubsection{Lotka-volterra model}
  \begin{itemize}
    \item	Write a script file that solves the following system of ODEs
      \footnote{ The Lotka-Volterra models are a famous type of models that
      either describe predator-prey interactions or competitive interactions
      between two species.
      A.J. Lotka and V. Volterra formulated the original model in the 1920's
      almost simultaneously (\citep{Lotka25}, \citep{Volterra26}).}:
  \[
    \begin{array}{l}
     \frac{{dx}}{{dt}} = a \cdot x \cdot (1 - \frac{x}{K}) - b \cdot x \cdot y \\
     \frac{{dy}}{{dt}} = g \cdot b \cdot x \cdot y - e \cdot y \\
    \end{array}
  \]

      for initial values x=300,y=10 and parameter values: a=0.05, K=500,
      b=0.0002, g=0.8, e=0.03
    \item	Make three plots, one for x and one for y as a function of time,
      and one plot expressing y as a function of x (this is called a
      phase-plane plot).  Arrange these plots in 2 rows and 2 columns.
    \item	Now run the model with other initial values (x=200, y=50); add
      the (x,y) trajectories to the phase-plane plot
  \end{itemize}

\subsubsection{Butterfly}
  The Lorenz equations \citep{Lorenz63} were the first chaotic system
  of differential equations to be discovered. They are three differential
  equations that were derived to represent idealized behavior of the
  earth’s atmosphere.
  \[
    \begin{array}{l}
     \frac{{dx}}{{dt}} =  - \frac{8}{3} \cdot x + y \cdot z \\
     \frac{{dy}}{{dt}} =  - 10 \cdot (y - z) \\
     \frac{{dz}}{{dt}} =  - x \cdot y + 28y - z \\
    \end{array}
  \]

  \begin{itemize}
    \item	It takes about 10 lines of R-code to generate the solutions
      and plot them.
    \item	Function \code{scatterplot3d} from the package scatterplot3d
      generates 3-D scatterplots. Can you recreate the following
      "butterfly" ? Use as initial conditions x=y=z=1; create output for
      a time sequence ranging from 0 to 100, and with a time step of 0.005.
  \end{itemize}

\begin{center}
<<fig=TRUE,echo=FALSE>>=
require(scatterplot3d)
model<-function(t,state,parameters)
  {
  with(as.list(c(state)),{

    dx     <- -8/3*x+y*z
    dy     <- -10*(y-z)
    dz     <- -x*y+28*y-z

    list(c(dx,dy,dz))            })

 }  # end of model

state <-c(x=1, y=1, z=1)
times <-seq(0,100,0.005)
out   <-as.data.frame(lsoda(state,times,model,0))
scatterplot3d(out$x,out$y,out$z,type="l",
        main="Lorenz butterfly",ylab="",grid=FALSE,box=FALSE)
@
\end{center}


\section{Finally}

\subsection{The questions}
  These lecture notes have been generated with LaTeX and making use of
  R-package \pkg{Sweave} \citep{Leisch02}, which allows to merge LaTeX
  with R-code.

  If you do not like the layout a PDF version ("ScientificComputing.pdf")
  made with WORD (Microsoft) can be found in the /inst/lecture subdirectory
  of package \pkg{marelacTeaching}.

\subsection{The answers}
  The answers to the questions in this course are present as an R-vignette
  in package "marelacTeaching".
  From within R, type:
\begin{verbatim}
> vignette("Answers")
\end{verbatim}
  Or, you can find the file "Answers.pdf" in the /inst/doc subdirectory
  of package "marelacTeaching".


\bibliography{vignettes}

\end{document}
