picard <- function(x, ...)
{
UseMethod("picard")
}
picard.default <- function(x,P,Kmax, pos=NULL, ...)
{
if (ncol(x) !=2)
stop("not yet implemented for more or less than two variables")
if (!is.null(pos)) {
if (length(pos)!=nrow(x))
stop("pos should have the same length as x")
}
x <- as.matrix(x)
xor <- x
cons <- apply(x, 1, function(y) all(!is.na(y)))
nna.places <- c(1:nrow(x))[cons]
x <- x[nna.places,]
Linc <- matrix(-Inf,nrow=Kmax,ncol=1)
n <- dim(x)[2]
param <- list()
Kmin <- P
if (P==1){
rupt <- c(1, n)
phi <- list(mu=matrix(rowMeans(x),ncol=1),
sigma <- matrix(apply(x,1,sd),ncol=1),
prop=1)
Linc[Kmin:Kmax] <- logdens_simultanee(x,phi)
param <- list(phi=phi,rupt=rupt, tau=NA)
resu <- list(Linc=Linc, param=param)
} else {
for (i in 1:ncol(x))
x[,i] <- as.double(x[,i])
resu <- .Call("hybridSimultanee", x, as.integer(P), as.integer(Kmax), PACKAGE="segTraj")
names(resu) <- c("Linc","param")
nonnul <- c(1:length(resu$param))[!sapply(resu$param, is.null)]
for (i in nonnul) {
names(resu$param[[i]]) <- c("phi","rupt","tau")
resu$param[[i]]$rupt <- resu$param[[i]]$rupt+1L
names(resu$param[[i]]$phi) <- c("mu","sigma","prop")
}
}
attr(resu, "P") <- P
attr(resu, "Kmax") <- Kmax
attr(resu, "seq.ori") <- xor
attr(resu, "pos") <- pos
attr(resu, "nna.places") <- nna.places
class(resu) <- "picard"
attr(resu, "typeseg") <- "default"
return(resu)
}
print.picard <- function(x, ...)
{
cat("=========\n\nSegmentation of a series using the method of Picard\n",
"based on P =", attr(x, "P"),"different models",
"with Kmax = ", attr(x, "Kmax"), ".\n\n=========\n",
"This object is a list containing\n",
"$Linc: the likelihood for the best segmentation with K-segment\n",
"$param: a list containing the parameters of the segmentation\n",
" with one element per K; each element is itself a list containing:\n",
"..$phi: the parameters of the models. A list containing:\n",
"......$mu: the mean of each variable (column) for each model (row)\n",
"......$sigma: the standard deviation of each variable (column) for each model (row)\n",
"......$prop: the proportion of observations generated by each model\n",
"..$rupt: the indices of the limits (start, end) of the K segments (one per row)\n",
"..$tau: the likelihood of each model (row) for each segment (column)\n\n")
}
picard.ltraj <- function(x, P, Kmax, which=c("x","y"), ...)
{
if (!inherits(x, "ltraj"))
stop("ltraj should be of class \"ltraj\"")
if (length(x)>1)
stop("only implemented for one trajectory")
if (length(which)!=2)
stop("for the moment, only implemented for two variables")
y <- x[[1]]
if (!is.null(attr(y, "infolocs")))
y <- cbind(y, attr(y, "infolocs"))
if (!all(which%in%names(y)))
stop("The variables in which were not found")
coin <- y[,which]
cons <- apply(coin, 1, function(x) all(!is.na(x)))
nna.places <- c(1:nrow(coin))[cons]
coin <- coin[nna.places,]
pic <- picard(x=coin, P=P, Kmax=Kmax)
attr(pic, "nna.places") <- nna.places
attr(pic, "seq.ori") <- x
attr(pic, "which") <- which
attr(pic, "typeseg") <- "ltraj"
return(pic)
}
plot.picard <- function(x, axes = FALSE, number=NULL,
addparam=FALSE, bg=TRUE,...)
{
if (!inherits(x, "picard"))
stop("x should be of class picard")
typex <- attr(x, "typeseg")
nna <- attr(x, "nna.places")
## prepare the two series and the x-coordinate
if (typex == "ltraj") {
xk <- attr(x, "seq.ori")
y <- xk[[1]]
pos <- y$date
if (!is.null(attr(y, "infolocs")))
y <- cbind(y, attr(y, "infolocs"))
xk <- y[nna,attr(x,"which")]
serie <- pos[nna]
} else {
xk <- attr(x, "seq.ori")
pos <- attr(x, "pos")
if (!is.null(pos)) {
serie <- pos[nna]
} else {
serie <- c(1:nrow(xk))[nna]
}
xk <- xk[nna,]
}
## ylim
moxk <- apply(xk, 2, mean)
sdxk <- apply(xk, 2, sd)
xk <- scale(xk)
ra <- apply(xk,2,range)
ra <- c(min(ra[1,]),max(ra[2,]))
## Number of plots
nonnul <- c(1:length(x$param))[!sapply(x$param, is.null)]
if (axes) {
opar <- par(mfrow = c(length(nonnul),1), mar=c(3,3,2,1))
} else {
opar <- par(mfrow = c(length(nonnul),1), mar=c(0,0,2,0))
}
on.exit(par(opar))
## Number of models
P <- attr(x, "P")
colo <- grey((1:P)/(P+1))
for (i in nonnul) {
plot(serie, xk[,1], type="n", axes=FALSE, main=paste("K =", i), xlab="", ylab="",
ylim = ra)
rupt <- x$param[[i]]$rupt
## the most likely model
modt <- apply(x$param[[i]]$tau,1,which.max)
if (bg) {
tmp <- sapply(1:nrow(rupt), function(j) {
## the polygon
pol <- cbind(c(serie[rupt[j,1]],serie[rupt[j,2]],
serie[rupt[j,2]],serie[rupt[j,1]]),
c(ra[1], ra[1], ra[2], ra[2]))
m <- modt[j]
polygon(pol[,1],pol[,2], col=colo[m])
})
}
## the two series
if (!is.null(number)){
if (addparam) {
m <- sapply(1:nrow(rupt), function(j) {
mo <- (x$param[[i]]$phi$mu[modt[j],number]-moxk[number])/sdxk[number]
si <- x$param[[i]]$phi$sigma[modt[j],number]/sdxk[number]
segments(serie[rupt[j,1]], mo,
serie[rupt[j,2]], mo, col="black", lwd=2)
pol <- cbind(c(serie[rupt[j,1]],serie[rupt[j,2]],
serie[rupt[j,2]],serie[rupt[j,1]]),
c(mo-si, mo-si, mo+si,mo+si))
polygon(pol[,1],pol[,2], density=10)
})
}
lines(serie, xk[,number], lwd=3, col="black")
lines(serie, xk[,number], col="white")
} else {
if (addparam) {
for (numberb in 1:2) {
m <- sapply(1:nrow(rupt), function(j) {
mo <- (x$param[[i]]$phi$mu[modt[j],numberb]-moxk[numberb])/sdxk[numberb]
si <- x$param[[i]]$phi$sigma[modt[j],numberb]/sdxk[numberb]
segments(serie[rupt[j,1]], mo,
serie[rupt[j,2]], mo, col="black", lwd=2)
pol <- cbind(c(serie[rupt[j,1]],serie[rupt[j,2]],
serie[rupt[j,2]],serie[rupt[j,1]]),
c(mo-si, mo-si, mo+si,mo+si))
polygon(pol[,1],pol[,2], density=c(20,10)[numberb],
col=c("black","red")[numberb])
})
}
}
lines(serie, xk[,1], lwd=3, col="black")
lines(serie, xk[,1], col="white")
lines(serie, xk[,2], lwd=3, col="red")
lines(serie, xk[,2], col="white")
}
## date:
if (axes) {
axis(1)
axis(2)
}
box()
}
return(invisible(NULL))
}
testConvexHull <- function(y)
{
n <- length(y)
print(y)
print(n)
resu <- .Call("testConvexHull", as.numeric(y), as.integer(n), PACKAGE="segTraj")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.