# R Stuff, Chapter 1

1. Assignment .

The real way of doing assignment is with a <- operator.
x <- 4
or
4 -> x

Shortcut
x=4
(note = is shortcut for <-, and I use it always)

2. First graph (probability mass)

x=c(0,1,2) #assigns a vector; "c" is for concatenate
f=c(.09,.42,.49)

plot(x,f,typ='h')

#other settings
plot(x,f,typ='l') #lines
plot(x,f,typ='p') #points
plot(x,f,typ='n') #set up axes, nothing else

#ok first graph
plot(x,f,typ='h')
points(x,f)

#fancy graph
par(mar=c(4,4,1,1),mgp=c(2.5,1,0),las=1,cex=1.2) #see help(par)
xlim=c(-.5,2.5),ylim=c(0,.5),axes=F)
axis(1,at=c(0,1,2),label=c(0,1,2))
axis(2)
box()
points(x,f,pch=21,bg='red',cex=1.1)

3. Binomial Distribution
y=0:20 #assigns x to 0,1,2,..,20
f=dbinom(y,20,.7) #probability mass function
plot(y,f,type='h')
points(y,f)

4. Expected Value
Operations in R are "element-wise" by default
Here is the expected value code in R.
x=1:5
f=c(.05,.15,.25,.25,.2)
terms=x*f
ev=sum(terms)

5. Realizations
rbinom(1,1,.7) #one flip of toast
rbinom(1, 20, .7) #one expt, #successes in 20 flips
rbinom(5, 20, .7) #five experiments of 20 flips
y=rbinom(200,20,.7) #200 experiments of 20 flips
hist(y,breaks=seq(-.5, 20.5, 1),freq=T)

freqs=table(y) #frequencies of realization values
props=freqs/200 # proportions of realization values
plot(props, xlab='Value', ylab='Relative Frequency',
typ='h',xlim=c(0,20))
points(0:20,dbinom(0:20,20,.7),pch=21,bg='orchid') #nor great covergence

#redo for 20,000 to show convergence
y=rbinom(20000,20,.7) #200 experiments of 20 flips
freqs=table(y) #frequencies of realization values
props=freqs/20000 # proportions of realization values
plot(props, xlab='Value', ylab='Relative Frequency',
typ='h',xlim=c(0,20))
points(0:20,dbinom(0:20,20,.7),pch=21,bg='orchid') #nor great covergence

#look at p-hat
p.hat=y/20
freqs=table(p.hat) #frequencies of realization values
props=freqs/20000 # proportions of realization values
plot(props, xlab='Value', ylab='Relative Frequency',
typ='h',xlim=c(0,1))

Estimators of p
y=rbinom(10000,20,.7)
p.hat=y/20
mean(p.hat) #sample mean
var(p.hat) #sample variance (N-1 in denominator)
sd(p.hat) #sample std. deviation (N-1 in denominator), std error

#define function
bsms=function(m,n,p)
{
z=rbinom(m,n,p)
mean(z)
}
#call function
bsms(10,20,.7)

M=10000
bsms.realization=1:M #define the vector ppes.realization
for(m in 1:M) bsms.realization[m]=bsms(10,20,.7)

bsms.props=table(bsms.realization)/M
plot(bsms.props, xlab="Estimate of Expected Value (Sample Mean)",
ylab="Relative Frequency", type='h')