rasmusab

---
title: "A Bayesian Model to Calculate Whether My Wife is Pregnant or Not"
author: "Rasmus Bååth"
output: html_document
---

```{r echo = FALSE}
library(knitr)
options(digits=2)
```

rasmusab

# A Bayesian model that calculates a probability that a couple is fertile 
# and pregnant. Please use this for fun only, not for any serious purpose 
# like *actually* trying to figure out whether you are pregnant.

# Enter your own period onsets here:
period_onset <- as.Date(c("2014-07-02", "2014-08-02", "2014-08-29", "2014-09-25", 
                          "2014-10-24", "2014-11-20", "2014-12-22", "2015-01-19"))
# If you have no dates you can just set days_between_periods to c() instead like:
# days_between_periods <- c()
days_between_periods <- as.numeric(diff(period_onset))

rasmusab

0   CLS
10  REM SETTING UP THE VARIABLES AND DATA
20 
30  REM NO OF PUPS PER DEN FROM A SAMPLE OF 16 WOLF DENS
40  DATA 5, 8, 7, 5, 3, 4, 3, 9, 5, 8, 5, 6, 5, 6, 4, 7 
50 
60  NDATA = 16
70  NSAMPLES = 1000
75  NBURNIN = 250
80  NPARAMS = 2

rasmusab

---
title: "Easy Bayesian Bootstrap in R"
author: "Rasmus"
date: "04/06/2015"
output:
  html_document: 
    keep_md: yes
    self_contained: no
---

rasmusab

# From https://aschinchon.wordpress.com/2015/06/23/the-2d-harmonograph-in-shiny/
# All code by Antonio S. Chinchón (@aschinchon), just rearranged so that you can directly
# run it by copy and pasting it into an R console.

library(shiny)
CreateDS = function () {
  f=jitter(sample(c(2,3),4, replace = TRUE))
  d=runif(4,0,1e-02)
  p=runif(4,0,pi)
  xt = function(t) exp(-d[1]*t)*sin(t*f[1]+p[1])+exp(-d[2]*t)*sin(t*f[2]+p[2])

rasmusab

library(ggplot2)
library(plyr)

n_games <- readRDS("n_games.Rdata")
n_positions <- readRDS("n_positions.Rdata")
fullmoves_white <- readRDS("fullmoves_white.Rdata")

fits <- readRDS("chess_fits.Rdata")
fit <- fits$fit
fit <- fit 

rasmusab

# Takes the matrix created by this script: https://gist.github.com/rasmusab/fb98cced046d4c675d74
# and calculates some statistics and fits some models to it. All this is pretty memory heavy so saving
# the interesting parts using saveRDS so the script only need to be run once.

set.seed(123)
games <- as.data.frame(readRDS("milionbase_matrix.rds"))

# Saving some info
n_games <- max(games[,"game_id"])
fullmoves_white <- games[ games[, "fullmoves"] < 100 & games[,"active_player"] == 1, "fullmoves"]

rasmusab

# Takes a json file describing chess games and produces a matrix with one row
# per turn showing how many pieces are left, one column per piece.

### Don't run this in R studio because it will take up twice the RAM 
### as R will make copies instead of references.

library(jsonlite)
library(stringi)

# Path to your json file as produced by this script: https://gist.github.com/rasmusab/07f1823cb4bd0bc7352d

rasmusab

# Converts the millionbase chess PGN database (http://www.top-5000.nl/pgn.htm) to json
# with one json dictionary per row. (That is, the resulting file is contain multiple json objects, 
# not just one large).

import json
import chess.pgn # From python-chess https://github.com/niklasf/python-chess

pgn = open("millionbase-2.22.pgn") # Or where you have put it
fout = open('milionbase-2.22.json', 'w') # Or where you want it

rasmusab

library(microbenchmark)
library(arm)
library(rstan)
library(bbmle)

log_post <- function(par, y, x) {
  sigma <- exp(par[1])
  intercept <- par[2]
  beta <- as.matrix(par[-c(1,2)])
  mu <- intercept + x %*% beta