Advertisement
Guest User

Untitled

a guest
Apr 23rd, 2017
63
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.19 KB | None | 0 0
  1. # chapter 2 code
  2.  
  3. library(rethinking)
  4.  
  5. col1 <- "black"
  6. col2 <- col.alpha( "black" , 0.7 )
  7.  
  8. #####
  9. # introducing the golem
  10. # show posterior as data comes in
  11.  
  12. # 1 indicates 'water'; 0 indicates 'land'
  13. d <- c(1,0,1,1,1,0,1,0,1) # length 9
  14. l <- ifelse( d==1 , "W" , "L" )
  15.  
  16. # prior uniform - "Laplace" prior
  17. a <- 1
  18. b <- 1
  19.  
  20. # draw first plot
  21. i <- 0
  22. curve( dbeta(x,a,b) , from=0 , to=1 , xlab="proportion water" , ylab="" , col=col1 , lwd=2 , ylim=c(0,2.75) , xaxt="n" , yaxt="n" )
  23. axis( 1 , at=c(0,0.5,1) , labels=c(0,0.5,1) )
  24. title( ylab="plausibility" , mgp=c(1,1,0) )
  25. text( 0.15 , 2.5 , paste("n =",i) )
  26. mtext( paste(l) , at=seq(from=0,to=1,length.out=length(l)) , col=c( rep( "darkgray" ,length(l)) ) )
  27.  
  28. # execute this code once for each subsequent plot
  29. i <- i + 1
  30. prior.a <- a
  31. prior.b <- b
  32. a <- a + d[i]
  33. b <- b + 1 - d[i]
  34. curve( dbeta(x,a,b) , from=0 , to=1 , xlab="probability of water" , ylab="" , col=col1 , lwd=2 , ylim=c(0,2.75) , xaxt="n" , yaxt="n" )
  35. axis( 1 , at=c(0,0.5,1) , labels=c(0,0.5,1) )
  36. text( 0.15 , 2.5 , paste("n =",i) )
  37. curve( dbeta(x,prior.a,prior.b) , add=TRUE , col=col1 , lty=2 )
  38. mtext( paste(l) , at=seq(from=0,to=1,length.out=length(l)) , col=c( rep(col1,i) , rep( "darkgray" ,length(l)-i) ) )
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement