Nothing
##############################################################
#' Calculate Negative-binomial on vector of parameters
#'
#' Computes the negative-binomial posterior predictive density
#' from input parameter vectors
#' corresponding to each possible run length for the current
#' time point. Outputs a vector of probabilities for use
#' in the accompanying poisson functions.
#'
#' @param x the current data point
#'
#' @param a matrix of alpha params
#'
#' @param b matrix of beta params
#'
#' @return Matrix of negative binomial pdf values corresponding
#' to each possible run length, for use in accompanying poisson
#' probability functions.
#'
#' @import stats
#'
#' @docType methods
#'
#' @export
##############################################################
negbinpdf <-function(x, a, b) {
nbpdfmat<- matrix(nrow=nrow(a), ncol = ncol(a))
for(xi in 1:nrow(a)){
for(dimi in 1:ncol(a)){
nbpdfmat[xi, dimi]<-dnbinom(x=x[dimi], size=a[xi, dimi],
prob=b[xi, dimi], log = TRUE)
}
}
return(nbpdfmat)
}
##############################################################
#' Initialize vectors for poisson probability functions
#'
#' Takes in the desired initialization parameters,
#' initializes the vectors needed for the poisson probability
#' function \code{poisson_update}
#'
#' @param init_params The list of parameters to be used for
#' initialization
#' @param dims the dimensionality of the data
#'
#' @return List of vectors to be used in the iteratively
#' updating algorithm of parameters describing the
#' underlying gaussian distribution of the data.
#'
#' @docType methods
##############################################################
poisson_init<- function(init_params = list(a=1, b=1), dims){
beta0 <- matrix(rep(init_params$b,times=dims),ncol=dims)
alpha0 <- matrix(rep(init_params$a, times=dims), ncol=dims)
return(list(alphaT=alpha0, betaT=beta0, muT=alpha0/beta0))
}
##############################################################
#' Update the poisson parameters
#'
#' Updates the parameters of the poissons based on each possible
#' run length, after taking into consideration the most recent
#' data point
#'
#' @param update_params0 The initialization parameters,
#' corresponding to predicting a changepoint (run length=0)
#' @param update_paramsT The vectors of parameters corresponding
#' to each possible run length, updated with each incoming data point
#' @param datapt the current data point
#' @param Rlength the length of the current vector of possible run
#' lengths
#' @param skippt If the current point should be skipped in the updating
#' because it was missing, and missPts was set to skip
#'
#' @return The list of the parameters for gaussians corresponding
#' to each possible runlength up to the current data point.
#' Lengths of vectors should correspond the length of the
#' R vector ("run length vector")
#'
#' @docType methods
##############################################################
poisson_update<- function(datapt, update_params0,
update_paramsT, Rlength, skippt = FALSE){
dims<- length(datapt)
alphaT<- update_paramsT$alphaT[1:Rlength,]
betaT<-update_paramsT$betaT[1:Rlength,]
muT<-update_paramsT$muT[1:Rlength,]
currpt<- datapt
currpt<-as.vector(datapt)
betaT2<-as.matrix(betaT) # if the pt was missing, dont update the count by adding 1
if(!skippt) betaT2<-betaT2+matrix(1, ncol=dims, nrow=Rlength) # add 1 to all betas
if(ncol(as.matrix(betaT2))<ncol(as.matrix(update_params0$betaT))) betaT2<- t(as.matrix(betaT2))
betaT <- rbind(update_params0$betaT,betaT2) # append beta0
alphaT2<-t(t(alphaT)) # dont add the point if it should be skipped
if(!skippt) alphaT2<-t(t(alphaT)+datapt) # add x to alpha
if(ncol(alphaT2)<ncol(update_params0$alphaT)) alphaT2<- t(alphaT2)
alphaT <- rbind(update_params0$alphaT,alphaT2) # append alpha0
muT<-alphaT/betaT # muT is not needed by the model, just for keeping track of this param
return(list(alphaT=alphaT, betaT=betaT, muT=muT))
}
##############################################################
#' Compute predictive probabilities based on Poisson
#'
#' Compute the probability of observing the current point,
#' given the current parameters of the poisson for each
#' possible run length. Returns a vector of predictive
#' probabilities from each possible run length, the
#' parameters of the poisson, the most likely lambda of the
#' current poisson, and the current point.
#'
#' @param update_params0 The initialization parameters,
#' corresponding to predicting a changepoint (run length=0)
#' @param update_paramsT The vectors of parameters corresponding
#' to each possible run length, updated with each incoming data point
#' @param datapt the current data point
#' @param Rlength the length of the current vector of possible run
#' lengths
#' @param missPts the method set to handle missing points
#' @param time the number of time points passed so far
#' @param cps the current most likely list of changepoints
#' @param skippt If the current point should be skipped in the updating
#' because it was missing, and missPts was set to skip
#'
#' @return Returns a vector of predictive
#' probabilities from each possible run length, the
#' parameters of the gaussian, the most likely mean of the
#' current gaussian, and the current point.
#'
#' @docType methods
##############################################################
poissonProb<- function(update_params0, update_paramsT, datapt,
time, cps, missPts, Rlength, skippt=FALSE){
prevRlength<- nrow(as.matrix(update_paramsT$alphaT))
alphaT<- as.matrix(update_paramsT$alphaT[1:(Rlength),])
if(ncol(alphaT)<length(datapt)) alphaT<- t(alphaT)
betaT<-as.matrix(update_paramsT$betaT[1:(Rlength),])
if(ncol(betaT)<length(datapt)) betaT<- t(betaT)
muT<-as.matrix(update_paramsT$muT[1:(Rlength),])
if(ncol(muT)<length(datapt)) muT<- t(muT)
currpt<- datapt
# get the most likely mu for the current point
prevmut<- time - cps[length(cps)]+1 # the length of the most likely current run = time - most recent cp
# prevmut should stay the same because its in reference to the start of the list, the end is the part thats truncated
if(prevmut<2) prevmut<-2
if(time<2) prevmut<-1
if(prevmut>(Rlength)){
prevmut<- length(Rlength) # this could happen due to truncation, previous cp could be past the trunc lim
}
currpt2<- as.vector(currpt)
if(currpt2<0) stop("poisson distribution is counts data, must be positive numbers, current point is", currpt2)
# make sure its integers
currpt2<- as.integer(currpt2)
# the truncated sized param vectors used in student pdf
predProbs <- negbinpdf(currpt2,alphaT,betaT/(betaT+1)) # step 3: predictive probability
if(any(!is.finite(predProbs))) stop("poission predictive probability not finite, check input is integers")
# the vectors that will be used in the next round should be at most one longer since the run length could grow by at most one
update_paramsT<-poisson_update(datapt=currpt, update_params0, update_paramsT, (Rlength),
skippt=skippt)
tryCatch({currmu<-as.matrix(update_paramsT$muT)[prevmut,]}, error = function(error_condition) {
print("trying to update currmu")
print(paste("dims muT: ", dim(as.matrix(update_paramsT$muT))))
print(paste("prevmut: ", prevmut))
print(paste(error_condition, "\n"))
return(list(update_paramsT=update_paramsT, predProb =predProbs, currpt= currpt, currmu= currmu ))})
return(list(update_paramsT=update_paramsT, predProb =predProbs, currpt= currpt, currmu= currmu )) # going to have tp list and unlist the current mu
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.