Nothing
init.theta.miss.MSAR <-
function(data,...,M,order,regime_names=NULL,nh.emissions=NULL,nh.transitions=NULL,label=NULL,ncov.emis = 0,ncov.trans=0,cl.init="mean"
) {
if (missing(M) || is.null(M) || M==0) {print("Need at least one regime : M=1"); M <- 1 }
if (missing(order)) { order <- 0 } # AR order
if (missing(label)) {label = 'HH'}
T = dim(data)[1]
N.samples = dim(data)[2]
d = dim(data)[3]
if (is.null(d) ) {d=1}
else if ( is.na(d)) {d=1}
#if (d==1) {data <- as.matrix(data)}
if (length(dim(data))<3) {d = 1}
else {d = dim(data)[3]}
if (cl.init == "mean") {
d.data = matrix(0,(T-1)*N.samples,d)
data.mat = matrix(0,T*N.samples,d)
for (id in 1:d) {
d.data[,id] = c(data[2:T,,id]-data[1:(T-1),,id])
data.mat[,id] = c(data[,,id])
}
if (is.na(sum(data))) {
w = which(is.na(apply(d.data,1,sum)))
d.data = matrix(d.data[-w,],(T-1)*N.samples-length(w),d)
w = which(is.na(apply(data.mat,1,sum)))
data.mat = matrix(data.mat[-w,],(T)*N.samples-length(w),d)
}
if (M>1) {
centers=d.data[sample(1:dim(d.data)[1],M),]
if (any(duplicated(centers))) {
cnt=0
while (any(duplicated(centers)) & cnt<5){
cnt = cnt+1
centers = d.data[sample(1:dim(d.data)[1],M),]
}
if (any(duplicated(centers))) {print("Is there a lot of equal data?")}
}
mc = kmeans(d.data,centers=centers)
class = mc$cluster
}
else {class = 1}
}
else {
d.var = array(0,c((T-1),N.samples,d))
ht = 3
for (id in 1:d) {
for (t in 1:(T-1)) {
ii = max(t-ht,1):min(T-1,t+ht)
d.var[t,,id] = apply(data[ii+1,,id]-data[ii,,id],2,var)
}
}
qv = quantile(d.var,probs=(1:M)/M)
d.var = matrix(d.var,(T-1)*N.samples,d)
if (M>1) {
mc = kmeans(d.var,centers=d.var[sample(1:dim(d.var)[1],M),])
class = mc$cluster
}
}
if (order>0) {
#print(paste("order =",order))
res = Mstep.classif(data,array(class,c(T,N.samples,1)),order=order)
A = res$A
A0 = res$A0
sigma = res$sigma
} else {
A = sigma = list()
A0 = matrix(NA,M,d)
for (m in 1:M){
w = which(class==m)
A0[m,] = apply(matrix(data.mat[w,],length(w),d),2,mean,na.rm =TRUE)
sigma[[m]] = cov(matrix(data.mat[w,],length(w),d),use="complete.obs")
}
}
#for (){if (det(sigma)<0) {sigma = sigma+1e-4*diag(1,d)}}
#TRANSITION PROBABILITY MATRIX / INITIAL DISTRIBUTION
if (M>1) {
prior <- normalise(runif(M))
transmat <- mk_stochastic(diag(1,M)+matrix(runif(M*M),M)) # homogeneous Markov chain
} else {
prior = 1
transmat = 1
}
prior=matrix(prior,M,1)
transmat=matrix(transmat,M,M)
n_par= M*(M-1)+(order*d^2+d+d*(d+1)/2)*M
#EMISSION PARAMETERS
emis.linear = FALSE
if (substr(label,2,2)=="N") {
par.emis <- list()
if (missing(nh.emissions)) {nh.emissions = 'linear'}
if (!is.function(nh.emissions)){
if (nh.emissions == 'linear') {
emis.linear = TRUE
if (ncov.emis<1) {ncov.emis=1}
nh.emissions <- function(covar,par.emis){
d <- dim(par.emis)[1]
if(is.null(d) || is.na(d)){d <- 1}
f <- matrix(0,d,dim(covar)[1]) # if covar is a vector?
for(i in 1:d){
f[i,] <- par.emis[i,1:dim(par.emis)[2]]%*%t(covar)
}
return(f)
}
}
}
for(i in 1:M){
par.emis[[i]] <- matrix(0,d,ncov.emis)
}
n_par <- n_par+ncov.emis
}
#TRANSITION PARAMETERS
if (substr(label,1,1)=="N") {
par.trans <- array(1,c(M,max(2,ncov.trans+1)))
if (missing(nh.transitions)) { nh.transition = 'gauss'}
if (!is.function(nh.transitions)){if (nh.transitions=='gauss') { nh.transitions <- function(covar,par.trans,transmat){
T <- dim(covar)[1]
N.samples = dim(covar)[2]
ncov = dim(covar)[3]
M = dim(transmat)[1]
f <- array(0,c(M,M,T))
p.t = array(0,c(T,N.samples,ncov))
for (j in 1:M) {
for (nc in 1:ncov) {p.t[,,nc] = par.trans[j,nc]}
xx = (covar-p.t)^2
sxx = apply(xx,1,sum)
temp=exp( -.5*sxx/par.trans[j,ncov+1]^2)/par.trans[j,ncov+1]
for (i in 1:M) {f[i,j,] = transmat[i,j]*temp }
}
f.sum = apply(f , c(1,3), sum)
for (i in 1:M) {f[i,,] = f[i,,]/t(matrix(f.sum[i,],T,M))}
return(f)
}
}
else if (nh.transitions=='logistic') {
nh.transitions <- function(covar,par.trans,transmat){
eps = 1e-10
T <- dim(covar)[1]
nc <- dim(covar)[3]
covar = matrix(covar,T,nc)
M = dim(transmat)[1]
f <- array(0,c(M,M,T))
for (m in 1:M) {
f[m,m,] = eps+(1-2*eps)/(1 + exp(par.trans[m,1] + par.trans[m,2:(nc+1)] %*% t(covar)))
}
f[1,2,] = 1-f[1,1,]
f[2,1,] = 1-f[2,2,]
return(f)
}
}
else if (nh.transitions=='VM') {
nh.transitions <- function(covar,par.trans,transmat){
T <- dim(covar)[1]
N.samples = dim(covar)[2]
ncov = dim(covar)[3]
M = dim(transmat)[1]
par.trans = repmat(t(par.trans[,2])*exp(1i*par.trans[,1]),M,1)%*%(matrix(1,M,M)-diag(1,nrow=M))
f <- array(0,c(M,M,T))
for (j in 1:M) {
for (i in 1:M) {
f[i,j,] = transmat[i,j]*abs(exp(par.trans[i,j]*exp(-1i*covar))) ;
}
}
f.sum = apply(f , c(1,3), sum)
for (i in 1:M) {f[i,,] = f[i,,]/t(matrix(f.sum[i,],T,M))}
return(f)
}
} else if (nh.transitions=='probit') {
nh.transitions <- function(covar,par.trans,transmat){
T <- dim(covar)[1]
N.samples = dim(covar)[2]
ncov = dim(covar)[3]
M = dim(transmat)[1]
f <- array(0,c(M,M,T))
for (j in 1:M) {
temp=pnorm(par.trans[j,1]+par.trans[j,2:(ncov+1)]%*%t(as.matrix(covar)))
for (i in 1:M) {
f[i,j,] = transmat[i,j]*temp ; # On perd la monotonie...
#f[i,j,] = temp ;
}
}
f.sum = apply(f , c(1,3), sum)
for (i in 1:M) {f[i,,] = f[i,,]/t(matrix(f.sum[i,],T,M))}
return(f)
}
}
}
n_par <- n_par+max(2,ncov.trans+1)
}
if (order==0 && label=='HH'){theta=list(A0,sigma,prior,transmat)}
if (order==0 && label=='HN'){theta=list(A0,sigma,prior,transmat,par.emis)}
if (order==0 && label=='NH'){theta=list(A0,sigma,prior,transmat,par.trans)}
if (order==0 && label=='NN'){theta=list(A0,sigma,prior,transmat,par.trans,par.emis)}
else if (order>0 && label=='HH'){theta=list(A,A0,sigma,prior,transmat)}
else if (order>0 && label=='HN'){theta=list(A,A0,sigma,prior,transmat,par.emis)}
else if (order>0 && label=='NH'){theta=list(A,A0,sigma,prior,transmat,par.trans)}
# else if (order>0 && label=='NN'){theta=list(A,A0,sigma,prior,transmat,par.emis,par.trans)}
else if (order>0 && label=='NN'){theta=list(A,A0,sigma,prior,transmat,par.trans,par.emis)}
attr(theta,'NbComp') <- d
attr(theta,'NbRegimes') <- M
attr(theta,'order') <- order
attr(theta,'label') <- label
attr(theta,'nh.emissions') <- nh.emissions
attr(theta,'nh.transitions') <- nh.transitions
attr(theta,'n_par') <- n_par
attr(theta,'emis.linear') <- emis.linear
theta=as.thetaMSAR(theta,label=label,ncov.emis=ncov.emis,ncov.trans=ncov.trans)
#class(theta) <- "MSAR"
#theta$call <- match.call()
return(theta)
}
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.