# location-scale t-distrubution functions
# not provided in the standard t-distribution in R
# mu is the location, a is the scale
dt_ls <- function(x, df, mu, a) 1/a * dt((x - mu)/a, df)
pt_ls <- function(x, df, mu, a) pt((x - mu)/a, df)
qt_ls <- function(prob, df, mu, a) qt(prob, df)*a + mu
rt_ls <- function(n, df, mu, a) rt(n,df)*a + mu
studentpdf <- function(x, mu, var, nu) {
	c <- exp(log(gamma(nu/2 + 0.5)) - log(gamma(nu/2))) * (nu * pi * var)^(-0.5)
	p <- c * (1 + (1 / (nu*var)) * (x-mu)^2)^(-(nu+1/2))
}

hazard_func <- function(r, lam) {
	odds <- 1/lam * rep(1, r)
	return (1/lam)
}

obc <- setRefClass("obc",
	fields = list(
		mu 		= "numeric",
		lambda 	= "numeric",
		alpha 	= "numeric",
		beta 	= "numeric",
		T 		= "numeric",
		R 		= "data.frame",
		maxes 	= "numeric",
		ct 		= "numeric"
	),
	methods = list(
		initialize = function(mu0=1, lambda0=lambda, alpha0=1, beta0=1) {
			mu <<- mu0
			lambda <<- lambda0
			alpha <<- alpha0
			beta <<- beta0
			ct <<- 0

			T <<- 1000
			R <<- as.data.frame(matrix(0, nrow=T+1, ncol=T+1))
			R[1,1] <<- 1
			maxes <<- numeric(T+1)
		},
		calculate = function( x ) {
			ct <<- ct + 1
			pred_post <- dt_ls(x, 2*alpha, mu, beta*(lambda+1)/(alpha*lambda))
			H <- hazard_func(ct, lambda)
			R[2:(ct+1), ct+1] <<- R[1:ct, ct] * pred_post * (1-H)
			R[1, ct+1] <<- sum( R[1:ct, ct] * pred_post * H )
		  	R[,ct+1] <<- R[,ct+1] / sum(R[,ct+1])

		  	mu_new 		<- c(mu[1], (lambda*mu + x) / (lambda + 1))
		  	lambda_new 	<- c(lambda[1], lambda + 1)
		  	alpha_new 	<- c(alpha[1], alpha + 0.5)
			beta_new 	<- c(beta[1], beta + (lambda *(x-mu)^2) / (2*(lambda+1)))

			mu <<- mu_new
			lambda <<- lambda_new
			alpha <<- alpha_new
			beta <<- beta_new
			maxes[ct] <<- which(R[,ct] == max(R[,ct]))
		}
	)
)

# data <- c(rnorm(333, 5, .1), rnorm(333, 2, .1), rnorm(334, 8, .1))
# o <- obc$new(mu0 = 5, lambda0=10, alpha0=10, beta0=10)