Nothing
require(glasso)
der.scad <- function(A,a,lambda) {
d = dim(A)[1]
dpen = matrix(1,d,d)
wp = which(abs(A)>lambda)
lwp = length(wp)
dpen[wp] = apply(matrix(c(a*lambda-abs(A[wp]),rep(0,lwp)),lwp,2),1,max)/(a-1)/lambda
return(dpen)
}
ll_reg <- function(A,S,Cxx,Cxy,Cyy){ # -ll
d = dim(S)[1]
A = matrix(A,d,d)
Si = solve(S)
return(sum(diag((Cyy-Cxy%*%t(A)-A%*%t(Cxy)+A%*%Cxx%*%t(A))%*%Si)))
}
grad.ll_reg <- function(A,S,Cxx,Cxy,Cyy){ # -grad.ll
d = dim(S)[1]
A = matrix(A,d,d)
Si = solve(S)
pdt1 = 2*Si%*%t(Cxy)
pdt2 = t(2*Cxx%*%t(A)%*%Si)
return(pdt2-pdt1)
}
hess.ll_reg = function(A,S,Cxx,Cxy,Cyy){ # -hess.ll
d = dim(S)[1]
Si = solve(S)
hess = matrix(0,d*d,d*d)
cnt=0
for (ii in 1:d){
for (jj in 1:d) {
cnt = cnt+1
hess[cnt,] = c(t(Cxx[,ii]%*%t(Si[jj,])))*2 # by symetry (jj-1)*d+ii
}
}
return(hess)
}
ll.pen2 <- function(A,S,Cxx,Cxy,Cyy,omega){ # -ll
d = dim(S)[1]
A = matrix(A,d,d)
Si = solve(S)
return(sum(diag((Cyy-Cxy%*%t(A)-A%*%t(Cxy)+A%*%Cxx%*%t(A))%*%Si)) + sum(diag(omega%*%A%*%t(A))))
}
grad.ll.pen2 <- function(A,S,Cxx,Cxy,Cyy,omega){ #-grad.ll
d = dim(S)[1]
A = matrix(A,d,d)
Si = solve(S)
pdt1 = 2*Si%*%Cxy
pdt2 = t(2*Cxx%*%t(A)%*%Si)
return(pdt2-pdt1+(omega+t(omega))%*%A)
}
# the unified algortihm Fan 2001
Mstep.hh.SCAD.MSAR <-
function(data,theta,FB,lambda1=.1,lambda2=.1,penalty="SCAD",par=NULL) {
Miter = 5
T=dim(data)[1]
N.samples = dim(as.array(data))[2]
d = dim(as.array(data))[3]
if(is.null(d)|is.na(d)) {d = 1}
M <- attributes(theta)$NbRegimes
p <- attributes(theta)$order
order <- max(p,1)
if (length(lambda1)==1) {lambda1 = matrix(lambda1,1,M)*matrix(theta$prior,1,M)}
if (length(lambda2)==1) {lambda2 = matrix(lambda2,1,M)*matrix(theta$prior,1,M)}
data = array(data,c(T,N.samples,d))
data2 = array(0,c(order*d,T-order+1,N.samples))
cpt=1
for (o in order:1) {
for (kd in 1:d) {
data2[cpt,,] = data[o:(T-order+o),,kd]
cpt =cpt+1
}
}
#browser()
T = length(o:(dim(data)[1]-order+o))
N = T*N.samples
exp_num_trans = 0
exp_num_visit = 0
exp_num_visits1 = 0
postmix = 0
m = matrix(0,d*order,M) ; m_1 = m ; c = matrix(0,M,1) ; s=c ;
op = array(0,c(d*order,d*order,M) ); op_1 = op ; op_2 = op_1 ;
Cxx = array(0,c(d*order,d*order,M) )
Cxy = array(0,c(d*order,d*order,M) )
Cyy = array(0,c(d*order,d*order,M) )
for (ex in 1:N.samples) {
obs = t(array(data2[,,ex],c(d*order,T)))
xit = array(FB$probSS[,,,ex],c(M,M,T-2))
gamma = matrix(FB$probS[ex,1:(T-1),],T-1,M)
exp_num_trans = exp_num_trans+apply(xit,c(1,2),sum)
exp_num_visits1 = exp_num_visits1+gamma[1,]
postmix = postmix+apply(gamma,2,sum)
for (j in 1:M) { # attention, ici on n'a qu'un ex
w = matrix(gamma[,j], 1,T-1)
wobs = obs[1:(T-1),] * repmat(t(w),1,d*order)
wobs_1= obs[2:T,] * repmat(t(w),1,d*order)
if (is.na(sum(obs))) {
wobs[is.na(wobs)] = 0
wobs_1[is.na(wobs_1)] = 0
obs[is.na(obs)] = 0
}
m[,j] = m[,j] + apply(wobs,2,sum)
m_1[,j] = m_1[,j] + apply(wobs_1, 2,sum)
op[,,j] = op[,,j] + t(wobs) %*% obs[1:(T-1),]
op_1[,,j] = op_1[,,j] + t(wobs) %*% obs[2:T,]
op_2[,,j] = op_2[,,j] + t(wobs_1) %*% obs[2:T,]
}
}
# Cxx = Cxx/N.samples
# Cxy = Cxy/N.samples
# Cyy = Cyy/N.samples
# Markov chain
prior = normalise(exp_num_visits1)
prior=matrix(prior,M,1)
transmat = mk_stochastic(exp_num_trans)
if (min(postmix)<1e-6) {stop("error : smoothing probabilities are to small, in one regime at least. You should revise initialisation.")}
moy <- array(0,c(M,d))
Areg=list()
Sreg = list()
sigma.inv=list()
a = 3.7
for (j in 1:M) { # attention, ici on n'a qu'un ex
Cxx[,,j] = postmix[j]*op[,,j] - m[,j]%*%t(m[,j])
Cxy[,,j] = postmix[j]*op_1[,,j] - m[,j]%*%t(m_1[,j])
Cyy[,,j] = postmix[j]*op_2[,,j] - m_1[,j]%*%t(m_1[,j])
if (det(Cxx[,,j])<1e-20){Cxx[,,j] = Cxx[,,j]+diag(1e-8,dim(Cxx[,,j])[1])}
#browser()
A2tmp = t(Cxy[,,j])%*%solve(Cxx[,,j])
A = theta$A[[j]][[1]]
S = theta$sigma[[j]]
S.emp = (Cyy[,,j]-Cxy[,,j]%*%t(A)-A%*%t(Cxy[,,j])+A%*%Cxx[,,j]%*%t(A))/postmix[j]^2
res.gl = list()
#S = S.emp
#res.gl$wi = solve(S)
if (!is.null(par$sigma.inv)) {res.gl$wi = par$sigma.inv[[j]]
} else {res.gl$wi = solve(S)}
w = matrix(0,d,d)
#browser()
if (lambda1[j]>0 & penalty=="SCAD"){
w.tmp = apply(matrix(c(a*lambda1[j]-abs(res.gl$wi),rep(0,d*d)),d*d,2),1,max)
w = (lambda1[j]*(as.numeric(abs(res.gl$wi)<=lambda1[j])+w.tmp/((a-1)*lambda1[j])*as.numeric(abs(res.gl$wi)>lambda1[j])))
} else {w = matrix(lambda1[j],d,d)}
if (lambda2[j]>0 & penalty=="SCAD") {
nzA = which(abs(A)>lambda2[j]) # a revoir?
#print(paste("Nombre de 0 dans A",sum(A==0)," et de petits coefficients",d*d-length(nzA),sep=" "))
iter = 0 ; err=1
while (iter<Miter & err>1e-3) {
A.old=A
pen.term1 = der.scad(A,a,lambda2[j])
pen.term2 = diag(pen.term1[nzA]/abs(A)[nzA])
H = hess.ll_reg(A,S,Cxx[,,j],Cxy[,,j],Cyy[,,j])[nzA,nzA]+N*pen.term2
if (length(H)==0) {stop(paste("In regime",j,", lambda2 = ",lambda2[j],"and is to strong."))}
#A[nzA] = c(A)[nzA]-solve(H) %*% matrix(c(grad.ll_reg(A,S,Cxx[,,j],Cxy[,,j],Cyy[,,j]))[nzA]+N*diag(pen.term2)*A[nzA],length(nzA),1)
A[nzA] = c(A)[nzA]-solve(H,matrix(c(grad.ll_reg(A,S,Cxx[,,j],Cxy[,,j],Cyy[,,j]))[nzA]+N*diag(pen.term2)*A[nzA],length(nzA),1))
err = norm(A-A.old)
iter = iter+1
}
A[-nzA] = 0
} else if (lambda2[j]>0 & penalty=="ridge") {
omega=matrix(lambda2[j],d,d)
omega = omega-diag(diag(omega))
res=ucminf(c(A), fn=ll.pen2,gr=grad.ll.pen2,S=S,Cxx=Cxx[,,j],Cxy=Cxy[,,j],Cyy=Cyy[,,j],omega=omega)
A = res$par
} else {A = A2tmp}
Areg[[j]] = list()
Areg[[j]][[1]] = A
if (lambda1[j]>0) {
w = N*w
res.gl = glasso(S, rho=matrix(w,d,d),penalize.diagonal=FALSE,maxit=10)
sigma.inv[[j]] = res.gl$wi
Sreg[[j]] = res.gl$w}
else {
tmp = (m_1[,j]-(A)%*%m[,j])/postmix[j]
op_1[,,j] = t(as.matrix(op_1[,,j]))
Sreg[[j]] = (op_2[, , j] + (A) %*% op[, , j] %*% t(A) -
((A) %*% t(op_1[, , j]) + t((A) %*% t(op_1[, , j]))))/postmix[j] - tmp %*% t(tmp)
sigma.inv[[j]] = solve(Sreg[[j]])
}
moy[j,] = (m_1[,j]-A%*%m[,j])/postmix[j]
# w = N*w
# res.gl = glasso(S, rho=matrix(w,d,d),penalize.diagonal=FALSE,maxit=10)
# sigma.inv[[j]] = res.gl$wi
}
if (p>0) {
list(A=Areg,A0=moy,sigma=Sreg,prior=prior,transmat=transmat,sigma.inv=sigma.inv)
} else {
list(sigma=Sreg,A0=moy,prior=prior,transmat=transmat,sigma.inv=sigma.inv)
}
}
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.