cboettig

##########################################################
# show updating process iteratively
# just paste all of this code into R to execute it

# generate data
n <- 20
y <- rnorm( n , mean=7 , sd=0.5 )
# assign prior and compute posterior as we add each y value to observations
prior.mu <- 3
k.sigma<-sd(y)

cboettig

# @file: abc.R
# @author: Carl Boettiger <cboettig@gmail.com>
# @date: 11 April, 2011
# Description: A trivial example of Approximate Bayesian Computing
# We will simulate under a Gaussian process to determine 

# The target data:
TrueMean <- 7     ## irrelevant note: "True" is rather a frequentist name for this
TrueSD <- 1
Npts <- 100

cboettig

## The observed data. 
X <- rnorm(1000, 2, 5)

loglik <- function(mu, sigma, X){
  ## The likelihood function
  sum( dnorm(X, mean=mu, sd=sigma, log=TRUE) )
}

## initial starting point 
pars <- c(mu=5, sigma=15)

cboettig

beta <- function(i, Delta_T=1){
  1/(1+Delta_T*(i-1))
}

step_fn <- function(pars, stepsizes = .02){
# Sequential random updating 
  j <- sample(1:length(pars), 1)
  pars[j] <- rnorm(1, pars[j], stepsizes)
  pars
}

cboettig

# image trees are named something like paup1.tre to paup100.tre
require(ape)
# could use a for loop, but sapply is faster and will parallelize
sapply(1:100, function(i){
  tree <- read.tree(paste("paup", i, ".tre", sep=""))
  pdf(file=paste("paup", i, ".pdf", sep="")) 
  plot(tree) # add fancy plotting options etc,
  dev.off()
})

cboettig

## this function f  will remember whatever options it is passed
f <- function(y, ...){
  out <- y + dnorm(...)
  list(out, opts=list(...))
}
## This is the way to call the function using its exact arguments
a <- f(1, x=4, sd=2, mean=4) 
b <- do.call(f, c(list(y=1), a$opt)) 
identical(a,b) # TRUE

cboettig

#####  socialR header info ##### 
require(socialR)
script <- "example.R" # Must specify the script name! 
gitcommit(script)     # Must commmit at start!

##### Optional but a good idea ##### 
#almost requisite, if not me running "repo"/"dir"
gitopts = list(user = "cboettig", dir = "demo", repo = "socialR") 
on.exit(system("git push")) #  For git links.  May prompt for pw,
tags <- "test"  ## multiple possible: space, delim, multiple items, etc.  

cboettig

data(faithful)
attach(faithful)

 ## EM Algorithm

W = waiting
s = c(0.5, 40, 90, 16, 16)

em = function(W,s) {
  Ep = s[1]*dnorm(W, s[2], sqrt(s[4]))/(s[1]*dnorm(W, s[2], sqrt(s[4])) +

cboettig

data(faithful)
attach(faithful) 
#Assumes two gaussian populations of eruptions
W = waiting
# Guess of the parameters: 
s = c(p=0.5, mu1=50, mu2=90, sigma1=30, sigma2=30) 
 
 
# EXPECTATION STEP
expectation_step <- function(observed, parameters){

cboettig

#segue_test.R
on.exit(stopCluster(myCluster))

library(segue)
setCredentials(getOption("amazon_key"), getOption("amazon_secret"))


# Jeffery Breen's example, adapted with "custom function" works fine
my_fun <- mean # a custom definition
myList <- NULL