BCPA/Model Specification
From QERM Wiki
(Difference between revisions)
m (moved BCPA/Model Selection to BCPA/Model Specification) |
|||
Line 206: | Line 206: | ||
} | } | ||
</pre> | </pre> | ||
+ | [[Category: BCPA]] |
Latest revision as of 02:03, 23 August 2011
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)) }