Nothing
invgaussMLE <- function (yi, ni = numeric(length(yi)) + 1,
si = numeric(length(yi)) + 1,
parameterization = "sigma2"
) {
## check if yi is a spikeTrain object if yes take the "diff"
if (inherits(yi,"spikeTrain")) yi <- diff(yi)
if (any(yi < 0))
stop("yi elements must be non-negative")
yi <- as.numeric(yi)
if (!identical(length(yi), length(ni)))
stop("yi and ni should have the same length")
if (any(ni < 0))
stop("ni elements must be non-negative")
if (!identical(class(ni), "integer") && !identical(ni, round(ni)))
stop("ni should be a vector of positive integers")
if (!identical(length(yi), length(si)))
stop("yi and si should have the same length")
if (any(si < 0))
stop("si elements must be non-negative")
if (!identical(si, round(si)))
stop("si should be a vector of positive integers")
if (any(si > ni))
stop("si elements should not be greater than ni elements")
minusLogLik <- function(p) {
if (missing(p)) {
txt <- paste("This function argument should be a 2 component vector:\n",
" component 1 is the log of the mu parameter,\n",
" component 2 is the log of the sigma2 parameter,\n",
"using the 'sigma2' parameterization of the inverse Gaussian.\n")
cat(txt)
} else {
mu <- exp(p[1])
sigma2 <- exp(p[2])
-(ifelse(s.dot > 0,
sum(dinvgauss(yi[si > 0], mu, sigma2, log = TRUE) * si[si > 0], 0)
) +
ifelse(c.dot > 0,
sum(pinvgauss(yi[ci > 0], mu, sigma2, lower.tail = FALSE,log.p = TRUE) * ci[ci > 0]), 0))
}
}
s.dot <- sum(si)
n.dot <- sum(ni)
ci <- ni - si
c.dot <- sum(ci)
if (s.dot == n.dot) {
mu.hat <- weighted.mean(yi, ni)
inv.mean <- weighted.mean(1/yi, ni)
sigma2.hat <- inv.mean - 1/mu.hat
estimate <- c(mu.hat, sigma2.hat)
observedI <- matrix(c(n.dot/(sigma2.hat * mu.hat^3),
0, 0, n.dot/sigma2.hat^2/2), nrow = 2, byrow = TRUE)
se <- sqrt(diag(solve(observedI)))
l <- -minusLogLik(log(estimate))
} else {
if (s.dot >= 10) {
mu.hat <- weighted.mean(yi, si)
inv.mean <- weighted.mean(1/yi, si)
sigma2.hat <- inv.mean - 1/mu.hat
} else {
mu.hat <- weighted.mean(yi, ni)
inv.mean <- weighted.mean(1/yi, ni)
sigma2.hat <- inv.mean - 1/mu.hat
}
mleFit <- optim(par=log(c(mu.hat, sigma2.hat)),
fn=minusLogLik,
method="BFGS",
hessian = TRUE)
estimate <- exp(mleFit$par)
newVar <- (1/estimate) %o% (1/estimate)
observedI <- mleFit$hessian * newVar
se <- sqrt(diag(solve(observedI)))
l <- -mleFit$value
}
if (parameterization == "sigma2") {
names(estimate) <- c("mu", "sigma2")
names(se) <- c("mu", "sigma2")
rFct <- function(mu, sigma2) -minusLogLik(log(c(mu, sigma2))) - l
} else {
boundary.hat <- (sigma2.hat)^(-0.5)
mu.hat <- mu.hat/boundary.hat
if (s.dot == n.dot) {
estimate <- c(mu.hat, boundary.hat)
observedI <- observedI * matrix(c(boundary.hat^2,
-2/boundary.hat^2, -2/boundary.hat^2, 4/boundary.hat^6),
nrow = 2, byrow = TRUE)
se <- se * c(1/boundary.hat, boundary.hat^3/2)
} else {
estimate <- c(estimate[1] * sqrt(estimate[2]), 1/sqrt(estimate[2]))
newVar <- newVar * (c(estimate[2], -2/estimate[2]^3) %o%
c(estimate[2], -2/estimate[2]^3))
observedI <- mleFit$hessian * newVar
se <- sqrt(diag(solve(observedI)))
}
names(estimate) <- c("mu", "boundary")
names(se) <- c("mu", "boundary")
rFct <- function(mu, boundary) -minusLogLik(log(c(mu * boundary, 1/boundary^2))) - l
}
result <- list(estimate = estimate,
se = se,
logLik = l,
r = rFct,
mll = minusLogLik,
call = match.call()
)
class(result) <- "durationFit"
return(result)
}
coef.durationFit <- function(object,...) object$estimate
logLik.durationFit <- function(object,...) object$logLik
gammaMLE <- function(yi,
ni = numeric(length(yi))+1,
si = numeric(length(yi))+1,
scale = TRUE
) {
## check if yi is a spikeTrain object if yes take the "diff"
if (inherits(yi,"spikeTrain")) yi <- diff(yi)
## check that yi elements are positive
if (any(yi < 0)) stop("yi elements must be non-negative")
## coerce yi to vector
yi <- as.numeric(yi)
## check that ni has the same length as yi
if (!identical(length(yi),length(ni))) stop("yi and ni should have the same length")
## check that the elements of ni are non negative and integer
if (any(ni < 0)) stop("ni elements must be non-negative")
if (!identical(ni, round(ni))) stop("ni should be a vector of positive integers")
if (!identical(length(yi),length(si))) stop("yi and si should have the same length")
## check that the elements of ni are non negative and integer
if (any(si < 0)) stop("si elements must be non-negative")
if (!identical(si, round(si))) stop("si should be a vector of positive integers")
if (any(si > ni)) stop("si elements should not be greater than ni elements")
## Create function returning the opposite of the log likelihood
minusLogLik <- function(p) {
if (missing(p)) {
txt <- paste("This function argument should be a 2 component vector:\n",
" component 1 is the log of the shape parameter,\n",
" component 2 is the log of the scale parameter.\n"
)
cat(txt)
} else {
shape <- exp(p[1])
scale <- exp(p[2])
-(ifelse(s.dot>0,sum(dgamma(yi[si>0],shape=shape,scale=scale,log=TRUE)*si[si>0],0))+
ifelse(c.dot>0,sum(pgamma(yi[ci>0],shape=shape,
scale=scale,lower.tail=FALSE,log.p=TRUE)*ci[ci>0]),0)
)
}
}
## Get the number of uncensored events
s.dot <- sum(si)
## Get the total number of events
n.dot <- sum(ni)
ci <- ni-si
c.dot <- sum(ci)
## Define function logProfileLik returning the log of the
## profiled likelihood for the shape parameter
logProfileLik <- function(shape) {
if (shape <= 0) return(-Inf)
term1 <- s.dot*shape*(log(shape) - log(scale.hat) -1)
term2 <- (shape - 1) * sum(ni * log(yi))
term3 <- - s.dot * lgamma(shape)
return(term1 + term2 + term3)
}
if (s.dot == n.dot) {
## no censored event
## Get the empirical mean which is the MLE of parameter scale
scale.hat <- weighted.mean(yi, ni)
s2 <- weighted.mean(yi^2, ni) - scale.hat^2
shape.hat <- scale.hat^2 / s2
shape.lower <- floor(shape.hat)
shape.upper <- ceiling(shape.hat)
shape.hat <- optimize(logProfileLik,
lower = shape.lower,
upper = shape.upper,
maximum = TRUE)$maximum
## mleFit <- nlm(minusLogLik,log(c(shape.hat,scale.hat)),hessian=TRUE)
mleFit <- optim(par=log(c(shape.hat,scale.hat)),
fn=minusLogLik,
method="BFGS",
hessian=TRUE)
} else {
## some censored events
## if more than 10 events are uncesored get inital guess from them
## otherwise use all events
if (s.dot >= 10) {
scale.hat <- weighted.mean(yi, si)
s2 <- weighted.mean(yi^2, si) - scale.hat^2
shape.hat <- scale.hat^2 / s2
shape.lower <- floor(shape.hat)
shape.upper <- ceiling(shape.hat)
shape.hat <- optimize(logProfileLik,
lower = shape.lower,
upper = shape.upper,
maximum = TRUE)$maximum
} else {
scale.hat <- weighted.mean(yi, ni)
s2 <- weighted.mean(yi^2, ni) - scale.hat^2
shape.hat <- scale.hat^2 / s2
}
## mleFit <- nlm(minusLogLik,log(c(shape.hat,scale.hat)),hessian=TRUE)
mleFit <- optim(par=log(c(shape.hat,scale.hat)),
fn=minusLogLik,
method="BFGS",
hessian=TRUE)
} ## End of conditional on s.dot == n.dot
## estimate <- exp(mleFit$estimate)
estimate <- exp(mleFit$par)
newVar <- (1/estimate) %o% (1/estimate)
observedI <- mleFit$hessian * newVar
se <- sqrt(diag(solve(observedI)))
## l <- -mleFit$minimum
l <- -mleFit$value
if (scale) {
names(estimate) <- c("shape","scale")
names(se) <- c("shape","scale")
rFct <- function(shape,scale) -minusLogLik(log(c(shape,scale))) - l
} else {
estimate <- c(estimate[1],1/estimate[2])
names(estimate) <- c("shape","rate")
se <- se * c(1,estimate[2]^2)
names(se) <- c("shape","rate")
rFct <- function(shape,rate) -minusLogLik(log(c(shape,1/rate))) - l
}
result <- list(estimate = estimate,
se = se,
logLik = l,
r = rFct,
mll = minusLogLik,
call = match.call()
)
class(result) <- "durationFit"
return(result)
}
llogisMLE <- function(yi,
ni = numeric(length(yi))+1,
si = numeric(length(yi))+1
) {
## check if yi is a spikeTrain object if yes take the "diff"
if (inherits(yi,"spikeTrain")) yi <- diff(yi)
## check that yi elements are positive
if (any(yi <= 0)) stop("yi elements must be positive or null")
## coerce yi to vector
yi <- as.numeric(yi)
## check that ni has the same length as yi
if (!identical(length(yi),length(ni))) stop("yi and ni should have the same length")
## check that the elements of ni are non negative and integer
if (any(ni < 0)) stop("ni elements must be non-negative")
if (!identical(ni, round(ni))) stop("ni should be a vector of positive integers")
if (!identical(length(yi),length(si))) stop("yi and si should have the same length")
## check that the elements of ni are non negative and integer
if (any(si < 0)) stop("si elements must be non-negative")
if (!identical(si, round(si))) stop("si should be a vector of positive integers")
if (any(si > ni)) stop("si elements should not be greater than ni elements")
## Create function returning the opposite of the log likelihood
minusLogLik <- function(p) {
if (missing(p)) {
txt <- paste("This function argument should be a 2 component vector:\n",
" component 1 is the location parameter,\n",
" component 2 is the log of the scale parameter.\n"
)
cat(txt)
} else {
location <- p[1]
scale <- exp(p[2])
-(ifelse(s.dot>0,sum(dllogis(yi[si>0],location,scale,log=TRUE)*si[si>0],0))+
ifelse(c.dot>0,sum(pllogis(yi[ci>0],location,scale,lower.tail=FALSE,log.p=TRUE)*ci[ci>0]),0)
)
}
}
## Get the number of uncensored events
s.dot <- sum(si)
## Get the total number of events
n.dot <- sum(ni)
ci <- ni-si
c.dot <- sum(ci)
lyi <- log(yi)
if (s.dot == n.dot) {
## no censored event
## Get the empirical mean which is the MLE of parameter mu
location.moment <- weighted.mean(lyi, si)
scale.moment <- sqrt(3*(weighted.mean(lyi^2, si) - location.moment^2))/pi
} else {
## some censored events
## if more than 10 events are uncesored get inital guess from them
## otherwise use all events
if (s.dot >= 10) {
location.moment <- weighted.mean(lyi, si)
scale.moment <- sqrt(3*(weighted.mean(lyi^2, si) - location.moment^2))/pi
} else {
location.moment <- weighted.mean(lyi, ni)
scale.moment <- sqrt(3*(weighted.mean(lyi^2, ni) - location.moment^2))/pi
}
}
## mleFit <- nlm(minusLogLik,c(location.moment,log(scale.moment)),hessian=TRUE)
mleFit <- optim(fn=minusLogLik,
par=c(location.moment,log(scale.moment)),
method="BFGS",
hessian=TRUE)
## estimate <- c(mleFit$estimate[1],exp(mleFit$estimate[2]))
estimate <- c(mleFit$par[1],exp(mleFit$par[2]))
newVar <- c(1,1/estimate[2]) %o% c(1,1/estimate[2])
se <- sqrt(diag(solve(mleFit$hessian * newVar)))
## l <- -mleFit$minimum
l <- -mleFit$value
names(estimate) <- c("location","scale")
names(se) <- c("location","scale")
rFct <- function(location,scale) -minusLogLik(c(location,log(scale))) - l
result <- list(estimate = estimate,
se = se,
logLik = l,
r = rFct,
mll = minusLogLik,
call = match.call()
)
class(result) <- "durationFit"
return(result)
}
lnormMLE <- function(yi,
ni = numeric(length(yi))+1,
si = numeric(length(yi))+1
) {
## check if yi is a spikeTrain object if yes take the "diff"
if (inherits(yi,"spikeTrain")) yi <- diff(yi)
## check that yi elements are positive
if (any(yi <= 0)) stop("yi elements must be strictly positive")
## coerce yi to vector
yi <- as.numeric(yi)
## check that ni has the same length as yi
if (!identical(length(yi),length(ni))) stop("yi and ni should have the same length")
## check that the elements of ni are non negative and integer
if (any(ni < 0)) stop("ni elements must be non-negative")
if (!identical(ni, round(ni))) stop("ni should be a vector of positive integers")
if (!identical(length(yi),length(si))) stop("yi and si should have the same length")
## check that the elements of ni are non negative and integer
if (any(si < 0)) stop("si elements must be non-negative")
if (!identical(si, round(si))) stop("si should be a vector of positive integers")
if (any(si > ni)) stop("si elements should not be greater than ni elements")
## Get the number of uncensored events
s.dot <- sum(si)
## Get the total number of events
n.dot <- sum(ni)
ci <- ni-si
c.dot <- sum(ci)
## Create function returning the opposite of the log likelihood
minusLogLik <- function(p) {
if (missing(p)) {
txt <- paste("This function argument should be a 2 component vector:\n",
" component 1 is the meanlog parameter,\n",
" component 2 is the log of the sdlog parameter.\n"
)
cat(txt)
} else {
mean <- p[1]
sd <- exp(p[2])
-(ifelse(s.dot>0,sum(dlnorm(yi[si>0],mean,sd,log=TRUE)*si[si>0],0))+
ifelse(c.dot>0,sum(plnorm(yi[ci>0],mean,sd,lower.tail=FALSE,log.p=TRUE)*ci[ci>0]),0)
)
}
}
## Get the empirical mean which is the MLE of parameter mu
if (s.dot >= 2) {
mu.hat <- weighted.mean(log(yi), si)
s2 <- weighted.mean(log(yi)^2, si) - mu.hat^2
} else {
mu.hat <- weighted.mean(log(yi), ni)
s2 <- weighted.mean(log(yi)^2, ni) - mu.hat^2
}
if (s.dot == n.dot) {
## all events are uncensored
estimate <- c(mu.hat,sqrt(s2))
se <- c(sqrt(s2/n.dot),sqrt(s2/(2*n.dot)))
l <- -minusLogLik(c(estimate[1],log(estimate[2])))
} else {
## numerical fit needed
## mleFit <- nlm(minusLogLik,c(mu.hat,0.5*log(s2)),hessian=TRUE)
mleFit <- optim(fn=minusLogLik,
par=c(mu.hat,0.5*log(s2)),
method="BFGS",
hessian=TRUE)
## estimate <- c(mleFit$estimate[1],exp(mleFit$estimate[2]))
estimate <- c(mleFit$par[1],exp(mleFit$par[2]))
newVar <- c(1,1/estimate[2]) %o% c(1,1/estimate[2])
se <- sqrt(diag(solve(mleFit$hessian * newVar)))
## l <- -mleFit$minimum
l <- -mleFit$value
}
names(estimate) <- names(formals(dlnorm))[2:3]
names(se) <- names(formals(dlnorm))[2:3]
result <- list(estimate = estimate,
se = se,
logLik = l,
r = function(meanlog,sdlog) -minusLogLik(c(meanlog,log(sdlog))) - l,
mll = minusLogLik,
call = match.call()
)
class(result) <- "durationFit"
return(result)
}
rexpMLE <- function(yi,
ni = numeric(length(yi))+1,
si = numeric(length(yi))+1
) {
## check if yi is a spikeTrain object if yes take the "diff"
if (inherits(yi,"spikeTrain")) yi <- diff(yi)
## check that yi elements are positive
if (any(yi <= 0)) stop("yi elements must be positive or null")
## coerce yi to vector
yi <- as.numeric(yi)
## check that ni has the same length as yi
if (!identical(length(yi),length(ni))) stop("yi and ni should have the same length")
## check that the elements of ni are non negative and integer
if (any(ni < 0)) stop("ni elements must be non-negative")
if (!identical(ni, round(ni))) stop("ni should be a vector of positive integers")
if (!identical(length(yi),length(si))) stop("yi and si should have the same length")
## check that the elements of ni are non negative and integer
if (any(si < 0)) stop("si elements must be non-negative")
if (!identical(si, round(si))) stop("si should be a vector of positive integers")
if (any(si > ni)) stop("si elements should not be greater than ni elements")
## Create function returning the opposite of the log likelihood
minusLogLik <- function(p) {
if (missing(p)) {
txt <- paste("This function argument should be a 2 component vector:\n",
" component 1 is the log of the rate parameter,\n",
" component 2 is the log of the rp (refractory period) parameter.\n"
)
cat(txt)
} else {
rate <- exp(p[1])
rp <- exp(p[2])
-(ifelse(s.dot>0,sum(drexp(yi[si>0],rate=rate,rp=rp,log=TRUE)*si[si>0],0))+
ifelse(c.dot>0,sum(prexp(yi[ci>0],rate=rate,
rp=rp,lower.tail=FALSE,log.p=TRUE)*ci[ci>0]),0)
)
}
}
## Get the number of uncensored events
s.dot <- sum(si)
## Get the total number of events
n.dot <- sum(ni)
ci <- ni-si
c.dot <- sum(ci)
if (s.dot == 0) stop("No uncensored events")
rp.hat <- min(yi[si==1]) - .Machine$double.eps
mu.hat <- weighted.mean(yi-rp.hat, si)
rate.hat <- 1/mu.hat
estimate <- c(rate.hat,rp.hat)
se <- c(rate.hat/sqrt(s.dot),NA)
l <- -minusLogLik(log(c(rate.hat,rp.hat)))
names(estimate) <- c("rate","rp")
names(se) <- c("rate","rp")
rFct <- function(rate,rp) -minusLogLik(log(c(rate,rp))) - l
result <- list(estimate = estimate,
se = se,
logLik = l,
r = rFct,
mll = minusLogLik,
call = match.call()
)
class(result) <- "durationFit"
return(result)
}
weibullMLE <- function(yi,
ni = numeric(length(yi))+1,
si = numeric(length(yi))+1,
shape.min = 0.05,
shape.max = 5) {
#####################################################
## define internal function makeProfileWeibullLogLik
#####################################################
makeProfileWeibullLogLik <- function(data) {
## check that data is strictly positive
if ( any(data <= 0) )
stop("data should be a vector of positive numbers")
## if data is somthing else than a vector (eg a matrix)
## coerce it to a vector
dim(data) <- NULL
## get the sample size
nDot <- length(data)
## get the sum of the logs of the variate
sumLog <- sum(log(data))
logLik <- function(shape) {
## check that shape is strictly positive
if (shape <= 0) {
return(-Inf)
} else {
term1 <- nDot * log(shape)
term2 <- (shape - 1) * sumLog
term3 <- - nDot * (1 + log(sum(data^shape)/nDot))
return( term1 + term2 + term3)
}
}
return(logLik)
}
###############################################################
## End of internal function makeProfileWeibullLogLik definition
###############################################################
## check if yi is a spikeTrain object if yes take the "diff"
if (inherits(yi,"spikeTrain")) yi <- diff(yi)
## check that yi elements are positive
if (any(yi <= 0)) stop("yi elements must be strictly positive")
## coerce yi to vector
yi <- as.numeric(yi)
## check that ni has the same length as yi
if (!identical(length(yi),length(ni))) stop("yi and ni should have the same length")
## check that the elements of ni are non negative and integer
if (any(ni < 0)) stop("ni elements must be non-negative")
if (!identical(ni, round(ni))) stop("ni should be a vector of positive integers")
if (!identical(length(yi),length(si))) stop("yi and si should have the same length")
## check that the elements of ni are non negative and integer
if (any(si < 0)) stop("si elements must be non-negative")
if (!identical(si, round(si))) stop("si should be a vector of positive integers")
if (any(si > ni)) stop("si elements should not be greater than ni elements")
## Create function returning the opposite of the log likelihood
minusLogLik <- function(p) {
if (missing(p)) {
txt <- paste("This function argument should be a 2 component vector:\n",
" component 1 is the log of the shape parameter,\n",
" component 2 is the log of the scale parameter.\n"
)
cat(txt)
} else {
shape <- exp(p[1])
scale <- exp(p[2])
-(ifelse(s.dot>0,sum(dweibull(yi[si>0],shape=shape,scale=scale,log=TRUE)*si[si>0],0))+
ifelse(c.dot>0,sum(pweibull(yi[ci>0],shape=shape,
scale=scale,lower.tail=FALSE,log.p=TRUE)*ci[ci>0]),0)
)
}
}
## Get the number of uncensored events
s.dot <- sum(si)
## Get the total number of events
n.dot <- sum(ni)
ci <- ni-si
c.dot <- sum(ci)
if (s.dot == n.dot) {
## no censored event
mu.hat <- weighted.mean(yi, si)
s2 <- weighted.mean(yi^2, si) - mu.hat^2
stat1 <- mu.hat^2/(s2 + mu.hat^2)
shapeFct <- function(shape) (gamma(1 + 1/shape))^2 / gamma(1 + 2/shape) - stat1
myRoot <- uniroot(shapeFct, lower = shape.min, upper = shape.max)$root
shapeInterval <- c(floor(myRoot) - .Machine$double.eps,ceiling(myRoot) + .Machine$double.eps)
data <- rep(yi[si>0],each=ni[si>0])
myLogLik <- makeProfileWeibullLogLik(data)
shape.hat <- optimize(myLogLik, interval = shapeInterval, maximum = TRUE)$maximum
scale.hat <- (sum(data^shape.hat)/s.dot)^(1/shape.hat)
} else {
## some censored events
## if more than 10 events are uncesored get inital guess from them
## otherwise use all events
if (s.dot >= 10) {
mu.hat <- weighted.mean(yi, si)
s2 <- weighted.mean(yi^2, si) - mu.hat^2
stat1 <- mu.hat^2/(s2 + mu.hat^2)
shapeFct <- function(shape) (gamma(1 + 1/shape))^2 / gamma(1 + 2/shape) - stat1
myRoot <- uniroot(shapeFct, lower = shape.min, upper = shape.max)$root
shapeInterval <- c(floor(myRoot) - .Machine$double.eps,ceiling(myRoot) + .Machine$double.eps)
data <- rep(yi[si>0],each=ni[si>0])
myLogLik <- makeProfileWeibullLogLik(data)
shape.hat <- optimize(myLogLik, interval = shapeInterval, maximum = TRUE)$maximum
scale.hat <- (sum(data^shape.hat)/s.dot)^(1/shape.hat)
} else {
mu.hat <- weighted.mean(yi, ni)
s2 <- weighted.mean(yi^2, ni) - mu.hat^2
stat1 <- mu.hat^2/(s2 + mu.hat^2)
shapeFct <- function(shape) (gamma(1 + 1/shape))^2 / gamma(1 + 2/shape) - stat1
myRoot <- uniroot(shapeFct, lower = shape.min, upper = shape.max)$root
shapeInterval <- c(floor(myRoot) - .Machine$double.eps,ceiling(myRoot) + .Machine$double.eps)
data <- rep(yi,each=ni)
myLogLik <- makeProfileWeibullLogLik(data)
shape.hat <- optimize(myLogLik, interval = shapeInterval, maximum = TRUE)$maximum
scale.hat <- (sum(data^shape.hat)/n.dot)^(1/shape.hat)
}
} ## End of conditional on s.dot == n.dot
## mleFit <- nlm(minusLogLik,log(c(shape.hat,scale.hat)),hessian=TRUE)
mleFit <- optim(fn=minusLogLik,
par=log(c(shape.hat,scale.hat)),
method="BFGS",
hessian=TRUE)
## estimate <- exp(mleFit$estimate)
estimate <- exp(mleFit$par)
newVar <- (1/estimate) %o% (1/estimate)
observedI <- mleFit$hessian * newVar
se <- sqrt(diag(solve(observedI)))
## l <- -mleFit$minimum
l <- -mleFit$value
names(estimate) <- c("shape","scale")
names(se) <- c("shape","scale")
rFct <- function(shape,scale) -minusLogLik(log(c(shape,scale))) - l
result <- list(estimate = estimate,
se = se,
logLik = l,
r = rFct,
mll = minusLogLik,
call = match.call()
)
class(result) <- "durationFit"
return(result)
}
is.durationFit <- function(obj) {
if (!("durationFit" %in% class(obj))) return(FALSE)
expectedNames <- c("estimate",
"se",
"logLik",
"r",
"mll",
"call")
if (!all(expectedNames %in% names(obj))) return(FALSE)
TRUE
}
qqDuration <- function(durationFit,
CI=c(0.95,0.99),
type="l",
xlab,
ylab,
main,
sub,
ylim,
dataLwd=2,
ablineCol=2,
...
) {
## Check that durationFit is a durationFit object
if (!is.durationFit(durationFit))
stop(paste(deparse(substitute(durationFit)),
"should be a durationFit object."
)
)
## check CI
if (!is.null(CI)) {
## make sure that CI is at most of length 2 otherwise
## keep the first 2 components
if (length(CI) > 2) CI <- CI[1:2]
## Check that each component of CI is between 0 and 1
if (any(CI>=1 | CI<=0))
stop(paste(deparse(substitute(CI)),
"components should be in (0,1)")
)
} ## End of conditional on !is.null(CI)
## get the uncensored events
## The next code line causes an unimportant warning during
## automatic code check by "R CMD check"
x <- evalq(rep(yi[si>0],each=si[si>0]),envir=environment(durationFit$r))
## get the y values of the QQ plot
y <- sort(x)
## get the theoretical x values
distribution <- switch(deparse(durationFit$call[[1]]),
invgaussMLE = qinvgauss,
gammaMLE = qgamma,
weibullMLE = qweibull,
lnormMLE = qlnorm,
rexpMLE = qrexp,
llogisMLE = qllogis)
if (is.null(distribution))
stop(paste("Unknown distribution in",
deparse(substitute(durationFit))
)
)
sampleSize <- length(y)
pp <- ppoints(sampleSize)
x.theo <- do.call(distribution,
c(list(p=pp),
as.list(durationFit$estimate)
)
)
if (!is.null(CI)) {
## Get the confidence intervals
ci <- sapply(1:sampleSize,
function(idx)
as.vector(sapply(CI,
function(l) {
thePs <- c((1-l)/2,1-(1-l)/2)
thePs <- qbeta(thePs,
idx,
sampleSize-idx+1)
do.call(distribution,c(list(p=thePs),
as.list(durationFit$estimate))
)
}
)
)
)
} ## End of conditional on !is.null(CI)
modelName <- switch(deparse(durationFit$call[[1]]),
invgaussMLE = "invgauss",
gammaMLE = "gamma",
weibullMLE = "weibull",
lnormMLE = "lnorm",
rexpMLE = "rexp",
llogisMLE = "llogis")
sampleName <- deparse(durationFit$call[["yi"]])
if (missing(xlab)) xlab <- paste("Quantiles of",
modelName)
if (missing(ylab)) ylab <- "Sample quantiles"
if (missing(main)) main <- paste("QQ plot of",
sampleName,
"vs fitted",
modelName)
if (missing(sub)) {
sub <- paste(sampleSize,
"uncensored intervals used.")
if (!is.null(CI))
sub <- paste(sub,
"CI at:",
paste(CI,collapse=",")
)
}
if (missing(ylim)) ylim <- c(min(min(x.theo),min(y)),
max(max(x.theo),max(y))
)
plot(x.theo,y,
xlab=xlab,
ylab=ylab,
main=main,
sub=sub,
ylim=ylim,
type="n",
...
)
abline(a=0,b=1,col=ablineCol)
if (!is.null(CI)) apply(ci,1,function(x) lines(x.theo,x,lty=2))
lines(x.theo,y,type=type,lwd=dataLwd)
}
compModels <- function(yi,
ni = numeric(length(yi))+1,
si = numeric(length(yi))+1,
models = c("invgauss","lnorm","gamma","weibull","llogis","rexp"),
type = c("qq","s"),
log = TRUE,
plot=TRUE
) {
if (is.spikeTrain(yi)) yi <- diff(yi)
## Check that the models are known
knownModels <- c("invgauss",
"lnorm",
"gamma",
"weibull",
"llogis",
"rexp")
if (!all(models %in% knownModels))
stop("Some requested models are not implemented.")
theFits <- lapply(models,
function(m) switch(m,
invgauss = invgaussMLE(yi,ni,si),
lnorm = lnormMLE(yi,ni,si),
gamma = gammaMLE(yi,ni,si),
weibull = weibullMLE(yi,ni,si),
llogis = llogisMLE(yi,ni,si),
rexp = rexpMLE(yi,ni,si)
)
)
if (plot) {
myNrow <- (length(models)-1) %/% 2 + 1
layout(matrix(1:(myNrow*2),nrow=myNrow))
oldpar <- par(mar=c(4,3,2,1))
on.exit(par(oldpar))
if (type[1] == "qq") {
sapply(1:length(models),
function(idx) {
if (log) {
qqDuration(theFits[[idx]],
main="",
ylab="",
sub="",
log="xy")
} else {
qqDuration(theFits[[idx]],
main="",
ylab="",
sub="")
}
}
)
} else {
isi <- rep(yi[ni>0],each=ni[ni>0])
status <- unlist(lapply(1:length(ni),
function(idx) {
if (ni[idx]==0) return(numeric())
else {
result <- numeric(ni[idx])
if (si[idx]>0) result[1:si[idx]] <- 1
return(result)
}
}
)
)
kmFit <- survfit(Surv(isi,status) ~1)
sapply(1:length(models),
function(idx) {
plot(kmFit,log=log,mark.time=FALSE,main=models[idx])
survFct <- switch(models[idx],
invgauss = pinvgauss,
lnorm = plnorm,
gamma = pgamma,
weibull = pweibull,
llogis = pllogis,
rexp = prexp)
isi <- sort(isi)
y <- do.call(survFct,c(list(q=isi,lower.tail=FALSE),as.list(theFits[[idx]]$estimate)))
lines(isi,y,col=2)
}
)
}
} ## End of conditional on plot
aic <- sapply(theFits, function(f) -2*f$logLik+4)
ordered <- sort.int(aic,index.return=TRUE)
result <- ordered$x
names(result) <- models[ordered$ix]
result
}
dinvgauss <- function(x,
mu = 1,
sigma2 = 1,
boundary = NULL,
log = FALSE
){
if( any(x <= 0) ) stop("y must contain positive values")
if( any(mu <= 0) ) stop("mu must be positive")
if (all(!is.null(sigma2)))
if( any(sigma2 <= 0) ) stop("sigma2 must be positive")
if (all(!is.null(boundary)))
if( any(boundary <= 0) ) stop("boundary must be positive")
if ( all(!is.null(sigma2)) && all(!is.null(boundary)) )
stop("One of sigma2 or boundary must be specified, not both")
if (all(!is.null(boundary))) {
## We work with the parameterization in term of boundary (ie, sigma2 = 1)
## We convert it in term of mu and sigma2
sigma2 <- (1/boundary)^2
mu <- boundary*mu
}
tmp <- -(x-mu)^2/(2*x*sigma2*mu^2)-(log(2*pi*sigma2)+3*log(x))/2
if(!log) tmp <- exp(tmp)
tmp
}
dllogis <- function(x,
location = 0,
scale = 1,
log = FALSE) {
## Check that x elements are strictly positive
if (any(x <= 0)) stop("x elements must be strictly positive")
if (!log) dlogis(log(x),location,scale)/x
else dlogis(log(x),location,scale,log=TRUE) - log(x)
}
drexp <- function(x,
rate = 10,
rp = 0.005,
log = FALSE) {
## Check that rate is strictly positive
if (rate <= 0) stop("rate should be strictly positive")
## Check that rp is positive or null
if (rp < 0) stop("rp should be positive or null")
x <- x-rp
dexp(x,rate=rate,log=log)
}
hgamma <- function(x,
shape,
rate = 1,
scale = 1/rate,
log = FALSE) {
if (!log) dgamma(x,shape,,scale)/pgamma(x,shape,,scale,lower.tail=FALSE)
else dgamma(x,shape,,scale,log=TRUE) -
pgamma(x,shape,,scale,lower.tail=FALSE,log.p=TRUE)
}
hinvgauss <- function(x,
mu = 1,
sigma2 = 1,
boundary = NULL,
log = FALSE) {
if( any(mu <= 0) ) stop("mu must be positive")
if (all(!is.null(sigma2)))
if( any(sigma2 <= 0) ) stop("sigma2 must be positive")
if (all(!is.null(boundary)))
if( any(boundary <= 0) ) stop("boundary must be positive")
if ( all(!is.null(sigma2)) && all(!is.null(boundary)) )
stop("One of sigma2 or boundary must be specified, not both")
if (all(!is.null(boundary))) {
## We work with the parameterization in term of boundary (ie, sigma2 = 1)
## We convert it in term of mu and sigma2
sigma2 <- (1/boundary)^2
mu <- boundary*mu
}
t <- x/mu
v <- sqrt(x*sigma2)
cutOff <- 1e-12
bigEnough <- pinvgauss(x, mu, sigma2, lower.tail = FALSE) > cutOff
logIntensity <- -( (t[bigEnough]-1)^2/(x[bigEnough]*sigma2) + log(2*sigma2*pi*x[bigEnough]^3) )/2 -
log(1-pnorm((t[bigEnough]-1)/v[bigEnough])-exp(2/(mu*sigma2))*pnorm(-(t[bigEnough]+1)/v[bigEnough])
)
if (sum(bigEnough) < length(x)) {
cutOffQ <- qinvgauss(1 - cutOff, mu, sigma2)
cutOffValue <- dinvgauss(cutOffQ, mu, sigma2, log = TRUE) - log(cutOff)
logIntensity <- c(logIntensity,rep(cutOffValue,length(x) - sum(bigEnough)))
}
if (!log) return(exp(logIntensity))
else return(logIntensity)
}
hllogis <- function(x,
location = 0,
scale = 1,
log = FALSE) {
## Check that x elements are strictly positive
if (any(x <= 0)) stop("x elements must be strictly positive")
if (!log) pllogis(x,location,scale)/(scale*x)
else pllogis(x,location,scale,log.p=TRUE) - log(scale) - log(x)
}
hlnorm <- function(x,
meanlog = 0,
sdlog = 1,
log = FALSE) {
## Check that y elements are strictly positive
if (any(x <= 0))
stop(paste("The elements of",
deparse(substitute(x)),
"should be strictly positive.")
)
if (!log) dlnorm(x,meanlog,sdlog)/plnorm(x,meanlog,sdlog,lower.tail=FALSE)
else dlnorm(x,meanlog,sdlog,log=TRUE)-
plnorm(x,meanlog,sdlog,lower.tail=FALSE,log.p=TRUE)
}
hrexp <- function(x,
rate = 10,
rp = 0.005,
log = FALSE) {
## Check that rate is strictly positive
if (rate <= 0) stop("rate should be strictly positive")
## Check that rp is positive or null
if (rp < 0) stop("rp should be positive or null")
if (!log) ifelse(x >= rp, rate, 0)
else ifelse(x >= rp, log(rate), -Inf)
}
hweibull <- function(x,
shape,
scale = 1,
log = FALSE) {
if (!log) dweibull(x, shape, scale) / pweibull(x,shape,scale,lower.tail=FALSE)
else dweibull(x,shape,scale,log=TRUE) - pweibull(x,shape,scale,lower.tail=FALSE,log.p=TRUE)
}
isiHistFit <- function(spikeTrain,
model,
nbins = 10,
CI=0.95,
...) {
spikeTrainName <- deparse(substitute(spikeTrain))
if (!is.spikeTrain(spikeTrain)) spikeTrain <- as.spikeTrain(spikeTrain)
isi <- diff(spikeTrain)
knownModels <- c("invgauss", "lnorm", "gamma", "weibull", "llogis", "rexp")
model <- model[1]
if (!model %in% knownModels)
stop(paste(deparse(substitute(model)),
"is not implemented.")
)
CI <- CI[1]
if (CI <= 0 | CI >= 1)
stop(paste("A CI of", CI, "does not make sense."))
theFit <- switch(model,
invgauss = invgaussMLE(isi),
lnorm = lnormMLE(isi),
gamma = gammaMLE(isi),
weibull = weibullMLE(isi),
llogis = llogisMLE(isi),
rexp = rexpMLE(isi)
)
pSeq <- (1:(nbins-1))/nbins
theEstimates <- theFit$estimate
qFct <- switch(model,
invgauss = function(p) qinvgauss(p,theEstimates[1],theEstimates[2]),
lnorm = function(p) qlnorm(p,theEstimates[1],theEstimates[2]),
gamma = function(p) qgamma(p,theEstimates[1],theEstimates[2]),
weibull = function(p) qweibull(p,theEstimates[1],theEstimates[2]),
llogis = function(p) qllogis(p,theEstimates[1],theEstimates[2]),
rexp = function(p) qrexp(p,theEstimates[1],theEstimates[2])
)
breaks <- c(0,qFct(pSeq),max(isi)+.Machine$double.eps)
mids <- breaks[-(nbins+1)] + diff(breaks)/2
main <- paste("Isi histogram and fitted",
model,
"distribution for",
spikeTrainName
)
sub <- paste("CI at ",
CI*100,
"%. Sample size: ",
length(isi),".",
sep=""
)
hist(isi,
breaks=breaks,
main=main,
sub=sub,
xlab="isi (s)",
...)
X <- seq(.Machine$double.eps,breaks[nbins+1],length.out=501)
dFct <- switch(model,
invgauss = function(x) dinvgauss(x,theEstimates[1],theEstimates[2]),
lnorm = function(x) dlnorm(x,theEstimates[1],theEstimates[2]),
gamma = function(x) dgamma(x,theEstimates[1],theEstimates[2]),
weibull = function(x) dweibull(x,theEstimates[1],theEstimates[2]),
llogis = function(x) dllogis(x,theEstimates[1],theEstimates[2]),
rexp = function(x) drexp(x,theEstimates[1],theEstimates[2])
)
Y <- dFct(X)
lines(X,Y,col=2)
ci <- qbinom(c((1-CI)/2,1-(1-CI)/2),size=length(isi),prob=1/nbins)
l <- ci[1]/(length(isi)*diff(breaks))
u <- ci[2]/(length(isi)*diff(breaks))
invisible(sapply(1:nbins,function(idx) segments(mids[idx],l[idx],mids[idx],u[idx],col=2)))
}
pinvgauss <- function(q,
mu = 1,
sigma2 = 1,
boundary = NULL,
lower.tail = TRUE,
log.p = FALSE
){
if( any(q <= 0) ) stop("q must contain positive values")
if( any(mu <= 0) ) stop("mu must be positive")
if (all(!is.null(sigma2)))
if( any(sigma2 <= 0) ) stop("sigma2 must be positive")
if (all(!is.null(boundary)))
if( any(boundary <= 0) ) stop("boundary must be positive")
if ( all(!is.null(sigma2)) && all(!is.null(boundary)) )
stop("One of sigma2 or boundary must be specified, not both")
if (all(!is.null(boundary))) {
## We work with the parameterization in term of boundary (ie, sigma2 = 1)
## We convert it in term of mu and sigma2
sigma2 <- (1/boundary)^2
mu <- boundary*mu
}
t <- q/mu
v <- sqrt(q*sigma2)
## Use Eq. 4 of Whitemore GA and Seshadri V (1987)
## The American Statistician 41:280-281
if (lower.tail & !log.p)
return(pnorm((t-1)/v)+exp(2/(mu*sigma2))*pnorm(-(t+1)/v))
if (!lower.tail & !log.p)
return(1 - (pnorm((t-1)/v)+exp(2/(mu*sigma2))*pnorm(-(t+1)/v)))
if (lower.tail & log.p)
return(log(pnorm((t-1)/v)+exp(2/(mu*sigma2))*pnorm(-(t+1)/v)))
if (!lower.tail & log.p)
return(log(1 - (pnorm((t-1)/v)+exp(2/(mu*sigma2))*pnorm(-(t+1)/v))))
}
pllogis <- function(q,
location = 0,
scale = 1,
lower.tail = TRUE,
log.p = FALSE) {
plogis(log(q),location,scale,lower.tail,log.p)
}
prexp <- function(q,
rate = 10,
rp = 0.005,
lower.tail = TRUE,
log.p = FALSE) {
## Check that rate is strictly positive
if (rate <= 0) stop("rate should be strictly positive")
## Check that rp is positive or null
if (rp < 0) stop("rp should be positive or null")
q <- q-rp
pexp(q,rate,lower.tail,log.p)
}
qinvgauss <- function(p,
mu = 1,
sigma2 = 1,
boundary = NULL
){
if( any(p < 0 | p > 1) ) stop("p must lie between 0 and 1")
if( any(mu <= 0) ) stop("mu must be positive")
if (all(!is.null(sigma2)))
if( any(sigma2 <= 0) ) stop("sigma2 must be positive")
if (all(!is.null(boundary)))
if( any(boundary <= 0) ) stop("boundary must be positive")
if ( all(!is.null(sigma2)) && all(!is.null(boundary)) )
stop("One of sigma2 or boundary must be specified, not both")
if (all(!is.null(boundary))) {
## We work with the parameterization in term of boundary (ie, sigma2 = 1)
## We convert it in term of mu and sigma2
sigma2 <- (1/boundary)^2
mu <- boundary*mu
}
len <- max(length(p),length(mu),length(sigma2))
if(length(p) != len) {
if(length(p) == 1) p <- rep(p,len)
else stop("length of p incorrect")
}
if(length(mu) != len) {
if(length(mu) == 1) mu <- rep(mu,len)
else stop("length of m incorrect")
}
if(length(sigma2) != len) {
if(length(sigma2) == 1) sigma2 <- rep(sigma2,len)
else stop("length of sigma2 incorrect")
}
## Use Whitemore and Yalovky (1978, Technometrics, 20:207-208)
## approximation to get starting value for the numerical
## inversion of the cumulative distribution function.
theta <- 1/mu/sigma2
approx <- mu * exp(qnorm(p)*sqrt(1/theta)-0.5/theta)
sapply(1:len, function(idx) {
if (identical(p[idx],0)) return(0)
if (identical(p[idx],1)) return(Inf)
interval <- approx[idx]*c(0.95,1.05)
h <- function(q) pinvgauss(q, mu[idx], sigma2[idx]) - p[idx]
while (h(interval[1])*h(interval[2]) > 0)
interval <- interval*c(0.9,1.1)
uniroot(h,interval)$root
}
)
}
qllogis <- function(p,
location = 0,
scale = 1,
lower.tail = TRUE,
log.p = FALSE) {
exp(qlogis(p,location,scale,lower.tail,log.p))
}
qrexp <- function(p,
rate = 10,
rp = 0.005,
lower.tail = TRUE,
log.p = FALSE) {
## Check that rate is strictly positive
if (rate <= 0) stop("rate should be strictly positive")
## Check that rp is positive or null
if (rp < 0) stop("rp should be positive or null")
rp + qexp(p,rate,lower.tail,log.p)
}
rinvgauss <- function(n = 1,
mu = 1,
sigma2 = 1,
boundary = NULL
){
if( any(mu <= 0) ) stop("mu must be positive")
if (all(!is.null(sigma2)))
if( any(sigma2 <= 0) ) stop("sigma2 must be positive")
if (all(!is.null(boundary)))
if( any(boundary <= 0) ) stop("boundary must be positive")
if ( all(!is.null(sigma2)) && all(!is.null(boundary)) )
stop("One of sigma2 or boundary must be specified, not both")
if (all(!is.null(boundary))) {
## We work with the parameterization in term of boundary (ie, sigma2 = 1)
## We convert it in term of mu and sigma2
sigma2 <- (1/boundary)^2
mu <- boundary*mu
}
## Use method of Michael JR, Schucany WR and Haas RW (1976)
## The American Statistician, 30:88-90
v0 <- rchisq(n,1)
x1 <- mu + 0.5*mu^2*sigma2*v0 - 0.5*mu*sigma2*sqrt(4*mu*v0/sigma2+mu^2*v0^2)
ifelse(rbinom(length(x1),1,mu/(mu+x1)) == 1,x1,mu^2/x1)
}
rllogis <- function(n,
location = 0,
scale = 1) {
exp(rlogis(n,location,scale))
}
rrexp <- function(n,
rate = 10,
rp = 0.005) {
## Check that rate is strictly positive
if (rate <= 0) stop("rate should be strictly positive")
## Check that rp is positive or null
if (rp < 0) stop("rp should be positive or null")
rp+rexp(n,rate)
}
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.