BCPA/Model Specification

From QERM Wiki
Jump to: navigation, search
	M0 <- function(x,t,tbreak,K=2)
	# null model: all mus, s's, rhos the same
	{
		rhoLL <- GetRho(x,t)
		LL <- rhoLL[2]
		
		bic <- K*LL + 3*log(length(x))
		
		rho1<-rhoLL[1]
		rho2<-rho1
		
		mu1<-mean(x)
		mu2<-mu1
		
		s1<-sd(x)
		s2<-s1
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}
	
	M1 <- function(x,t,tbreak,K=2)
	# mus different, all else the same
	{
		
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1<-mean(x1)
		mu2<-mean(x2)
		
		xprime <- c(x1-mu1,x2-mu2)
		s1<-sd(xprime)
		s2<-s1
		
		rho1<-as.numeric(GetRho(xprime,t)[1])
		rho2<-rho1
		
		LL1<-GetLL(x1,t1,mu1,s1,rho1)
		LL2<-GetLL(x2,t2,mu2,s2,rho2)
		LL<-LL1+LL2
		bic <- K*LL + 5*log(length(x))
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}
	
	M2 <- function(x,t,tbreak,K=2)
	# sds different, all else same
	{
		
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1<-mean(x)
		mu2<-mu1
		
		s1<-sd(x1)
		s2<-sd(x2)
		
		xprime <- c( (x1-mu1)/s1 , (x2-mu2)/s2 )
		rho1 <- as.numeric(GetRho(xprime,t)[1])
		rho2 <- rho1
		
		LL1<-GetLL(x1,t1,mu1,s1,rho1)
		LL2<-GetLL(x2,t2,mu2,s2,rho2)
		LL<-LL1+LL2
		bic <- K*LL + 5*log(length(x))
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}	
		
	M3 <- function(x,t,tbreak,K=2)
	# rhos different, all else same
	{
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1 <- mean(x)
		mu2 <- mu1
		
		s1 <- sd(x)
		s2 <- s1
		
		rho1<-as.numeric(GetRho(x1,t1)[1])
		rho2<-as.numeric(GetRho(x2,t2)[1])
				
		LL1<-GetLL(x1,t1,mu1,s1,rho1)
		LL2<-GetLL(x2,t2,mu2,s2,rho2)
		LL<-LL1+LL2
		bic <- K*LL + 5*log(length(x))
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}		
	
	M4 <- function(x,t,tbreak,K=2)
	# mu and sigma different, rho same
	{
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1<-mean(x1)
		mu2<-mean(x2)
		s1<-sd(x1)
		s2<-sd(x2)
		
		xprime <- c( (x1-mu1)/s1 , (x2-mu2)/s2 )
		rho1 <- as.numeric(GetRho(xprime,t)[1])
		rho2 <- rho1
			
		LL1<-GetLL(x1,t1,mu1,s1,rho1)
		LL2<-GetLL(x2,t2,mu2,s2,rho2)
		LL<-LL1+LL2
		bic <- K*LL + 6*log(length(x))
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}		

	M5 <- function(x,t,tbreak,K=2)
	# mu and rho different, sigma same
	{
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1<-mean(x1)
		mu2<-mean(x2)
		
		xprime <- c(x1-mu1, x2-mu2)
		s1<-sd(xprime)
		s2<-s1
				
		rho1<-as.numeric(GetRho(x1,t1)[1])
		rho2<-as.numeric(GetRho(x2,t2)[1])
		
		LL1<-GetLL(x1,t1,mu1,s1,rho1)
		LL2<-GetLL(x2,t2,mu2,s2,rho2)
		LL<-LL1+LL2
		bic <- K*LL+ 6*log(length(x))
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}		

	M6 <- function(x,t,tbreak,K=2)
	# sigma and rho different, mu same
	{
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1<-mean(x)
		mu2<-mean(x)
		s1<-sd(x1)
		s2<-sd(x2)
		
		x1prime <- (x1-mu1)/s1 
		x2prime <- (x2-mu2)/s2
		
		rho1 <- as.numeric(GetRho(x1prime,t1)[1])
		rho2 <- as.numeric(GetRho(x2prime,t2)[1])
			
		LL1<-GetLL(x1,t1,mu1,s1,rho1)
		LL2<-GetLL(x2,t2,mu2,s2,rho2)
		
		LL<-LL1+LL2
		
		bic <- K*LL+ 6*log(length(x))
		
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}	

	M7 <- function(x,t,tbreak,K=2)
	# most "alternative" model: all mus, s's, rhos different
	{
		rhoLL1 <- GetRho(x[1:tbreak],t[1:tbreak])
		rhoLL2 <- GetRho(x[(tbreak+1):length(x)],t[(tbreak+1):length(x)])
		
		LL1 <- rhoLL1[2]
		LL2 <- rhoLL2[2]
		
		x1<-x[1:tbreak]
		x2<-x[(tbreak+1):length(x)]
		t1<-t[1:tbreak]
		t2<-t[(tbreak+1):length(x)]
		
		mu1<-mean(x1)
		mu2<-mean(x2)
		s1<-sd(x1)
		s2<-sd(x2)
		rho1 <- rhoLL1[1]
		rho2 <- rhoLL2[1]
		
		LL <- LL1+LL2
		bic <- K*LL + 7*log(length(x))
		return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
	}
Personal tools