Why statistics?

From Motulsky (2014) Chapter 1:

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).

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==