## R code 0.1
print( "All models are wrong, but some are useful." )
## [1] "All models are wrong, but some are useful."
## R code 0.2
( log( 0.01^200 ) )
## [1] -Inf
( 200 * log(0.01) )
## [1] -921.034
## R code 0.3
# Load the data:
# car braking distances in feet paired with speeds in km/h
# see ?cars for details
data(cars)
# fit a linear regression of distance on speed
m <- lm( dist ~ speed , data=cars )
# estimated coefficients from the model
coef(m)
## (Intercept) speed
## -17.579095 3.932409
# 95% confidence intervals of the estimates
confint(m)
## 2.5 % 97.5 %
## (Intercept) -31.167850 -3.990340
## speed 3.096964 4.767853
# plot residuals against speed
plot( resid(m) ~ speed , data=cars )
## R code 0.4
#install.packages(c("coda","mvtnorm","devtools"))
#library(devtools)
#devtools::install_github("rmcelreath/rethinking")
## R code 2.1
ways <- c( 0 , 3 , 8 , 9 , 0 )
ways/sum(ways)
## [1] 0.00 0.15 0.40 0.45 0.00
## R code 2.2
dbinom( 6 , size=9 , prob=0.5 )
## [1] 0.1640625
## R code 2.3
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- rep( 1 , 20 )
# compute likelihood at each value in grid
likelihood <- dbinom( 6 , size=9 , prob=p_grid )
# compute product of likelihood and prior
unstd.posterior <- likelihood * prior
# standardize the posterior, so it sums to 1
posterior <- unstd.posterior / sum(unstd.posterior)
## R code 2.4
plot( p_grid , posterior , type="b" ,
xlab="probability of water" , ylab="posterior probability" )
mtext( "20 points" )
## R code 2.5
prior <- ifelse( p_grid < 0.5 , 0 , 1 )
prior <- exp( -5*abs( p_grid - 0.5 ) )
## R code 2.6
library(rethinking)
## Loading required package: rstan
## Loading required package: Rcpp
## Loading required package: ggplot2
## rstan (Version 2.8.0, packaged: 2015-09-19 14:48:38 UTC, GitRev: 05c3d0058b6a)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
## Loading required package: parallel
## rethinking (Version 1.57)
globe.qa <- map(
alist(
w ~ dbinom(9,p) , # binomial likelihood
p ~ dunif(0,1) # uniform prior
) ,
data=list(w=6) )
# display summary of quadratic approximation
precis( globe.qa )
## Mean StdDev 5.5% 94.5%
## p 0.67 0.16 0.42 0.92
## R code 2.7
# analytical calculation
w <- 6
n <- 9
curve( dbeta( x , w+1 , n-w+1 ) , from=0 , to=1 )
# quadratic approximation
curve( dnorm( x , 0.67 , 0.16 ) , lty=2 , add=TRUE )
