Nothing
####################################################
### Authors: Simone Padoan and Moreno Bevilacqua.
### Emails: simone.padoan@unibocconi.it,
### moreno.bevilacqua@uv.cl
### Institutions: Department of Decision Sciences,
### University Bocconi of Milan and
### Departamento de Estadistica
### Universidad de Valparaiso
### File name: Distribution.r
### Description:
### This file contains a set of procedures
### that define several distributions and
### related functions.
### Last change: 28/03/2013.
####################################################
### Procedures are in alphabetical order.
FitGev <- function(data, method='Nelder-Mead', start, varest=FALSE)
{
param <- NULL
varcov <- NULL
stderr <- NULL
numdata <- length(data)
if(missing(data) || numdata == 0 || !is.numeric(data))
stop('insert a numeric vector of data')
### Starting values:
momest <- MomEst(data, numdata)
location <- momest$location
scale <- momest$scale
shape <- momest$shape
if(shape >= 0) shape <- .1
else shape <- -.1
param <- list(location=location, scale=scale, shape=shape)
parnames <- c('location', 'scale', 'shape')
if(!missing(start))
{
if(!is.list(start))
stop('start need to be named list of starting values')
for(i in 1 : length(start))
{
namestart <- names(start[i])
if(!any(namestart == parnames))
stop('insert location, scale and shape as a named list')
param[namestart] <- start[namestart]
}
if(param$scale < 0)
param$scale <- scale
}
param <- unlist(param)
### fit the log-likelihood of the GEV model:
fitted <- optim(param, GevLogLik, data=data, method=method, numdata=numdata,
control=list(fnscale=-1, reltol=1e-14, maxit=1e8), hessian=varest)
if(varest)
{
varcov <- - try(solve(fitted$hessian), silent = TRUE)
if(!is.matrix(varcov))
{
varcov <- 'none'
stderr <- 'none'
}
else
{
stderr <- diag(varcov)
if(any(stderr < 0))
stderr <- 'none'
else
stderr <- sqrt(stderr)
}
}
return(list(param=fitted$par, varcov=varcov, stderr=stderr))
}
Dist2Dist <- function(data, from='Gev', to='sFrechet', loc=NULL, scale=NULL, shape=NULL)
{
Dist2Dist <- NULL
if(missing(data) || !is.numeric(data))
stop('insert a numeric vector of data')
dimdata <- dim(data)
if(is.null(dimdata))
{
numdata <- length(data)
if(numdata==0)
stop('insert a numeric vector of data')
numcoord <- 1
dimdata <- c(numdata, numcoord)
dim(data) <- dimdata
numparam <- numdata
}
else
{
if(length(dimdata)==2)
{
numdata <- dimdata[1]
numcoord <- dimdata[2]
}
if(length(dimdata)==3)
{
numdata <- dimdata[3]
numcoord <- prod(dimdata[1:2])
}
numparam <- numcoord
}
if(is.null(loc)) loc <- double(numcoord)
else if(!is.numeric(loc) || !length(loc)==numparam)
stop('insert a numeric vector of locations\n')
if(is.null(scale)) scale <- double(numcoord) + 1
else if(!is.numeric(scale) || !length(scale)==numparam)
stop('insert a numeric vector of scales\n')
if(is.null(shape)) shape <- double(numcoord)
else if(!is.numeric(shape) || !length(shape)==numparam)
stop('insert a numeric vector of shapes\n')
# Define the fitting procedure to compute the ML estimates:
Mle <- function(x)
return(as.double(FitGev(x)$param))
# Transform from GEV to other distributions:
if(from=='Gev')
{
type <- switch(to,
Uniform=0,
sFrechet=1,
sGumbel=2,
sWeibull=3,
Gev=4)
if(is.null(type))
stop('insert a valid distribution name\n')
param <- apply(data, 2, Mle) # ML estimates
# Margins transformation:
Dist2Dist <- .C('Dist2Dist', as.double(data), as.double(param[1,]), as.double(param[2,]),
as.double(param[3,]), as.integer(numdata), as.integer(numcoord),
as.double(loc), as.double(scale), as.double(shape), as.integer(type),
res=double(numdata * numcoord), DUP=TRUE, NAOK=TRUE)$res
}
# Transform from standard Frechet to the GEV distribution:
if(from=='sFrechet')
{
type <- switch(to, Gev=5)
if(is.null(type))
stop('insert a valid distribution name\n')
param <- rep(1, 3)
# Margins transformation:
Dist2Dist <- .C('Dist2Dist', as.double(data), as.double(param[1]), as.double(param[2]),
as.double(param[3]), as.integer(numdata), as.integer(numcoord),
as.double(loc), as.double(scale), as.double(shape), as.integer(type),
res=double(numdata * numcoord), DUP=TRUE, NAOK=TRUE)$res
}
# Transform from standard Gumbel to the GEV distribution:
if(from=='sGumbel')
{
type <- switch(to, Gev=6)
if(is.null(type))
stop('insert a valid distribution name\n')
param <- rep(1, 3)
# Margins transformation:
Dist2Dist <- .C('Dist2Dist', as.double(data), as.double(param[1]), as.double(param[2]),
as.double(param[3]), as.integer(numdata), as.integer(numcoord),
as.double(loc), as.double(scale), as.double(shape), as.integer(type),
res=double(numdata * numcoord), DUP=TRUE, NAOK=TRUE)$res
}
if(from=='sWeibull')
{
type <- switch(to, Gev=7)
if(is.null(type))
stop('insert a valid distribution name\n')
param <- rep(1, 3)
# Margins transformation:
Dist2Dist <- .C('Dist2Dist', as.double(data), as.double(param[1]), as.double(param[2]),
as.double(param[3]), as.integer(numdata), as.integer(numcoord),
as.double(loc), as.double(scale), as.double(shape), as.integer(type),
res=double(numdata * numcoord), DUP=TRUE, NAOK=TRUE)$res
}
dim(Dist2Dist) <- dimdata
return(Dist2Dist)
}
GevLogLik <- function(data, numdata, param)
{
### Compute the log-likelihood:
result <- .C('GevLogLik', as.double(data), as.integer(numdata), as.double(param),
res=double(1), DUP=TRUE, NAOK=TRUE)$res
return(result)
}
MomEst <- function(data, n)
{
# Scale estimate:
scale <- sqrt(6 * var(data)) / pi
# Location estimate:
location <- mean(data) - 0.58 * scale
k <- round(.5 * n)
mind <- min(data)
if(mind < 0) data <- data - mind
data <- sort(data)
delta <- log(data[n-0:(k-1)] / data[n-k])
M1 <- sum(delta) / k
M2 <- sum(delta^2) / k
# Shape estimate:
shape <- M1 + 1 - .5 * (M2 / (M2 - M1^2))
return(list(location=location, scale=scale, shape=shape))
}
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.