Why statistics?

From Motulsky (2014) Chapter 1:

We tend to jump to conclusiong
We tend to be overconfident
We see patterns in random data
Coincidences are common
Incorrect intuition about probability
…

What is Probability?

Long-term Frequency — Frequentist

(# of head)/(# of coin toss)
predicted from model
estimated from data
probability of Higgs boson exists?
- We need to create the universe many times…

Strength of Belief — Bayesian

subjective probability
quantification of uncertainty/ignorance (Cox, 1946)
- sum rule: $prob(X) + prob(\bar{X}) = 1$
- product rule: $prob(X,Y) = prob(X) \times prob(Y)$

Probability and Statistics

\[\begin{matrix} & \mbox{Population (Model)}& \\ \mbox{Probability} & \downarrow \quad \uparrow & \mbox{Statistics} \\ & \mbox{Sample (Data)} & \end{matrix}\]

R Language

R was developed by R. Ihaka and R. Gentleman (1996) as a successor of S language (Chambers & Hastie, 1988).

Pros:
- Free and open-source, for Linux, Mac, and Windows
- Easy to start using
- Sample data sets
- Programming for custom use
- Many packages contributed by statisticians
Cons:
- Fortran legacy
- Object-oriented like, but not really
- Classis standard graphics

Running R

First, download R package for your computer from R Project site and install as instructed there.
The basic way is to run a command “r” on your computer terminal:

$ r

It is more convenient to use an integrated develpment environment (IDE) like RStudio
In RStudito, you can
- type in commands in the “Console” window
- run scripts in “Terminal” window
- run code chunks in R Markdown files within the editor window

R Markdown

R Markdown combines explanation, code, and result in a same file.

“File >New File >R Markdown…” menu to create a new file
You can embed LaTeX-like mathematical formula by Mathjax
-  or $ $ for inline, e.g., $\sum_{i=1}^n x_i$ or $\int_0^\infty f(x)dx$.
- \[ \] or $$ $$ for a separate line, e.g., \[\sum_{i=1}^n x_i\] or \[\int_0^\infty f(x)dx.\]
To embed your R code, “Option+Command+I” or “Insert” icon to make a code chunk

To run a code chunk:
- “Shift+Comman+Enter” or the triangle in the uppper-right corner
- “Comman+Enter” for each line
- “Run All” to run all the code chunks in the file
“File >Knit Document” menu, or “Knit” icon, to produce an html file.

R Notebook

R Notebook a new type of R Markdown for embedding R Markdown code in an html file. This file is an example of an R Notebook.

“File >New File >R Notebook” menu to create a new file, which includes
output: html_notebook
in the file header.
“Prevew >Preview Notebook” button to show a page preview.
Everytime you save your “file.Rmd”, R Studio creates “file.nb.html”.

R as a calculator

You can just type in numbers and operators in the Console window.

1+1

[1] 2

2*3

[1] 6

4/5

[1] 0.8

Either ** or ^ for power.

2**3

[1] 8

2^3

[1] 8

Many math functions and some constants are pre-defined.

exp(1)

[1] 2.718282

sin(pi/2)

[1] 1

Variables

You can use <- or = to assign a value to a variable.

x <- 1
y = 2
x+y

[1] 3

x/y

[1] 0.5

z = x+y
z

[1] 3

Getting information

ls() lists objects in your work space.
rm() removes objects.

ls()

 [1] "a"        "ai"       "alpha"   
 [4] "av"       "b"        "b0"      
 [7] "b0hat"    "b1"       "b1hat"   
[10] "c"        "CI"       "d2norm"  
[13] "dmn"      "eps"      "ess"     
[16] "f"        "F"        "fpbc"    
[19] "fpfd"     "fpt"      "fptt"    
[22] "g"        "i"        "interval"
[25] "j"        "k"        "L"       
[28] "lambda"   "m"        "M"       
[31] "miss"     "mu"       "n"       
[34] "nb"       "nu"       "ones"    
[37] "p"        "phat"     "pp"      
[40] "pval"     "q"        "r2"      
[43] "rss"      "s"        "S"       
[46] "SE"       "sig"      "sigma"   
[49] "t"        "T"        "t1"      
[52] "tss"      "tstar"    "tt"      
[55] "V"        "w"        "W"       
[58] "x"        "X"        "x1"      
[61] "X1"       "x2"       "X2"      
[64] "xbar"     "xq"       "XY"      
[67] "y"        "ybar"     "z"

rm(x,y)
ls()

 [1] "a"        "ai"       "alpha"   
 [4] "av"       "b"        "b0"      
 [7] "b0hat"    "b1"       "b1hat"   
[10] "c"        "CI"       "d2norm"  
[13] "dmn"      "eps"      "ess"     
[16] "f"        "F"        "fpbc"    
[19] "fpfd"     "fpt"      "fptt"    
[22] "g"        "i"        "interval"
[25] "j"        "k"        "L"       
[28] "lambda"   "m"        "M"       
[31] "miss"     "mu"       "n"       
[34] "nb"       "nu"       "ones"    
[37] "p"        "phat"     "pp"      
[40] "pval"     "q"        "r2"      
[43] "rss"      "s"        "S"       
[46] "SE"       "sig"      "sigma"   
[49] "t"        "T"        "t1"      
[52] "tss"      "tstar"    "tt"      
[55] "V"        "w"        "W"       
[58] "X"        "x1"       "X1"      
[61] "x2"       "X2"       "xbar"    
[64] "xq"       "XY"       "ybar"    
[67] "z"

getwd() to get or set the working directory

getwd()

[1] "/Users/doya/Dropbox (OIST)/R/StatisticalMethods"

help() or ? for help. ?? for related topics.

?sin
??sin

You can also use “Help” panel of RStudio.

Vectors

c( ) cancatenates objects to make a vector.

a = c(1,2,3)
a

[1] 1 2 3

You can access each component by [ ].

a[1]

[1] 1

a[c(2,3)]

[1] 2 3

Operators work component-wise. There are many functions that work with vectors.

a+a

[1] 2 4 6

3*a

[1] 3 6 9

a*a

[1] 1 4 9

mean(a)

[1] 2

exp(a)

[1]  2.718282  7.389056 20.085537

Regular vectors can be created by : or seq().

b = 1:5
b

[1] 1 2 3 4 5

c = seq(0,10,2)
c

[1]  0  2  4  6  8 10

You can also create an empty array and later assign values.

d = array(dim=10)
d[1] = 1
d

 [1]  1 NA NA NA NA NA NA NA NA NA

Matrix

You can create a new matrix by

A = matrix(0, 2, 3)
A

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0

or by combining vectors

B = rbind(1:3, 4:6)
B

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6

C = cbind(1:3, 4:6)
C

     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

You can reshaping a vector into a matrix.

D = 1:6
dim(D) = c(2, 3)
D

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

You can access the components by [ , ]

D[1,2]

[1] 3

D[1:2,2:3]

     [,1] [,2]
[1,]    3    5
[2,]    4    6

D[1,]

[1] 1 3 5

D[,1]

[1] 1 2

Most operators and functions work component-wise.

A = rbind(1:2, 3:4)
A

     [,1] [,2]
[1,]    1    2
[2,]    3    4

A+A

     [,1] [,2]
[1,]    2    4
[2,]    6    8

A*A

     [,1] [,2]
[1,]    1    4
[2,]    9   16

sin(A)

         [,1]       [,2]
[1,] 0.841471  0.9092974
[2,] 0.141120 -0.7568025

To perform matrix multiplication, we must use %*%

A %*% A

     [,1] [,2]
[1,]    7   10
[2,]   15   22

Inverse is given by solve() function

Ainv = solve(A)
Ainv

     [,1] [,2]
[1,] -2.0  1.0
[2,]  1.5 -0.5

A %*% Ainv

     [,1]         [,2]
[1,]    1 1.110223e-16
[2,]    0 1.000000e+00

Use t() for a transpose matrix.

t(A)

     [,1] [,2]
[1,]    1    3
[2,]    2    4

List

A list is a collection of components of arbitrary types, e.g., numbers, matrices, characters, or functions.

L = list(name="OIST", address=c("Tancha","Onna","Okinawa"), age=6)
L[1]

$name
[1] "OIST"

a component can be specified by its index or name

L[[1]]

[1] "OIST"

L[["name"]]

[1] "OIST"

L$name

[1] "OIST"

L$add

[1] "Tancha"  "Onna"    "Okinawa"

Script

You can store multiple lines of R commands as a script for repeated use.

Choose “File >New File >R Script” to make a new file. Type in, for example:

# bern.R: Bernoulli distribution of 0 and 1
p = 0.5  # probability of 1
if (runif(1) < p){  # sample a random number uniformly in 0 and 1
    x = 1  # 1 if that falls below p
  }else{ 
    x = 0  # 0 otherwise
  }
print(x)

# starts a comment, which is skipped by R but helpful for humans, including yourself, to understand what the code does and what the parameter means.

Choose “File >Save As…” and name it ‘bern.R’

Go to the directory where the file was saved, start R, and type in:

source("bern.R")

[1] 0

Defining your own function

For repeating a procedure with different parameters, you would define your own function, in the following convention:

bern <- function(p){  # p: probability of 1
  if (runif(1) < p){  # sample a uniform random number
    x = 1    # 1 if that falls below p
  }else{
    x = 0    # 0 otherwise
  }
  return(x)  
}

You can call it with any parameter:

bern(0.2)

[1] 0

bern(0.8)

[1] 1

Loops

Try making bionomial distribution by adding 1s of n bernoulli samples using for loop:

binom = function(n, p) { # n:samples, p:probability
  y = 0   # prepare a counter
  for (i in 1:n) {
    x = bern(p)   # 1 or 0
    # print(x)    # uncomment to check each x
    y = y + x
  }
  return(y)
}

You can test it with different parameters:

binom(5, 0.5)

[1] 2

binom(10, 0.3)

[1] 4

Plotting

Let us take samples from a binomial distribution

m = 1000  # number of samples
n = 10
p = 0.5
x = array(dim=m)   # prepare an empty array
for (i in 1:m) {
  x[i] = binom(n, p)  # between 0 to n
}

You can visualize the raw data by plot() function.

plot(x)

or make a histogram by hist().

hist(x)

You can plot any function by specifying x and y for plot() function.

X = seq(0, n)
plot(X, dbinom(X, n, p))

To overlay a curve on an existing plot, use lines() function.

hist(x, freq=FALSE)  # normalize to sum up to 1
lines(X, dbinom(X, n, p))

3D Plots

R has a simple 3D plot function persp().

x = seq(-5, 5, 0.1)
y = seq(0.5, 3, 0.1)
f = function(x,y){exp(-x**2/(2*y**2))/y}
z = outer(x, y, f)
persp(x, y, z, theta=45, phi=30)

There are fancier 3D packages like plot3D and rgl.

To be covered in later chapters

probability distributions
data frame
pre-loaded data sets
reading data from a file

LS0tCnRpdGxlOiAiMS4gSW50cm9kdWN0aW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIyBXaHkgc3RhdGlzdGljcz8KCkZyb20gTW90dWxza3kgKDIwMTQpIENoYXB0ZXIgMToKCiogV2UgdGVuZCB0byBqdW1wIHRvIGNvbmNsdXNpb25nICAKKiBXZSB0ZW5kIHRvIGJlIG92ZXJjb25maWRlbnQgIAoqIFdlIHNlZSBwYXR0ZXJucyBpbiByYW5kb20gZGF0YSAgCiogQ29pbmNpZGVuY2VzIGFyZSBjb21tb24gIAoqIEluY29ycmVjdCBpbnR1aXRpb24gYWJvdXQgcHJvYmFiaWxpdHkgIAoqIC4uLiAgCgojIyBXaGF0IGlzIFByb2JhYmlsaXR5PwoKIyMjIExvbmctdGVybSBGcmVxdWVuY3kg4oCUIEZyZXF1ZW50aXN0CgoqICgjIG9mIGhlYWQpLygjIG9mIGNvaW4gdG9zcykgIAoqIHByZWRpY3RlZCBmcm9tIG1vZGVsICAKKiBlc3RpbWF0ZWQgZnJvbSBkYXRhICAKKiBwcm9iYWJpbGl0eSBvZiBIaWdncyBib3NvbiBleGlzdHM/ICAKICAgICogV2UgbmVlZCB0byBjcmVhdGUgdGhlIHVuaXZlcnNlIG1hbnkgdGltZXMuLi4KICAgIAojIyMgU3RyZW5ndGggb2YgQmVsaWVmIOKAlCBCYXllc2lhbgoKKiBzdWJqZWN0aXZlIHByb2JhYmlsaXR5ICAKKiBxdWFudGlmaWNhdGlvbiBvZiB1bmNlcnRhaW50eS9pZ25vcmFuY2UgKENveCwgMTk0NikKICAgICsgc3VtIHJ1bGU6ICRwcm9iKFgpICsgcHJvYihcYmFye1h9KSA9IDEkICAKICAgICsgcHJvZHVjdCBydWxlOiAkcHJvYihYLFkpID0gcHJvYihYKSBcdGltZXMgcHJvYihZKSQgIAoKIyMgUHJvYmFiaWxpdHkgYW5kIFN0YXRpc3RpY3MKCiQkXGJlZ2lue21hdHJpeH0gJiBcbWJveHtQb3B1bGF0aW9uIChNb2RlbCl9JiBcXApcbWJveHtQcm9iYWJpbGl0eX0gJiBcZG93bmFycm93IFxxdWFkIFx1cGFycm93ICYgXG1ib3h7U3RhdGlzdGljc30gXFwKJiBcbWJveHtTYW1wbGUgKERhdGEpfSAmClxlbmR7bWF0cml4fSQkCgojIFIgTGFuZ3VhZ2UKClIgd2FzIGRldmVsb3BlZCBieSBSLiBJaGFrYSBhbmQgUi4gR2VudGxlbWFuICgxOTk2KSBhcyBhIHN1Y2Nlc3NvciBvZiBTIGxhbmd1YWdlIChDaGFtYmVycyAmIEhhc3RpZSwgMTk4OCkuCgoqIFByb3M6CiAgICAqIEZyZWUgYW5kIG9wZW4tc291cmNlLCBmb3IgTGludXgsIE1hYywgYW5kIFdpbmRvd3MKICAgICogRWFzeSB0byBzdGFydCB1c2luZwogICAgKiBTYW1wbGUgZGF0YSBzZXRzCiAgICAqIFByb2dyYW1taW5nIGZvciBjdXN0b20gdXNlCiAgICAqIE1hbnkgcGFja2FnZXMgY29udHJpYnV0ZWQgYnkgc3RhdGlzdGljaWFucwogICAgCiogQ29uczoKICAgICogRm9ydHJhbiBsZWdhY3kKICAgICogT2JqZWN0LW9yaWVudGVkIGxpa2UsIGJ1dCBub3QgcmVhbGx5IAogICAgKiBDbGFzc2lzIHN0YW5kYXJkIGdyYXBoaWNzCiAgICAKIyMgUnVubmluZyBSCgoqIEZpcnN0LCBkb3dubG9hZCBSIHBhY2thZ2UgZm9yIHlvdXIgY29tcHV0ZXIgZnJvbSBbUiBQcm9qZWN0IHNpdGVdKGh0dHBzOi8vd3d3LnItcHJvamVjdC5vcmcpCmFuZCBpbnN0YWxsIGFzIGluc3RydWN0ZWQgdGhlcmUuCgoqIFRoZSBiYXNpYyB3YXkgaXMgdG8gcnVuIGEgY29tbWFuZCAiciIgb24geW91ciBjb21wdXRlciB0ZXJtaW5hbDoKYGBgCiQgcgpgYGAKCiogSXQgaXMgbW9yZSBjb252ZW5pZW50IHRvIHVzZSBhbiBpbnRlZ3JhdGVkIGRldmVscG1lbnQgZW52aXJvbm1lbnQgKElERSkgbGlrZSBbUlN0dWRpb10oaHR0cHM6Ly93d3cucnN0dWRpby5jb20pCgoqIEluIFJTdHVkaXRvLCB5b3UgY2FuCiAgICAqIHR5cGUgaW4gY29tbWFuZHMgaW4gdGhlICJDb25zb2xlIiB3aW5kb3cKICAgICogcnVuIHNjcmlwdHMgaW4gIlRlcm1pbmFsIiB3aW5kb3cKICAgICogcnVuIGNvZGUgY2h1bmtzIGluIFIgTWFya2Rvd24gZmlsZXMgd2l0aGluIHRoZSBlZGl0b3Igd2luZG93CgojIyBSIE1hcmtkb3duCgpbUiBNYXJrZG93bl0oaHR0cHM6Ly9ybWFya2Rvd24ucnN0dWRpby5jb20pIGNvbWJpbmVzIGV4cGxhbmF0aW9uLCBjb2RlLCBhbmQgcmVzdWx0IGluIGEgc2FtZSBmaWxlLiAgCgoqICJGaWxlID5OZXcgRmlsZSA+UiBNYXJrZG93bi4uLiIgbWVudSB0byBjcmVhdGUgYSBuZXcgZmlsZQoKKiBZb3UgY2FuIGVtYmVkIExhVGVYLWxpa2UgbWF0aGVtYXRpY2FsIGZvcm11bGEgYnkgW01hdGhqYXhdKGh0dHBzOi8vd3d3Lm1hdGhqYXgub3JnLykgIAoKICAgICogYFwoIFwpYCBvciBgJCAkYCBmb3IgaW5saW5lLCBlLmcuLCAKICAgIFwoXHN1bV97aT0xfV5uIHhfaVwpIG9yICRcaW50XzBeXGluZnR5IGYoeClkeCQuICAgCiAgICAKICAgICogYFxbIFxdYCBvciBgJCQgJCRgIGZvciBhIHNlcGFyYXRlIGxpbmUsIGUuZy4sIAogICAgXFtcc3VtX3tpPTF9Xm4geF9pXF0gb3IgJCRcaW50XzBeXGluZnR5IGYoeClkeC4kJCAgIAoKKiBUbyBlbWJlZCB5b3VyIFIgY29kZSwgIk9wdGlvbitDb21tYW5kK0kiIG9yICJJbnNlcnQiIGljb24gdG8gbWFrZSBhICpjb2RlIGNodW5rKgoKYGBge3J9CgpgYGAKCiogVG8gcnVuIGEgY29kZSBjaHVuazoKICAgICogIlNoaWZ0K0NvbW1hbitFbnRlciIgb3IgdGhlIHRyaWFuZ2xlIGluIHRoZSB1cHBwZXItcmlnaHQgY29ybmVyICAgIAogICAgKiAiQ29tbWFuK0VudGVyIiBmb3IgZWFjaCBsaW5lCiAgICAqICJSdW4gQWxsIiB0byBydW4gYWxsIHRoZSBjb2RlIGNodW5rcyBpbiB0aGUgZmlsZQoKKiAiRmlsZSA+S25pdCBEb2N1bWVudCIgbWVudSwgb3IgIktuaXQiIGljb24sIHRvIHByb2R1Y2UgYW4gaHRtbCBmaWxlLgoKIyMgUiBOb3RlYm9vawoKUiBOb3RlYm9vayBhIG5ldyB0eXBlIG9mIFIgTWFya2Rvd24gZm9yIGVtYmVkZGluZyBSIE1hcmtkb3duIGNvZGUgaW4gYW4gaHRtbCBmaWxlLiBUaGlzIGZpbGUgaXMgYW4gZXhhbXBsZSBvZiBhbiBSIE5vdGVib29rLgoKKiAiRmlsZSA+TmV3IEZpbGUgPlIgTm90ZWJvb2siIG1lbnUgdG8gY3JlYXRlIGEgbmV3IGZpbGUsIHdoaWNoIGluY2x1ZGVzICAKYG91dHB1dDogaHRtbF9ub3RlYm9va2AgIAppbiB0aGUgZmlsZSBoZWFkZXIuICAKCiogIlByZXZldyA+UHJldmlldyBOb3RlYm9vayIgYnV0dG9uIHRvIHNob3cgYSBwYWdlIHByZXZpZXcuCgoqIEV2ZXJ5dGltZSB5b3Ugc2F2ZSB5b3VyICJmaWxlLlJtZCIsIFIgU3R1ZGlvIGNyZWF0ZXMgImZpbGUubmIuaHRtbCIuCgojIyBSIGFzIGEgY2FsY3VsYXRvcgoKWW91IGNhbiBqdXN0IHR5cGUgaW4gbnVtYmVycyBhbmQgb3BlcmF0b3JzIGluIHRoZSBDb25zb2xlIHdpbmRvdy4gIAoKCmBgYHtyfQoxKzEKMiozCjQvNQpgYGAKCkVpdGhlciBgKipgIG9yIGBeYCBmb3IgcG93ZXIuCgpgYGB7cn0KMioqMwoyXjMKYGBgCgpNYW55IG1hdGggZnVuY3Rpb25zIGFuZCBzb21lIGNvbnN0YW50cyBhcmUgcHJlLWRlZmluZWQuCgpgYGB7cn0KZXhwKDEpCnNpbihwaS8yKQpgYGAKCiMjIFZhcmlhYmxlcwoKWW91IGNhbiB1c2UgYDwtYCBvciBgPWAgdG8gYXNzaWduIGEgdmFsdWUgdG8gYSB2YXJpYWJsZS4KCmBgYHtyfQp4IDwtIDEKeSA9IDIKeCt5CngveQp6ID0geCt5CnoKYGBgCgojIyBHZXR0aW5nIGluZm9ybWF0aW9uCgpgbHMoKWAgbGlzdHMgb2JqZWN0cyBpbiB5b3VyIHdvcmsgc3BhY2UuICAKYHJtKClgIHJlbW92ZXMgb2JqZWN0cy4gIAoKYGBge3IgaW5mb30KbHMoKQpybSh4LHkpCmxzKCkKYGBgCgpgZ2V0d2QoKWAgdG8gZ2V0IG9yIHNldCB0aGUgd29ya2luZyBkaXJlY3RvcnkKCmBgYHtyfQpnZXR3ZCgpCmBgYAoKYGhlbHAoKWAgb3IgYD9gIGZvciBoZWxwLiBgPz9gIGZvciByZWxhdGVkIHRvcGljcy4KCmBgYHtyfQo/c2luCj8/c2luCmBgYAoKWW91IGNhbiBhbHNvIHVzZSAiSGVscCIgcGFuZWwgb2YgUlN0dWRpby4KCiMjIFZlY3RvcnMKCmBjKCApYCBjYW5jYXRlbmF0ZXMgb2JqZWN0cyB0byBtYWtlIGEgdmVjdG9yLiAgCgpgYGB7ciB2ZWN0b3J9CmEgPSBjKDEsMiwzKQphCmBgYAoKWW91IGNhbiBhY2Nlc3MgZWFjaCBjb21wb25lbnQgYnkgYFsgXWAuICAKCmBgYHtyfQphWzFdCmFbYygyLDMpXQpgYGAKCk9wZXJhdG9ycyB3b3JrIGNvbXBvbmVudC13aXNlLiBUaGVyZSBhcmUgbWFueSBmdW5jdGlvbnMgdGhhdCB3b3JrIHdpdGggdmVjdG9ycy4KCmBgYHtyfQphK2EKMyphCmEqYQptZWFuKGEpCmV4cChhKQpgYGAKClJlZ3VsYXIgdmVjdG9ycyBjYW4gYmUgY3JlYXRlZCBieSBgOmAgb3IgYHNlcSgpYC4KCmBgYHtyfQpiID0gMTo1CmIKYyA9IHNlcSgwLDEwLDIpCmMKYGBgCgpZb3UgY2FuIGFsc28gY3JlYXRlIGFuIGVtcHR5IGFycmF5IGFuZCBsYXRlciBhc3NpZ24gdmFsdWVzLgoKYGBge3J9CmQgPSBhcnJheShkaW09MTApCmRbMV0gPSAxCmQKYGBgCgojIyBNYXRyaXgKCllvdSBjYW4gY3JlYXRlIGEgbmV3IG1hdHJpeCBieQoKYGBge3J9CkEgPSBtYXRyaXgoMCwgMiwgMykKQQpgYGAKCm9yIGJ5IGNvbWJpbmluZyB2ZWN0b3JzCgpgYGB7cn0KQiA9IHJiaW5kKDE6MywgNDo2KQpCCkMgPSBjYmluZCgxOjMsIDQ6NikKQwpgYGAKCllvdSBjYW4gcmVzaGFwaW5nIGEgdmVjdG9yIGludG8gYSBtYXRyaXguCgpgYGB7cn0KRCA9IDE6NgpkaW0oRCkgPSBjKDIsIDMpCkQKYGBgCgpZb3UgY2FuIGFjY2VzcyB0aGUgY29tcG9uZW50cyBieSBgWyAsIF1gCgpgYGB7cn0KRFsxLDJdCkRbMToyLDI6M10KRFsxLF0KRFssMV0KYGBgCgpNb3N0IG9wZXJhdG9ycyBhbmQgZnVuY3Rpb25zIHdvcmsgY29tcG9uZW50LXdpc2UuCgpgYGB7cn0KQSA9IHJiaW5kKDE6MiwgMzo0KQpBCkErQQpBKkEKc2luKEEpCmBgYAoKVG8gcGVyZm9ybSBtYXRyaXggbXVsdGlwbGljYXRpb24sIHdlIG11c3QgdXNlIGAlKiVgICAKCmBgYHtyfQpBICUqJSBBCmBgYAoKSW52ZXJzZSBpcyBnaXZlbiBieSBgc29sdmUoKWAgZnVuY3Rpb24KCmBgYHtyfQpBaW52ID0gc29sdmUoQSkKQWludgpBICUqJSBBaW52CmBgYAoKVXNlIGB0KClgIGZvciBhIHRyYW5zcG9zZSBtYXRyaXguCgpgYGB7cn0KdChBKQpgYGAKCiMjIExpc3QKCiogQSBgbGlzdGAgaXMgYSBjb2xsZWN0aW9uIG9mIGNvbXBvbmVudHMgb2YgYXJiaXRyYXJ5IHR5cGVzLCBlLmcuLCBudW1iZXJzLCBtYXRyaWNlcywgY2hhcmFjdGVycywgb3IgZnVuY3Rpb25zLiAgCgogICAgCmBgYHtyfQpMID0gbGlzdChuYW1lPSJPSVNUIiwgYWRkcmVzcz1jKCJUYW5jaGEiLCJPbm5hIiwiT2tpbmF3YSIpLCBhZ2U9NikKTFsxXQpgYGAKCiogYSBjb21wb25lbnQgY2FuIGJlIHNwZWNpZmllZCBieSBpdHMgaW5kZXggb3IgbmFtZQoKYGBge3J9CkxbWzFdXQpMW1sibmFtZSJdXQpMJG5hbWUKTCRhZGQKYGBgCgoKIyMgU2NyaXB0CgpZb3UgY2FuIHN0b3JlIG11bHRpcGxlIGxpbmVzIG9mIFIgY29tbWFuZHMgYXMgYSAqc2NyaXB0KiBmb3IgcmVwZWF0ZWQgdXNlLgoKKiBDaG9vc2UgIkZpbGUgPk5ldyBGaWxlID5SIFNjcmlwdCIgdG8gbWFrZSBhIG5ldyBmaWxlLiBUeXBlIGluLCBmb3IgZXhhbXBsZToKCmBgYAojIGJlcm4uUjogQmVybm91bGxpIGRpc3RyaWJ1dGlvbiBvZiAwIGFuZCAxCnAgPSAwLjUgICMgcHJvYmFiaWxpdHkgb2YgMQppZiAocnVuaWYoMSkgPCBwKXsgICMgc2FtcGxlIGEgcmFuZG9tIG51bWJlciB1bmlmb3JtbHkgaW4gMCBhbmQgMQogICAgeCA9IDEgICMgMSBpZiB0aGF0IGZhbGxzIGJlbG93IHAKICB9ZWxzZXsgCiAgICB4ID0gMCAgIyAwIG90aGVyd2lzZQogIH0KcHJpbnQoeCkKYGBgCgpcIyBzdGFydHMgYSBjb21tZW50LCB3aGljaCBpcyBza2lwcGVkIGJ5IFIgYnV0IGhlbHBmdWwgZm9yIGh1bWFucywgaW5jbHVkaW5nIHlvdXJzZWxmLCB0byB1bmRlcnN0YW5kIHdoYXQgdGhlIGNvZGUgZG9lcyBhbmQgd2hhdCB0aGUgcGFyYW1ldGVyIG1lYW5zLgoKKiBDaG9vc2UgIkZpbGUgPlNhdmUgQXMuLi4iIGFuZCBuYW1lIGl0ICdiZXJuLlInCgpHbyB0byB0aGUgZGlyZWN0b3J5IHdoZXJlIHRoZSBmaWxlIHdhcyBzYXZlZCwgc3RhcnQgUiwgYW5kIHR5cGUgaW46CgpgYGB7ciBzY3JpcHR9CnNvdXJjZSgiYmVybi5SIikKYGBgCgojIyBEZWZpbmluZyB5b3VyIG93biBmdW5jdGlvbgoKRm9yIHJlcGVhdGluZyBhIHByb2NlZHVyZSB3aXRoIGRpZmZlcmVudCBwYXJhbWV0ZXJzLCB5b3Ugd291bGQgZGVmaW5lIHlvdXIgb3duIGZ1bmN0aW9uLCBpbiB0aGUgZm9sbG93aW5nIGNvbnZlbnRpb246CgpgYGB7ciBmdW5jdGlvbn0KYmVybiA8LSBmdW5jdGlvbihwKXsgICMgcDogcHJvYmFiaWxpdHkgb2YgMQogIGlmIChydW5pZigxKSA8IHApeyAgIyBzYW1wbGUgYSB1bmlmb3JtIHJhbmRvbSBudW1iZXIKICAgIHggPSAxICAgICMgMSBpZiB0aGF0IGZhbGxzIGJlbG93IHAKICB9ZWxzZXsKICAgIHggPSAwICAgICMgMCBvdGhlcndpc2UKICB9CiAgcmV0dXJuKHgpICAKfQpgYGAKCllvdSBjYW4gY2FsbCBpdCB3aXRoIGFueSBwYXJhbWV0ZXI6CgpgYGB7ciBiZXJufQpiZXJuKDAuMikKYmVybigwLjgpCmBgYAoKIyMgTG9vcHMKClRyeSBtYWtpbmcgYmlvbm9taWFsIGRpc3RyaWJ1dGlvbiBieSBhZGRpbmcgMXMgb2YgbiBiZXJub3VsbGkgc2FtcGxlcyB1c2luZyBmb3IgbG9vcDoKCmBgYHtyIGxvb3B9CmJpbm9tID0gZnVuY3Rpb24obiwgcCkgeyAjIG46c2FtcGxlcywgcDpwcm9iYWJpbGl0eQogIHkgPSAwICAgIyBwcmVwYXJlIGEgY291bnRlcgogIGZvciAoaSBpbiAxOm4pIHsKICAgIHggPSBiZXJuKHApICAgIyAxIG9yIDAKICAgICMgcHJpbnQoeCkgICAgIyB1bmNvbW1lbnQgdG8gY2hlY2sgZWFjaCB4CiAgICB5ID0geSArIHgKICB9CiAgcmV0dXJuKHkpCn0KYGBgCgpZb3UgY2FuIHRlc3QgaXQgd2l0aCBkaWZmZXJlbnQgcGFyYW1ldGVyczoKCmBgYHtyIGJpbm9tfQpiaW5vbSg1LCAwLjUpCmJpbm9tKDEwLCAwLjMpCmBgYAoKIyMgUGxvdHRpbmcKCkxldCB1cyB0YWtlIHNhbXBsZXMgZnJvbSBhIGJpbm9taWFsIGRpc3RyaWJ1dGlvbgoKYGBge3IgYmlub21faGlzdH0KbSA9IDEwMDAgICMgbnVtYmVyIG9mIHNhbXBsZXMKbiA9IDEwCnAgPSAwLjUKeCA9IGFycmF5KGRpbT1tKSAgICMgcHJlcGFyZSBhbiBlbXB0eSBhcnJheQpmb3IgKGkgaW4gMTptKSB7CiAgeFtpXSA9IGJpbm9tKG4sIHApICAjIGJldHdlZW4gMCB0byBuCn0KYGBgCgpZb3UgY2FuIHZpc3VhbGl6ZSB0aGUgcmF3IGRhdGEgYnkgYHBsb3QoKWAgZnVuY3Rpb24uCgpgYGB7cn0KcGxvdCh4KQpgYGAKCm9yIG1ha2UgYSBoaXN0b2dyYW0gYnkgYGhpc3QoKWAuCgpgYGB7cn0KaGlzdCh4KQpgYGAKCllvdSBjYW4gcGxvdCBhbnkgZnVuY3Rpb24gYnkgc3BlY2lmeWluZyB4IGFuZCB5IGZvciBgcGxvdCgpYCBmdW5jdGlvbi4KCmBgYHtyfQpYID0gc2VxKDAsIG4pCnBsb3QoWCwgZGJpbm9tKFgsIG4sIHApKQpgYGAKClRvIG92ZXJsYXkgYSBjdXJ2ZSBvbiBhbiBleGlzdGluZyBwbG90LCB1c2UgYGxpbmVzKClgIGZ1bmN0aW9uLgoKYGBge3J9Cmhpc3QoeCwgZnJlcT1GQUxTRSkgICMgbm9ybWFsaXplIHRvIHN1bSB1cCB0byAxCmxpbmVzKFgsIGRiaW5vbShYLCBuLCBwKSkKYGBgCgojIyAzRCBQbG90cwoKUiBoYXMgYSBzaW1wbGUgM0QgcGxvdCBmdW5jdGlvbiBgcGVyc3AoKWAuCgpgYGB7cn0KeCA9IHNlcSgtNSwgNSwgMC4xKQp5ID0gc2VxKDAuNSwgMywgMC4xKQpmID0gZnVuY3Rpb24oeCx5KXtleHAoLXgqKjIvKDIqeSoqMikpL3l9CnogPSBvdXRlcih4LCB5LCBmKQpwZXJzcCh4LCB5LCB6LCB0aGV0YT00NSwgcGhpPTMwKQpgYGAKClRoZXJlIGFyZSBmYW5jaWVyIDNEIHBhY2thZ2VzIGxpa2UgYHBsb3QzRGAgYW5kIGByZ2xgLgoKIyMgVG8gYmUgY292ZXJlZCBpbiBsYXRlciBjaGFwdGVycwoKKiBwcm9iYWJpbGl0eSBkaXN0cmlidXRpb25zCiogZGF0YSBmcmFtZQoqIHByZS1sb2FkZWQgZGF0YSBzZXRzCiogcmVhZGluZyBkYXRhIGZyb20gYSBmaWxlCg==

1. Introduction