[Vinecopula-commits] r42 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mi Nov 27 16:18:08 CET 2013
Author: ben_graeler
Date: 2013-11-27 16:18:07 +0100 (Wed, 27 Nov 2013)
New Revision: 42
Modified:
pkg/DESCRIPTION
pkg/R/BiCopSelect.r
pkg/R/RVineCopSelect.r
pkg/man/BiCopSelect.Rd
pkg/man/RVineCopSelect.Rd
pkg/man/RVineStructureSelect.Rd
Log:
- version 1.2-1
- adopted BICs and AICs in BiCopSelect for larger family indicies
- added in docs: "families not listed in familyset might be included for limit case handling"
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-11-20 11:50:59 UTC (rev 41)
+++ pkg/DESCRIPTION 2013-11-27 15:18:07 UTC (rev 42)
@@ -1,8 +1,8 @@
Package: VineCopula
Type: Package
Title: Statistical inference of vine copulas
-Version: 1.2
-Date: 2013-09-03
+Version: 1.2-1
+Date: 2013-11-27
Author: Ulf Schepsmeier, Jakob Stoeber, Eike Christian Brechmann, Benedikt Graeler
Maintainer: Ulf Schepsmeier <schepsmeier at ma.tum.de>
Depends: R (>= 2.11.0), MASS, mvtnorm, igraph
@@ -10,4 +10,3 @@
Description: This package provides functions for statistical inference of vine copulas. It contains tools for bivariate exploratory data analysis, bivariate copula selection and (vine) tree construction. Models can be estimated either sequentially or by joint maximum likelihood estimation. Sampling algorithms and plotting methods are also included. Data is assumed to lie in the unit hypercube (so-called copula data). For C- and D-vines links to the package CDVine are provided.
License: GPL (>= 2)
LazyLoad: yes
-Packaged: 2013-02-07 14:29:57 UTC; Eike
Modified: pkg/R/BiCopSelect.r
===================================================================
--- pkg/R/BiCopSelect.r 2013-11-20 11:50:59 UTC (rev 41)
+++ pkg/R/BiCopSelect.r 2013-11-27 15:18:07 UTC (rev 42)
@@ -1,5 +1,4 @@
-BiCopSelect <- function(u1,u2,familyset=NA,selectioncrit="AIC",indeptest=FALSE,level=0.05,weights=NA)
-{
+BiCopSelect <- function(u1, u2, familyset=NA, selectioncrit="AIC", indeptest=FALSE, level=0.05, weights=NA) {
if(is.null(u1)==TRUE || is.null(u2)==TRUE) stop("u1 and/or u2 are not set or have length zero.")
if(length(u1)!=length(u2)) stop("Lengths of 'u1' and 'u2' do not match.")
if(length(u1)<2) stop("Number of observations has to be at least 2.")
@@ -22,14 +21,17 @@
}else{
- if(!is.na(familyset[1]) && (!any(c(1,2,5,23,24,26:30,33,34,36:40,104,114,204,214) %in% familyset) || !any(c(1:10,13,14,16:20,124,134,224,234) %in% familyset))) stop("'familyset' has to include at least one bivariate copula family for positive and one for negative dependence.")
+ if(!is.na(familyset[1]) && (!any(c(1,2,5,23,24,26:30,33,34,36:40,104,114,204,214) %in% familyset) || !any(c(1:10,13,14,16:20,124,134,224,234) %in% familyset)))
+ stop("'familyset' has to include at least one bivariate copula family for positive and one for negative dependence.")
+ if(is.na(familyset[1]))
+ familyset <- c(1:10,13,14,16:20,23,24,26:30,33,34,36:40,104,114,124,134,204,214,224,234)
emp_tau = fasttau(data1,data2,weights)
- if(indeptest == TRUE){
+ if(indeptest == TRUE) {
out$p.value.indeptest = BiCopIndTest(data1,data2)$p.value
- }else{
+ } else{
out$p.value.indeptest = NA
}
@@ -38,14 +40,13 @@
out$family = 0
out$par = c(0,0)
- }else{
-
+ } else {
start = list()
start[[1]] = sin(pi * emp_tau/2)
start[[2]] = c(sin(emp_tau*pi/2),10)
start[[3]] = start[[13]] = 2*abs(emp_tau)/(1-abs(emp_tau))
start[[4]] = start[[14]] = 1/(1-abs(emp_tau))
- if(5%in% familyset) start[[5]] = Frank.itau.JJ(emp_tau) else start[[5]] = 0
+ if(5 %in% familyset) start[[5]] = Frank.itau.JJ(emp_tau) else start[[5]] = 0
if(any(c(6,16)%in% familyset)) start[[6]] = start[[16]] = Joe.itau.JJ(abs(emp_tau)) else start[[6]] = start[[16]] = 0
start[[7]] = start[[17]] = c(0.5, 1.5)
start[[8]] = start[[18]] = c(1.5, 1.5)
@@ -58,22 +59,22 @@
start[[28]] = start[[38]] = c(-1.5, -1.5)
start[[29]] = start[[39]] = c(-1.5, -0.5)
start[[30]] = start[[40]] = c(-1.5,-0.5)
- #start[[41]] = start[[51]] = ipsA.tau2cpar(emp_tau)
- #start[[61]] = start[[71]] = -ipsA.tau2cpar(emp_tau)
- start[[104]] = start[[204]] = start[[114]] = start[[214]] = c(2,0.5)
- start[[124]] = start[[224]] = start[[134]] = start[[234]] = c(-2,0.5)
+ #start[[41]] = start[[51]] = ipsA.tau2cpar(emp_tau)
+ #start[[61]] = start[[71]] = -ipsA.tau2cpar(emp_tau)
+ start[[104]] = start[[204]] = start[[114]] = start[[214]] = c(2,0.5)
+ start[[124]] = start[[224]] = start[[134]] = start[[234]] = c(-2,0.5)
if(emp_tau < 0){
todo = c(1,2,5,23,24,26:30,33,34,36:40,124,134,224,234)
}else{
todo = c(1:10,13,14,16:20,104,114,204,214)
}
- if(!is.na(familyset[1])) todo = todo[which(todo %in% familyset)]
+ todo = todo[which(todo %in% familyset)]
optiout = list()
if(any(todo == 2)){
- optiout[[2]] = suppressWarnings(BiCopEst(data1,data2,family=2, max.df=30,weights=weights))
+ optiout[[2]] = suppressWarnings(BiCopEst(data1, data2, family=2, max.df=30, weights=weights))
optiout[[2]]$par=c(optiout[[2]]$par,optiout[[2]]$par2)
if(optiout[[2]]$par[2] >= 30){
todo[todo==2] = 1
@@ -285,80 +286,69 @@
optiout[[40]] = list()
}
}
-
-
-
+
+ cat(todo, "\n")
+
for(i in todo[!(todo%in%c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234))]){
optiout[[i]] = MLE_intern(cbind(data1,data2),start[[i]],i,weights=weights)
}
-
- for(i in todo[(todo%in%c(104,114,124,134,204,214,224,234))]){
+
+ for(i in todo[(todo%in%c(104,114,124,134,204,214,224,234))]){
optiout[[i]] = MLE_intern_Tawn(cbind(data1,data2),start[[i]],i)
}
-
if(selectioncrit == "AIC"){
+ AICs = rep(Inf, max(todo))
- AICs = rep(Inf,40)
-
- for(i in todo)
- {
- if(i %in% c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234))
- {
- if(any(is.na(weights))){
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2])))
- }else{
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2]))%*%weights)
- }
+ for(i in todo) {
+ if(i %in% c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234)) {
+ if(any(is.na(weights))) {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2])))
+ } else {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2]))%*%weights)
+ }
AICs[i] = -2*ll + 4
- }
- else
- {
- if(any(is.na(weights))){
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par)))
- }else{
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par))%*%weights)
- }
+ } else {
+ if(any(is.na(weights))){
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par)))
+ } else {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par))%*%weights)
+ }
AICs[i] = -2*ll + 2
}
}
out$family = todo[which.min(AICs[todo])]
- }else if(selectioncrit == "BIC"){
-
- BICs = rep(Inf,40)
-
- for(i in todo)
- {
- if(i %in% c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234))
- {
- if(any(is.na(weights))){
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2])))
- }else{
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2]))%*%weights)
- }
- BICs[i] = -2*ll + 2*log(length(data1))
+ } else {
+ if(selectioncrit == "BIC") {
+ BICs = rep(Inf, max(todo))
+
+ for(i in todo) {
+ if(i %in% c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234)) {
+ if(any(is.na(weights))) {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2])))
+ } else {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par[1], optiout[[i]]$par[2]))%*%weights)
+ }
+ BICs[i] = -2*ll + 2*log(length(data1))
+ } else {
+ if(any(is.na(weights))) {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par)))
+ } else {
+ ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par))%*%weights)
+ }
+ BICs[i] = -2*ll + log(length(data1))
+ }
}
- else
- {
- if(any(is.na(weights))){
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par)))
- }else{
- ll=sum(log(BiCopPDF(data1,data2,i, optiout[[i]]$par))%*%weights)
- }
- BICs[i] = -2*ll + log(length(data1))
- }
+
+ out$family = todo[which.min(BICs[todo])]
}
-
- out$family = todo[which.min(BICs[todo])]
-
}
-
- out$par = optiout[[out$family]]$par
- if(!(out$family%in%c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234)) ) out$par[2] = 0
-
+ out$par = optiout[[out$family]]$par
+ if(!(out$family %in% c(2,7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234)) )
+ out$par[2] = 0
}
}
@@ -366,5 +356,4 @@
out$par=out$par[1]
return(out)
-
}
Modified: pkg/R/RVineCopSelect.r
===================================================================
--- pkg/R/RVineCopSelect.r 2013-11-20 11:50:59 UTC (rev 41)
+++ pkg/R/RVineCopSelect.r 2013-11-27 15:18:07 UTC (rev 42)
@@ -1,23 +1,35 @@
-RVineCopSelect <- function(data,familyset=NA,Matrix,selectioncrit="AIC",indeptest=FALSE,level=0.05,trunclevel=NA){
+RVineCopSelect <- function(data,familyset=NA,Matrix,selectioncrit="AIC",indeptest=FALSE,level=0.05,trunclevel=NA) {
n = dim(data)[2]
N = nrow(data)
- if(dim(Matrix)[1]!=dim(Matrix)[2]) stop("Structure matrix has to be quadratic.")
- if(max(Matrix)>dim(Matrix)[1]) stop("Error in the structure matrix.")
+ if(dim(Matrix)[1]!=dim(Matrix)[2])
+ stop("Structure matrix has to be quadratic.")
+ if(max(Matrix)>dim(Matrix)[1])
+ stop("Error in the structure matrix.")
- if(N<2) stop("Number of observations has to be at least 2.")
- if(n<2) stop("Dimension has to be at least 2.")
- if(any(data>1) || any(data<0)) stop("Data has be in the interval [0,1].")
- if(!is.na(familyset[1])) for(i in 1:length(familyset)) if(!(familyset[i] %in% c(0,1:10,13,14,16:20,23,24,26:30,33,34,36:40,104,114,124,134,204,214,224,234))) stop("Copula family not implemented.")
- if(selectioncrit != "AIC" && selectioncrit != "BIC") stop("Selection criterion not implemented.")
- if(level < 0 & level > 1) stop("Significance level has to be between 0 and 1.")
+ if(N<2)
+ stop("Number of observations has to be at least 2.")
+ if(n<2)
+ stop("Dimension has to be at least 2.")
+ if(any(data>1) || any(data<0))
+ stop("Data has be in the interval [0,1].")
+ if(!is.na(familyset[1])) {
+ if(any(!(familyset %in% c(0,1:10,13,14,16:20,23,24,26:30,33,34,36:40,104,114,124,134,204,214,224,234))))
+ stop("Copula family not implemented.")
+ }
+ if(selectioncrit != "AIC" && selectioncrit != "BIC")
+ stop("Selection criterion not implemented.")
+ if(level < 0 & level > 1)
+ stop("Significance level has to be between 0 and 1.")
- if(is.na(trunclevel)) trunclevel = n
+ if(is.na(trunclevel))
+ trunclevel = n
types = familyset
- if(trunclevel == 0) types = 0
-
+ if(trunclevel == 0)
+ types = 0
+
M = Matrix
Mold = M
@@ -41,7 +53,6 @@
V$direct[n,,] = t(data[,n:1])
for(i in (n-1):1){
-
for(k in n:(i+1)){
m = MaxMat[k,i]
@@ -74,9 +85,10 @@
}
}
- varnames = rep(0,n)
- for(i in 1:n) varnames[i] = paste("V", i, sep="")
-print(Types)
+ varnames = paste("V", 1:n, sep="")
+
+ print(Types)
+
RVM = RVineMatrix(Mold, family=Types, par=Params, par2=Params2, names=varnames)
return(RVM)
Modified: pkg/man/BiCopSelect.Rd
===================================================================
--- pkg/man/BiCopSelect.Rd 2013-11-20 11:50:59 UTC (rev 41)
+++ pkg/man/BiCopSelect.Rd 2013-11-27 15:18:07 UTC (rev 42)
@@ -18,7 +18,7 @@
\item{u1,u2}{Data vectors of equal length with values in [0,1].}
\item{familyset}{Vector of bivariate copula families to select from
(the independence copula MUST NOT be specified in this vector, otherwise it will be selected).
- The vector has to include at least one bivariate copula family that allows for positive and one that allows for negative dependence.
+ The vector has to include at least one bivariate copula family that allows for positive and one that allows for negative dependence. Not listed copula families might be included to better handle limit cases.
If \code{familyset = NA} (default), selection among all possible families is performed.
Coding of bivariate copula families: \cr
\code{1} = Gaussian copula \cr
@@ -66,7 +66,7 @@
(default: \code{indeptest = FALSE}; see \code{\link{BiCopIndTest}}).
The independence copula is chosen if the null hypothesis of independence cannot be rejected.}
\item{level}{Numeric; significance level of the independence test (default: \code{level = 0.05}).}
- \item{weights}{Numerical; weights for each observation (opitional).}
+ \item{weights}{Numerical; weights for each observation (optional).}
}
\value{
Modified: pkg/man/RVineCopSelect.Rd
===================================================================
--- pkg/man/RVineCopSelect.Rd 2013-11-20 11:50:59 UTC (rev 41)
+++ pkg/man/RVineCopSelect.Rd 2013-11-27 15:18:07 UTC (rev 42)
@@ -16,7 +16,7 @@
\arguments{
\item{data}{An N x d data matrix (with uniform margins).}
\item{familyset}{An integer vector of pair-copula families to select from (the independence copula MUST NOT be specified in this vector unless one wants to fit an independence vine!).
- The vector has to include at least one pair-copula family that allows for positive and one that allows for negative dependence.
+ The vector has to include at least one pair-copula family that allows for positive and one that allows for negative dependence. Not listed copula families might be included to better handle limit cases.
If \code{familyset = NA} (default), selection among all possible families is performed.
The coding of pair-copula families is shown below.}
\item{Matrix}{Lower triangular d x d matrix that defines the R-vine tree structure.}
Modified: pkg/man/RVineStructureSelect.Rd
===================================================================
--- pkg/man/RVineStructureSelect.Rd 2013-11-20 11:50:59 UTC (rev 41)
+++ pkg/man/RVineStructureSelect.Rd 2013-11-27 15:18:07 UTC (rev 42)
@@ -17,7 +17,7 @@
\arguments{
\item{data}{An N x d data matrix (with uniform margins).}
\item{familyset}{An integer vector of pair-copula families to select from (the independence copula MUST NOT be specified in this vector unless one wants to fit an independence vine!).
- The vector has to include at least one pair-copula family that allows for positive and one that allows for negative dependence.
+ The vector has to include at least one pair-copula family that allows for positive and one that allows for negative dependence. Not listed copula families might be included to better handle limit cases.
If \code{familyset = NA} (default), selection among all possible families is performed.
Coding of pair-copula families: \cr
\code{1} = Gaussian copula \cr
Mehr Informationen über die Mailingliste Vinecopula-commits