BayesianBlocks.Rmd

```{r}
library("repr")
options(warn=-1)# repr.plot.width=14, repr.plot.height=8) 
```

Algorithm

```{r}
bayesian_blocks <- function(xs, prior=7.61, nn_vec=rep(1,length(xs))){
    
    
    xs <- sort(xs)
    begin <- xs[1]
    end <- xs[length(xs)]
    
    xs <- xs[!(nn_vec==0)]
    nn_vec <- nn_vec[!(nn_vec==0)]
    
    N <- length(xs)
    start <- xs[1]
    stop <- xs[N]
  
    blocks <- c(start, 0.5*(xs[2:N]+xs[1:N-1]), stop);
    block_length <- stop - blocks
    # ---------------------------------------------
    # Start with first data cell; add one cell at
    # each iteration
    # ---------------------------------------------

    best <- NULL 
    last <- NULL
    supp <- NULL
    
    for (R in 1:N){
    # Compute fit_vec : fitness of putative last block (end at R)
        arg_log <- block_length[1:R] - block_length[R+1]
        #arg_log[arg_log <= 0] <- Inf

        nn_cum_vec <- cumsum(nn_vec[R:1])
        nn_cum_vec <- nn_cum_vec[R:1]
        
        fit_vec <- nn_cum_vec * (log(nn_cum_vec) - log(arg_log))
        fit_vec <- fit_vec - prior
        
        supp <- c(0, best) + fit_vec
        
        best <- c(best, max(supp))
        last <- c(last, which.max(supp))
         
}    
    
    #---------------------------------------------
    # Now find changepoints by iteratively peeling off the last block
    #---------------------------------------------
    index <- last[N]

    change_points <- NULL

    while(index > 1){
        change_points <- c(index, change_points)
        index <- last[index - 1]
    }
    
#     change_points <- c(change_points, blocks[1])
    
    return(c(begin, blocks[change_points], end))
}
```

Test dataset

```{r}
test <- c(rnorm(1000, 0, 1), rgamma(1000, 9, 2))
test <- test[(test > -5) & (test < 10)]
test <- sort(test)
```

```{r}
cp_h <- bayesian_blocks(test, prior=7.61)
cp_h
```

```{r}
options(repr.plot.width=6, repr.plot.height=4) 
#par(mfrow=c(1,2))

h_t  <- hist(test, breaks=200, freq=FALSE, col=rgb(0,0,1,0))

h_bb <- hist(test, breaks=cp_h, freq=FALSE, col=rgb(1,0,0,0), add=TRUE)

```

Test priors

```{r}
ncp_prior <- function(p0, N){ 4 - log(73.53 * p0 * N^(-0.478))}

p0 <- c(0.0001, 0.001, 0.01, 0.1)
N  <- length(test)

par(mfrow=c(1,2))

for (i in 1:length(p0)){
    cp <- bayesian_blocks(xs = test, prior=ncp_prior(p0[i], N))
    h  <- hist(test, breaks=cp, freq=FALSE, plot=FALSE)
    
    plot(h_t)
    plot(h, add=TRUE)
}

```

Import "real" dataset 

```{r}
data   <- read.table("./Data/B19036_AmCsCo_20180316.dat", skip=2)
counts <- data[,1]
bins   <- c(0:8191)

# Preliminary computation: start data at the first non-empty bin

idx    <- which(counts >0)[1]
counts <- counts[idx:length(counts)]
bins   <- bins[idx:length(bins)]
```

```{r}
options(repr.plot.width=6, repr.plot.height=4)  #to set graph size

plot(bins, counts, type='s')#, xlim=c(6520,6545))

```

Prior

```{r}
ncp_prior <- c(1, 1.3, 2, 2.5, 3.2, 4)
#ncp_prior <- 7.61

```

Algorithm to rebin

```{r}
#create histogram

rebin_bb <- function(bins, counts, change_points){
    
    rebin <- NULL
    y <- 0
    dn <- bins[2]-bins[1]
    
    for (i in 1:(length(change_points)-1)){
        ifelse( ( bins < change_points[i+1] ) & ( bins > change_points[i] ),
               { y <- sum(counts[bins < change_points[i+1] & bins > change_points[i]])
                 y <- y/(dn*length(counts[bins < change_points[i+1] & bins > change_points[i]]))
                 rebin <- c(rebin, y)},
               next
            )
    }
    return(c(0,rebin,0))
}
```

```{r}
cp <- bayesian_blocks(xs = bins, prior=1, nn_vec = counts)
rebin <- rebin_bb(bins = bins, counts = counts, change_points = cp)
```

```{r}
plot(bins, counts, col= 'grey', 
     type='s', lwd=0.1, log='y')
lines(cp, rebin, type='s', col='red')
#abline(v=cp, col='red')


```

```{r}

```