diff -Nru fcopulae-3042.82.1/ChangeLog fcopulae-4021.84/ChangeLog --- fcopulae-3042.82.1/ChangeLog 2014-09-16 14:49:22.000000000 +0000 +++ fcopulae-4021.84/ChangeLog 2022-07-18 13:58:20.000000000 +0000 @@ -2,6 +2,15 @@ ChangeLog Package fCopulae +2021-07-18 smith + + * DESCRIPTION: Updated version number and maintainer + * R/EllipticalCopulae.R method for random number generation from + normal copula calls fMultivar::rnorm2d rather then + fMultivar:::.rnorm2d. This breaks the reproducability of the + output from previous versions. + * inst/obsolete removed + 2014-09-16 setz * ChangeLog, DESCRIPTION: Updated ChangeLog and DESCRIPTION files diff -Nru fcopulae-3042.82.1/debian/changelog fcopulae-4021.84/debian/changelog --- fcopulae-3042.82.1/debian/changelog 2020-05-31 06:04:40.000000000 +0000 +++ fcopulae-4021.84/debian/changelog 2022-07-20 10:39:56.000000000 +0000 @@ -1,8 +1,11 @@ -fcopulae (3042.82.1-1build1) groovy; urgency=medium +fcopulae (4021.84-1) unstable; urgency=medium - * No-change rebuild against r-api-4.0 + * New upstream release + + * debian/control: Set Standards-Version: to current version + * debian/control: Set Build-Depends: to current R version - -- Steve Langasek Sun, 31 May 2020 06:04:40 +0000 + -- Dirk Eddelbuettel Wed, 20 Jul 2022 05:39:56 -0500 fcopulae (3042.82.1-1) unstable; urgency=medium diff -Nru fcopulae-3042.82.1/debian/control fcopulae-4021.84/debian/control --- fcopulae-3042.82.1/debian/control 2020-05-31 06:04:40.000000000 +0000 +++ fcopulae-4021.84/debian/control 2022-07-20 10:39:56.000000000 +0000 @@ -1,10 +1,9 @@ Source: fcopulae Section: gnu-r Priority: optional -Maintainer: Ubuntu Developers -XSBC-Original-Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper (>= 10), r-base-dev (>= 3.6.3), dh-r, r-cran-sn, r-cran-timedate, r-cran-timeseries, r-cran-fbasics (>= 2100.78), r-cran-fmultivar, xvfb, xauth, xfonts-base -Standards-Version: 4.5.0 +Maintainer: Dirk Eddelbuettel +Build-Depends: debhelper (>= 10), r-base-dev (>= 4.2.1), dh-r, r-cran-sn, r-cran-timedate, r-cran-timeseries, r-cran-fbasics (>= 2100.78), r-cran-fmultivar, xvfb, xauth, xfonts-base +Standards-Version: 4.6.1 Vcs-Browser: https://salsa.debian.org/edd/r-cran-copulae Vcs-Git: https://salsa.debian.org/edd/r-cran-fcopulae.git Homepage: https://cran.r-project.org/package=fCopulae diff -Nru fcopulae-3042.82.1/DESCRIPTION fcopulae-4021.84/DESCRIPTION --- fcopulae-3042.82.1/DESCRIPTION 2020-03-07 11:06:14.000000000 +0000 +++ fcopulae-4021.84/DESCRIPTION 2022-07-20 08:20:02.000000000 +0000 @@ -1,11 +1,12 @@ Package: fCopulae Title: Rmetrics - Bivariate Dependence Structures with Copulae -Date: 2017-11-12 -Version: 3042.82.1 +Date: 2022-07-14 +Version: 4021.84 Author: Diethelm Wuertz [aut], - Tobias Setz [cre], - Yohan Chalabi [ctb] -Maintainer: Tobias Setz + Tobias Setz [aut], + Yohan Chalabi [ctb], + Paul Smith [cre] +Maintainer: Paul Smith Description: Provides a collection of functions to manage, to investigate and to analyze bivariate financial returns by Copulae. Included are the families of Archemedean, Elliptical, @@ -13,10 +14,12 @@ Depends: R (>= 2.15.1), timeDate, timeSeries, fBasics, fMultivar Imports: grDevices, graphics, stats Suggests: methods, RUnit, tcltk, mvtnorm, sn -LazyData: yes License: GPL (>= 2) URL: https://www.rmetrics.org -NeedsCompilation: no -Packaged: 2020-03-07 10:25:19 UTC; hornik Repository: CRAN -Date/Publication: 2020-03-07 11:06:14 UTC +Repository/R-Forge/Project: rmetrics +Repository/R-Forge/Revision: 6224 +Repository/R-Forge/DateTimeStamp: 2022-07-19 14:29:47 +Date/Publication: 2022-07-20 08:20:02 UTC +NeedsCompilation: no +Packaged: 2022-07-19 14:52:09 UTC; rforge diff -Nru fcopulae-3042.82.1/inst/obsolete/R/biv-binning.R fcopulae-4021.84/inst/obsolete/R/biv-binning.R --- fcopulae-3042.82.1/inst/obsolete/R/biv-binning.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/biv-binning.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,377 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received A copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: DESCRIPTION: -# squareBinning Square binning of irregularly spaced points -# plot S3 Method for plotting square binned points -# FUNCTION: DESCRIPTION: -# hexBinning Hexagonal binning of irregularly spaced points -# plot S3 Method for plotting hexagonal binned points -################################################################################ - - -################################################################################ -# FUNCTION: DESCRIPTION: -# squareBinning Square binning of irregularly spaced points -# plot S3 Method for plotting square binned points - - -squareBinning = -function(x, y = NULL, bins = 30) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Returns 2D Histogram Counts - - # Arguments: - # x, y - two vectors of coordinates of data. If y is NULL then x - # is assumed to be a two column matrix, where the first column - # contains the x data, and the second column the y data. - # 'timeSeries' objects are also allowed as input. - # bins - number of bins in each dimension, may be a scalar or a 2 - # element vector. The default value is 20. - - # Value: - # A list with three elements x, y, and z. x and y are vectors - # spanning the two dimensioanl grid and z the corresponding - # matrix. The output can directly serve as input to the - # plotting functions image, contour and persp. - - # Example: - # sB = squareBinning(x = rnorm(1000), y = rnorm(1000)); plot(sB) - - # Note: - # Partly copied from R Package gregmisc, function 'hist2d'. - - # FUNCTION: - - # 2D Histogram Counts: - if (is.null(y)) { - x = as.matrix(x) - y = x[, 2] - x = x[, 1] - } else { - x = as.vector(x) - y = as.vector(y) - } - data = cbind(x, y) - - # Bins: - n = bins - if (length(n) == 1) { - nbins = c(n, n) - } else { - nbins = n - } - - # Binning: - xo = seq(min(x), max(x), length = nbins[1]) - yo = seq(min(y), max(y), length = nbins[2]) - xvals = xo[-1] - diff(xo)/2 - yvals = yo[-1] - diff(yo)/2 - ix = findInterval(x, xo) - iy = findInterval(y, yo) - xcm = ycm = zvals = matrix(0, nrow = nbins[1], ncol = nbins[2]) - - for (i in 1:length(x)) { - zvals[ix[i], iy[i]] = zvals[ix[i], iy[i]] + 1 - xcm[ix[i], iy[i]] = xcm[ix[i], iy[i]] + x[i] - ycm[ix[i], iy[i]] = ycm[ix[i], iy[i]] + y[i] - } - - # Reduce to non-empty cells: - u = v = w = ucm = vcm = rep(0, times = nbins[1]*nbins[2]) - L = 0 - for (i in 1:(nbins[1]-1)) { - for (j in 1:(nbins[2]-1)) { - if (zvals[i, j] > 0) { - L = L + 1 - u[L] = xvals[i] - v[L] = yvals[j] - w[L] = zvals[i, j] - ucm[L] = xcm[i, j]/w[L] - vcm[L] = ycm[i, j]/w[L] - } - } - } - length(u) = length(v) = length(w) = L - length(ucm) = length(vcm) = L - - ans = list(x = u, y = v, z = w, xcm = ucm, ycm = vcm, bins = bins, - data = data) - class(ans) = "squareBinning" - - # Return Value: - ans - -} - - -# ------------------------------------------------------------------------------ - - -plot.squareBinning = -function(x, col = heat.colors(12), addPoints = TRUE, addRug = TRUE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Plot square binned data points - - # FUNCTION: - - # Binning: - X = x$x - Y = x$y - - # Plot Center Points: - plot(X, Y, type = "n", ...) - - # Create Hexagon Coordinates: - rx = min(diff(unique(sort(X))))/2 - ry = min(diff(unique(sort(Y))))/2 - u = c(-rx, rx, rx, -rx) - v = c( ry, ry, -ry, -ry) - - # Create Color Palette: - N = length(col) - Z = x$z - zMin = min(Z) - zMax = max(Z) - Z = (Z - zMin)/(zMax - zMin) - Z = trunc(Z*(N-1)+1) - - # Add Colored Hexagon Polygons: - for (i in 1:length(X)) { - polygon(u+X[i], v+Y[i], col = col[Z[i]], border = "white") - } - - # Add Center of Mass Points: - if (addPoints) { - points(x$xcm, x$ycm, pch = 19, cex = 1/3, col = "black") - } - - # Add rug: - if (addRug) { - rug(x$data[, 1], ticksize = 0.01, side = 3) - rug(x$data[, 2], ticksize = 0.01, side = 4) - } - - - # Return Value: - invisible(NULL) -} - - -################################################################################ -# FUNCTION: DESCRIPTION: -# hexBinning Hexagonal binning of irregularly spaced points -# plot S3 Method for plotting hexagonal binned points - - -hexBinning = -function(x, y = NULL, bins = 30) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Does a hexagonal binning of data points - - # Arguments: - # x, y - two vectors of coordinates of data. If y is NULL then x - # is assumed to be a two column matrix, where the first column - # contains the x data, and the second column the y data. - # 'timeSeries' objects are also allowed as input. - # bins - number of bins in each dimension, may be a scalar or a 2 - # element vector. The default value is 20. - - # Example: - # hB = hexBinning(x = rnorm(10000), y = rnorm(10000)); plot(hB) - - # FUNCTION: - - # Extract Series: - if (is.null(y)) { - x = as.matrix(x) - y = x[, 2] - x = x[, 1] - } else { - x = as.vector(x) - y = as.vector(y) - } - data = cbind(x, y) - - # Set Parameters: - shape = 1 - n = length(x) - xbnds = range(x) - ybnds = range(y) - jmax = floor(bins + 1.5001) - c1 = 2 * floor((bins *shape)/sqrt(3) + 1.5001) - imax = trunc((jmax*c1 -1)/jmax + 1) - lmax = jmax * imax - cell = cnt = xcm = ycm = rep(0, times = max(n, lmax)) - xmin = xbnds[1] - ymin = ybnds[1] - xr = xbnds[2] - xmin - yr = ybnds[2] - ymin - c1 = bins/xr - c2 = bins*shape/(yr*sqrt(3.0)) - jinc = jmax - lat = jinc + 1 - iinc = 2*jinc - con1 = 0.25 - con2 = 1.0/3.0 - - # Count Bins: - for ( i in 1:n ) { - sx = c1 * (x[i] - xmin) - sy = c2 * (y[i] - ymin) - j1 = floor(sx + 0.5) - i1 = floor(sy + 0.5) - dist1 = (sx-j1)^2 + 3.0*(sy-i1)^2 - if( dist1 < con1) { - L = i1*iinc + j1 + 1 - } else if (dist1 > con2) { - L = floor(sy)*iinc + floor(sx) + lat - } else { - j2 = floor(sx) - i2 = floor(sy) - test = (sx-j2 -0.5)^2 + 3.0*(sy-i2-0.5)^2 - if ( dist1 <= test ) { - L = i1*iinc + j1 + 1 - } else { - L = i2*iinc + j2 + lat - } - } - cnt[L] = cnt[L]+1 - xcm[L] = xcm[L] + (x[i] - xcm[L])/cnt[L] - ycm[L] = ycm[L] + (y[i] - ycm[L])/cnt[L] - } - - # Reduce to Non-Empty Cells: - nc = 0 - for ( L in 1:lmax ) { - if(cnt[L] > 0) { - nc = nc + 1 - cell[nc] = L - cnt[nc] = cnt[L] - xcm[nc] = xcm[L] - ycm[nc] = ycm[L] - } - } - bnd = c(imax, jmax) - bnd[1] = (cell[nc]-1)/bnd[2] + 1 - length(cell) = nc - length(cnt) = nc - length(xcm) = nc - length(ycm) = nc - if(sum(cnt) != n) warning("Lost counts in binning") - - # Compute Positions: - c3 = diff(xbnds)/bins - ybnds = ybnds - c4 = (diff(ybnds) * sqrt(3))/(2 * shape * bins) - cell = cell - 1 - i = cell %/% jmax - j = cell %% jmax - y = c4 * i + ybnds[1] - x = c3 * ifelse(i %% 2 == 0, j, j + 0.5) + xbnds[1] - - # Result: - ans = list(x = x, y = y, z = cnt, xcm = xcm, ycm = ycm, bins = bins, - data = data) - class(ans) = "hexBinning" - - # Return Value: - ans - -} - - -# ------------------------------------------------------------------------------ - - -plot.hexBinning = -function(x, col = heat.colors(12), addPoints = TRUE, addRug = TRUE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Plot hexagonal binned data points - - # Example: - # hexPlot(rnorm(1000), rnorm(1000), bins = 20) - - # FUNCTION: - - # Binning: - X = x$x - Y = x$y - - # Plot Center Points: - plot(X, Y, type = "n", ...) - - # Create Hexagon Coordinates: - rx = min(diff(unique(sort(X)))) - ry = min(diff(unique(sort(Y)))) - rt = 2*ry - u = c(rx, 0, -rx, -rx, 0, rx) - v = c(ry, rt, ry, -ry, -rt, -ry) / 3 - - # Create Color Palette: - N = length(col) - z = x$z - zMin = min(z) - zMax = max(z) - Z = (z - zMin)/(zMax - zMin) - Z = trunc(Z*(N-1)+1) - - # Add Colored Hexagon Polygons: - for (i in 1:length(X)) { - polygon(u+X[i], v+Y[i], col = col[Z[i]], border = "white") - } - - # Add Center of Mass Points: - if (addPoints) { - points(x$xcm, x$ycm, pch = 19, cex = 1/3, col = "black") - } - - # Add rug: - if (addRug) { - rug(x$data[, 1], ticksize = 0.01, side = 3) - rug(x$data[, 2], ticksize = 0.01, side = 4) - } - - # Return Value: - invisible(NULL) -} - - -################################################################################ - - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/biv-density.R fcopulae-4021.84/inst/obsolete/R/biv-density.R --- fcopulae-3042.82.1/inst/obsolete/R/biv-density.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/biv-density.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,251 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: DESCRIPTION: -# grid2d Returns from two vectors x-y grid coordinates -# density2d Returns 2D Kernel Density Estimates -# hist2d Returns 2D Histogram Counts -# integrate2d Integrates over a two dimensional unit square -################################################################################ - - -grid2d = -function(x = (0:10)/10, y = x) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Creates from two vectors x-y grid coordinates - - # Arguments: - # x, y - two numeric vectors defining the x and y coordinates. - - # Value: - # returns a list with two vectors named $x and $y spanning the - # grid defined by the coordinates x and y. - - # Example: - # > grid2d(1:3, 1:2) - # $x - # [1] 1 2 3 1 2 3 - # $y - # [1] 1 1 1 2 2 2 - - # FUNCTION: - - # Prepare for Input: - nx = length(x) - ny = length(y) - xoy = cbind(rep(x, ny), as.vector(matrix(y, nx, ny, byrow = TRUE))) - XY = matrix(xoy, nx * ny, 2, byrow = FALSE) - - # Return Value: - list(x = XY[, 1], y = XY[, 2]) -} - - -# ------------------------------------------------------------------------------ - - -density2d = -function (x, y = NULL, n = 20, h = NULL, limits = c(range(x), range(y))) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Returns 2D Kernel Density Estimates - - # Arguments: - # x, y - two vectors of coordinates of data. If y is NULL then x - # is assumed to be a two column matrix, where the first column - # contains the x data, and the second column the y data. - # n - Number of grid points in each direction. - # h - a vector of bandwidths for x and y directions. Defaults to - # normal reference bandwidth. - # limits - the limits of the rectangle covered by the grid. - - # Value: - # A list with three elements x, y, and z. x and y are vectors - # spanning the two dimensioanl grid and z the corresponding - # matrix. The output can directly serve as input to the - # plotting functions image, contour and persp. - - # Details: - # Two-dimensional kernel density estimation with an axis-aligned - # bivariate normal kernel, evaluated on a square grid. - - # Note: - # Partly copied from R Package MASS, function 'kde2d'. - - # Reference: - # Venables, W.N., Ripley, B. D. (2002); - # Modern Applied Statistics with S. - # Fourth edition, Springer. - - # FUNCTION: - - # Settings: - lims = limits - if (is.null(y)) { - y = x[, 2] - x = x[, 1] - } - - # Bandwidth: - .bandwidth.nrd = function (x) { - r = quantile(x, c(0.25, 0.75)) - h = (r[2] - r[1])/1.34 - 4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5) } - - # Kernel Density Estimator: - nx = length(x) - if (length(y) != nx) stop("Data vectors must be the same length") - gx = seq(lims[1], lims[2], length = n) - gy = seq(lims[3], lims[4], length = n) - if (is.null(h)) h = c(.bandwidth.nrd(x), .bandwidth.nrd(y)) - h = h/4 - ax = outer(gx, x, "-")/h[1] - ay = outer(gy, y, "-")/h[2] - z = matrix(dnorm(ax), n, nx) %*% t(matrix(dnorm(ay), n, - nx))/(nx * h[1] * h[2]) - - # Return Value: - list(x = gx, y = gy, z = z) -} - - -# ------------------------------------------------------------------------------ - - -hist2d = -function(x, y = NULL, n = c(20, 20)) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Returns 2D Histogram Counts - - # Arguments: - # x, y - two vectors of coordinates of data. If y is NULL then x - # is assumed to be a two column matrix, where the first column - # contains the x data, and the second column the y data. - # n - number of bins in each dimension, may be a scalar or a 2 - # element vector. The default value is 20. - - # Value: - # A list with three elements x, y, and z. x and y are vectors - # spanning the two dimensioanl grid and z the corresponding - # matrix. The output can directly serve as input to the - # plotting functions image, contour and persp. - - # Note: - # Partly copied from R Package gregmisc, function 'hist2d'. - - # FUNCTION: - - # 2D Histogram Counts: - if (is.null(y)) { - y = x[, 2] - x = x[, 1] - } - if (length(n) == 1) { - nbins = c(n, n) - } else { - nbins = n - } - nas = is.na(x) | is.na(y) - x.cuts = seq(from = min(x, y), to = max(x,y), length = nbins[1]+1) - y.cuts = seq(from = min(x, y), to = max(x,y), length = nbins[2]+1) - index.x = cut(x, x.cuts, include.lowest = TRUE) - index.y = cut(y, y.cuts, include.lowest = TRUE) - m = matrix(0, nrow=nbins[1], ncol = nbins[2], - dimnames = list( levels(index.x), levels(index.y) ) ) - for ( i in 1:length(index.x) ) { - m[index.x[i], index.y[i] ] = m[index.x[i], index.y[i] ] + 1 - } - xvals = x.cuts[1:nbins[1]] - yvals = y.cuts[1:nbins[2]] - - # Return Value: - list(x = xvals, y = yvals, z = m) -} - - -# ------------------------------------------------------------------------------ - - -integrate2d = function(fun, error = 1.0e-5, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # 2-dimension quadrature rule on [0,1]^2 - - # Arguments: - # fun - function to be integrated. The first argument requests - # the x values, the second the y values, and the remaining - # are reserved for additional parameters. - # ... - parameters passed to the function to be integrated - - # Details: - # see: Abramowitz and Stegun, p. 892 - - # FUNCTION: - - # Estimate a reasonable number of subintervals: - H = sqrt(sqrt(error)) - n = ceiling(1/H + 1) - blocks = ceiling(log(n+1)/log(2)) - n = 2^blocks-1 - h = 1/(n-1) - - # The error will be of order h^4: - error = h^4 - - # Create all grid coordinates: - x = y = h*seq(1, n-1, by = 2) - nx = ny = length(x) - xoy = cbind(rep(x, ny), as.vector(matrix(y, nx, ny, byrow = TRUE))) - XY = matrix(xoy, nx * ny, 2, byrow = FALSE) - - # The integration rule: - rule = function(x, h, ...) { - X = x[1] + h*c( 0, -1, -1, 1, 1, -1, 1, 0, 0) - Y = x[2] + h*c( 0, -1, 1, -1, 1, 0, 0, -1, 1) - W = c( 16, 1, 1, 1, 1, 4, 4, 4, 4)/36 - ans = sum( W * fun(X, Y, ...) ) - } - - # Result: - ans = (4*h^2)*sum(apply(XY, 1, rule, h = h, ...)) - - # Return Value: - list(value = ans, error = error, points = n) -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/biv-gridding.R fcopulae-4021.84/inst/obsolete/R/biv-gridding.R --- fcopulae-3042.82.1/inst/obsolete/R/biv-gridding.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/biv-gridding.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,133 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received A copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -# fEcofin::4A-BivariateGridding.R -################################################################################ -# FUNCTION: GRID DATA: -# gridData Generates grid data set -# persp.gridData Generates perspective plot from a grid data object -# contour.gridData Generates contour plot from a grid data object -################################################################################ - - -################################################################################ -# FUNCTION: GRID DATA: -# gridData Generates grid data set -# persp.gridData Generates perspective plot from a grid data object -# contour.gridData Generates contour plot from a grid data object - - -gridData = -function(x = (-10:10)/10, y = x, z = outer(x, y, function(x, y) (x^2+y^2)) ) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Generates a grid data set - - # Arguments: - # x, y - two numeric vectors of grid pounts - # z - a numeric matrix or any other rectangular object which can - # be transformed by the function 'as.matrix' into a matrix - # object. - - # Example: - # persp(as.gridData()) - - # FUNCTION: - - # Grid Data: - data = list(x = x, y = y, z = as.matrix(z)) - class(data) = "gridData" - - # Return Value: - data -} - - -# ------------------------------------------------------------------------------ - - -persp.gridData = -function(x, theta = -40, phi = 30, col = "steelblue", ticktype = "detailed", -...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # S3 method to generate a perspective plot from a grid data object - - # Example: - # x = y = seq(-10, 10, length = 30) - # z = outer(x, y, function(x, y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }) - # data = list(x = x, y = y, z = z) - # class(data) = "gridData" - # persp(data) - - # FUNCTION: - - # Grid Data: - class(x) = "default" - persp(x, theta = theta, phi = phi, col = col, ticktype = ticktype, ...) - - # Return Value: - invisible(NULL) -} - - -# ------------------------------------------------------------------------------ - - -contour.gridData = -function(x, addImage = TRUE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # S3 method to generate a contour plot from a grid data object - - # Example: - # x = y = seq(-10, 10, length = 30) - # z = outer(x, y, function(x, y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }) - # data = list(x = x, y = y, z = z) - # class(data) = "gridData" - # contour(data) - - # FUNCTION: - - # Grid Data: - class(x) = "default" - if (addImage) image(x, ...) - contour(x, add = addImage, ...) - box() - - # Return Value: - invisible(NULL) -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/builtin-adapt.R fcopulae-4021.84/inst/obsolete/R/builtin-adapt.R --- fcopulae-3042.82.1/inst/obsolete/R/builtin-adapt.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/builtin-adapt.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,162 +0,0 @@ - - -# Title: adapt -- multidimensional numerical integration -# Package: adapt -# Version: 1.0-4 -# Author: FORTRAN by Alan Genz, -# S by Mike Meyer, R by Thomas Lumley and Martin Maechler -# Description: Adaptive Quadrature in up to 20 dimensions -# Depends: -# License: Unclear (Fortran) -- code in Statlib's ./S/adapt -# Maintainer: Thomas Lumley -# Packaged: Fri Apr 20 11:38:07 2007; thomas - - -# [from Statlib's original http://lib.stat.cmu.edu/S/adapt ] -# This code contains an S function and supporting C and Fortran code for -# adaptive quadrature. The underlyling fortran code is purported to -# work in from 2 to 20 dimensions. The code is set up to dynamically -# load from a central library area. If you can not do dynamic loading, -# you may need to build a staticly loaded version. The adapt S function -# calls load.if.needed to do the dynamic loading. You will have to -# change the functions used here (probably to call library.dynam). -# S code written by Michael Meyer (mikem@andrew.cmu.edu). -# October, 1989. - - -# 2002-03-14 Martin Maechler -# * DESCRIPTION (Version): 1.0-3 --> CRAN -# * R/adapt.R (adapt): use defaults for minpts, maxpts, eps; -# more logical maxpts default (for ndim >= 7) using rulcls -# * man/adapt.Rd: extended example -# 2002-03-13 Martin Maechler -# * DESCRIPTION (Version): 1.0-2 -# * man/adapt.Rd: indentation, using \code{.}, etc; -# example also tries p=5 dimensions -# * R/adapt.R: clean up (spaces) -# 2002-01-09 Martin Maechler -# * R/adapt.R: do not use .Alias anymore -# 2001-06-29 Thomas Lumley -# * move (improved!) integrate() into base, using .Call() etc. - - -# Message-ID: <4AD7A74B.3020108@math.wsu.edu> -# Date: Thu, 15 Oct 2009 15:50:51 -0700 -# From: Alan Genz -# User-Agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.8.1.21) -# Gecko/20090402 SeaMonkey/1.1.16 -# MIME-Version: 1.0 -# To: Diethelm Wuertz -# CC: Alan C Genz -# Subject: Re: adapt -# References: <4AD3032B.4090801@itp.phys.ethz.ch> -# In-Reply-To: <4AD3032B.4090801@itp.phys.ethz.ch> -# Content-Type: text/plain; charset=ISO-8859-1; format=flowed -# Content-Transfer-Encoding: 7bit -# Status: O -# Dear Prof. Wuertz, -# Thank you for your message and your interest in my adaptive integration -# Fortran code. I would be pleased if you included my code in your open -# source R fCopulae package under the Gnu GPL2 license. You have my -# permission to do this. -# Sincerely, -# Alan Genz - - -################################################################################ - - -adapt <- function (ndim, lower, upper, minpts = 100, maxpts = NULL, - functn, eps = 0.01, ...) -{ - keep.trying <- is.null(maxpts) - - if (ndim == 1) { ## fudge for 1-d functions - warning("Using integrate() from base package for 1-d integration") - if (keep.trying) maxpts <- minpts - return(integrate(functn,lower,upper,subdivisions=maxpts,rel.tol=eps,...)) - } - ## else ndim >= 2 : - - ## Check to make sure that upper and lower are reasonable lengths - ## Both the upper and lower limits should be at least of length ndim - if (length(lower) < ndim || length(upper) < ndim)#MM: dropped 'at least': - stop(paste("The lower and upper vectors need to have ndim elements\n", - "Your parameters are: ndim", ndim, ", length(lower)", - length(lower), ", length(upper)", length(upper), "\n")) - ff <- - if(length(list(...)) && length(formals(functn)) > 1) - function(x) functn(x, ...) - else functn # .Alias - rulcls <- 2^ndim + 2*ndim^2 + 6*ndim + 1 #-> ../src/adapt.f - - ## maxpts should be large enough. Prefer 10*rulclc, but use 2*rulclc. - if (keep.trying) - maxpts <- max(minpts, 500, 2 * rulcls) - else { - if (minpts >= maxpts) { - warning(paste("maxpts must be > minpts.\n", - "Maxpts has be increased to minpts + 1")) - maxpts <- minpts + 1 - } - ## - if (maxpts < 2 * rulcls) { - warning(paste( - "You have maxpts (= ", maxpts, ") too small\n", - "It needs to be at least 2 times 2^ndim + 2*ndim^2 + 6*ndim+1\n", - "It has been reset to ", 2 * rulcls, "\n", sep="")) - maxpts <- 2 * rulcls - } - } - - repeat { - lenwrk <- (2*ndim + 3)* (1 + maxpts/rulcls)/2# mandated in adapt source - - x <- .C("cadapt", - as.integer(ndim), - as.double(lower), - as.double(upper), - minpts = as.integer(minpts), - maxpts = as.integer(maxpts), - ## now pass ff and current environment - ff, rho = environment(), - as.double(eps), - relerr = double(1), - lenwrk = as.integer(lenwrk), - value = double(1), # will contain the value of the integral - ifail = integer(1), - PACKAGE = "fCopulae")[ - c("value", "relerr", "minpts", "lenwrk", "ifail")] - - if (x$ifail == 1 && keep.trying) - maxpts <- maxpts*2 - else - break - } - if(x$ifail) - warning(x$warn <- - c("Ifail=1, maxpts was too small. Check the returned relerr!", - paste("Ifail=2, lenwrk was too small. -- fix adapt() !\n", - "Check the returned relerr!"), - "Ifail=3: ndim > 20 -- rewrite the fortran code ;-) !", - "Ifail=4, minpts > maxpts; should not happen!", - "Ifail=5, internal non-convergence; should not happen!" - )[x$ifail]) - - class(x) <- "integration" - - x -} - - -# ------------------------------------------------------------------------------ - - -print.integration <- function(x, ...) { - print(noquote(sapply(x, format, ...)),...) - invisible(x) -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/bv-dcauchy.R fcopulae-4021.84/inst/obsolete/R/bv-dcauchy.R --- fcopulae-3042.82.1/inst/obsolete/R/bv-dcauchy.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/bv-dcauchy.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,128 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: BIVARIATE CAUCHY DISTRIBUTION: -# pcauchy2d Computes bivariate Cauchy probability function -# dcauchy2d Computes bivariate Cauchy density function -# rcauchy2d Generates bivariate Cauchy random deviates -################################################################################ - - -pcauchy2d = -function(x, y = x, rho = 0) -{ # A function Implemented by Diethelm Wuertz - - # Description: - # Computes bivariate Cauchy probability function - - # Arguments: - # x, y - two numeric values or vectors of the same length at - # which the probability will be computed. - - # Example: - # pt2d(rnorm(5), rnorm(5), 0.5, 5) - - # Value: - # returns a numeric vector of probabilities of the same length - # as the input vectors - - # FUNCTION: - - # Settings: - # Probaility: - ans = pt2d(x = x, y = y, rho = rho, nu = 1) - attr(ans, "control") = c(rho = rho) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -dcauchy2d = -function(x, y = x, rho = 0) -{ # A function implemented by Diethelm Wuertz - - # Arguments: - # n - number of random deviates to be generated - # rho - the linear correlation, a numeric value between - # minus one and one. - - # Description: - # Computes bivariate Cauchy density function - - # Note: - # Partly copied from contributed R package 'sn' - - # FUNCTION: - - # Density: - density = dt2d(x = x, y = y, rho = rho, nu = 1) - attr(density, "control") = c(rho = rho) - - # Return value: - density -} - - -# ------------------------------------------------------------------------------ - - -rcauchy2d = -function(n, rho = 0) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Generates bivariate Cauchy random deviates - - # Arguments: - # n - number of random deviates to be generated - # rho - the linear correlation, a numeric value between - # minus one and one. - - # Note: - # Partly copied from contributed R package 'mvtnorm' - # Author Friedrich Leisch - - # FUNCTION: - - # Random Deviates: - ans = rt2d(n = n, rho = rho) - attr(ans, "control") = c(rho = rho) - - # Return Value: - ans -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/bv-delliptical.R fcopulae-4021.84/inst/obsolete/R/bv-delliptical.R --- fcopulae-3042.82.1/inst/obsolete/R/bv-delliptical.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/bv-delliptical.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,313 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: ELLIPTICAL BIVARIATE DISTRIBUTIONS: -# delliptical2d Computes density for elliptical distributions -# .gfunc2d Generator Function for elliptical distributions -# .delliptical2dSlider Slider for bivariate densities -################################################################################ - - -delliptical2d = -function(x, y = x, rho = 0, param = NULL, type = c("norm", "cauchy", "t", -"logistic", "laplace", "kotz", "epower"), output = c("vector", "list")) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Density function for bivariate elliptical distributions - - # Arguments: - # x, y - two numeric vectors of the same length. - # rho - a anumeric value specifying the correlation. - # param - NULL, a numeric value, or a numeric vector adding - # additional parameters to the generator function. - # type - a character string denoting the type of distribution. - # This may be either - # "norm" for the normal distribution, or - # "cauchy" for the Cauchy distribution, or - # "t" for the Student-t distribution, or - # "logistic" for the logistic distribution, or - # "laplace" for the distribution, or - # "kotz" for the original Kotz distribution, or - # "epower" for the exponential power distribution - - # FUNCTION: - - # Type: - type = type[1] - - # Settings: - if (is.list(x)) { - y = x$y - x = x$x - } - if (is.matrix(x)) { - y = x[, 2] - x = x[, 2] - } - - # Add Default Parameters: - if (is.null(param)) { - if (type == "t") param = c(nu = 4) - if (type == "kotz") param = c(r = sqrt(2)) - if (type == "epower") param = c(r = sqrt(2), s = 1/2) - } - - # Density: - xoy = ( x^2 - 2*rho*x*y + y^2 ) / (1-rho^2) - lambda = .gfunc2d(param = param, type = type)[[1]] - density = lambda * .gfunc2d(x = xoy, param = param, type = type) / - sqrt(1 - rho^2) - - # Add attributes: - if (is.null(param)) { - attr(density, "control") = unlist(list(type = type, rho = rho)) - } else { - attr(density, "control") = unlist(list(type = type, rho = rho, - param = param)) - } - - # As List ? - if (output[1] == "list") { - N = sqrt(length(x)) - x = x[1:N] - y = matrix(y, ncol = N)[1, ] - density = list(x = x, y = y, z = matrix(density, ncol = N)) - } - - # Return Value: - density -} - - -# ------------------------------------------------------------------------------ - - -.gfunc2d = -function(x, param = NULL, type = c("norm", "cauchy", "t", "logistic", -"laplace", "kotz", "epower")) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Generator function for elliptical distributions - - # Note: - # A copy from fExtremes 'gfunc' - - # Arguments: - # x - a numeric vector - # param - NULL, a numeric value, or a numeric vector adding. - # additional parameters to the generator function. - # type - a character string denoting the type of distribution. - # This may be either - # "norm" for the normal distribution, or - # "cauchy" for the Cauchy distribution, or - # "t" for the Student-t distribution, or - # "logistic" for the logistic distribution, or - # "laplace" for the distribution, or - # "kotz" for the original Kotz distribution, or - # "epower" for the exponential power distribution - - # Value: - # Returns a numeric vector "g(x)" for the generator computed at - # the x values taken from the input vector. If x is missing, - # the normalizing constant "lambda" will be returned. - - # FUNCTION: - - # Handle Missing x: - if (missing(x)) { - x = NA - output = "lambda" - } else { - output = "g" - } - - # Get Type: - type = type[1] - - # Get Parameters: - # if (is.null(param)) param = .ellipticalParam$param - - # Create Generator: - if (type == "norm") { - g = exp(-x/2) - lambda = 1 / (2*pi) - param = NULL - } - if (type == "cauchy") { - g = ( 1 + x )^ (-3/2 ) - lambda = 1 / (2*pi) - param = NULL - } - if (type == "t") { - if (is.null(param)) { - nu = 4 - } else { - nu = param[[1]] - } - g = ( 1 + x/nu )^ ( -(nu+2)/2 ) - lambda = 1/(2*pi) - param = c(nu = nu) - } - if (type == "logistic"){ - g = exp(-x/2)/(1+exp(-x/2))^2 - # lambda: - # integrate(function(x) { exp(-x)/(1+exp(-x))^2}, 0, Inf, - # subdivision = 10000, rel.tol = .Machine$double.eps^0.8) - # 0.5 with absolute error < 2.0e-13 - lambda = 1 / pi - param = NULL - } - if (type == "laplace") { # or "double exponential" - # epower: - r = sqrt(2) - s = 1/2 - g = exp(-r*(x/2)^s) - lambda = s * r^(1/s) / ( 2 * pi * gamma(1/s) ) - param = NULL - } - if (type == "kotz") { - # epower: s = 1 - if (is.null(param)) { - r = sqrt(2) - } else { - r = param - } - g = exp(-r*(x/2)) - lambda = r/(2*pi) - param = c(r = r) - } - if (type == "epower") { - if (is.null(param)) { - r = sqrt(2) - s = 1/2 - } else { - r = param[[1]] - s = param[[2]] - } - g = exp(-r*(x/2)^s) - lambda = s * r^(1/s) / ( 2 * pi * gamma(1/s) ) - param = c(r = r, s = s) - } - - # Output: - output = output[1] - if (output == "g") { - ans = g - } else if (output == "lambda") { - ans = lambda - } - - # Add Control: - if (output == "g") { - attr(ans, "control") = c(type = type, lambda = as.character(lambda)) - } else if (output == "lambda") { - if (is.null(param)) { - attr(ans, "control") = unlist(list(type = type)) - } else { - attr(ans, "control") = unlist(list(type = type, param = param)) - } - } - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.delliptical2dSlider = -function(B = 10, eps = 1.e-3) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Displays interactively perspective plots of density - - #FUNCTION: - - # Graphic Frame: - par(mfrow = c(1, 1), cex = 0.7) - - # Internal Function: - refresh.code = function(...) - { - # Sliders: - Distribution = .sliderMenu(no = 1) - N = .sliderMenu(no = 2) - rho = .sliderMenu(no = 3) - nu = .sliderMenu(no = 4) - r = .sliderMenu(no = 5) - s = .sliderMenu(no = 6) - nlev = .sliderMenu(no = 7) - ncol = .sliderMenu(no = 8) - if (rho == +1) rho = rho - eps - if (rho == -1) rho = rho + eps - - # Title: - Names = c("- Normal", "- Cauchy", "- Student t", "- Logistic", - "- Laplace", "- Kotz", "- Exponential Power") - Title = paste("Elliptical Density No:", as.character(Distribution), - Names[Distribution], "\nrho = ", as.character(rho)) - if (Distribution == 3) Title = paste(Title, "nu =", as.character(nu)) - if (Distribution >= 6) Title = paste(Title, "r =", as.character(r)) - if (Distribution >= 7) Title = paste(Title, "s =", as.character(s)) - - # Plot: - xy= grid2d(x = seq(-5, 5, length = N)) - Type = c("norm", "cauchy", "t", "logistic", "laplace", "kotz", "epower") - param = NULL - if (Distribution == 3) param = nu - if (Distribution == 6) param = r - if (Distribution == 7) param = c(r, s) - D = delliptical2d(x = xy, rho = rho, param = param, - type = Type[Distribution], output = "list") - image(D, col = heat.colors(ncol), xlab = "x", ylab = "y" ) - contour(D, nlevels = nlev, add = TRUE) - title(main = Title) - - # Reset Frame: - par(mfrow = c(1, 1), cex = 0.7) - } - - # Open Slider Menu: - plot.names = c("Plot - levels", "... colors") - .sliderMenu(refresh.code, - names = c("Distribution", "N", "rho", "t: nu", "r", "s", plot.names), - minima = c( 1, 10, -1, 1, 0, 0, 10, 12), - maxima = c( 7, 100, +1, B, B, B, 100, 256), - resolutions = c( 1, 10, 0.1, 0.1, 0.1, 0.1, 10, 1), - starts = c( 1, 10, 0, 4, 1, 1, 10, 12)) -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/bv-dnorm.R fcopulae-4021.84/inst/obsolete/R/bv-dnorm.R --- fcopulae-3042.82.1/inst/obsolete/R/bv-dnorm.R 2013-02-19 16:15:18.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/bv-dnorm.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,364 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: BIVARIATE NORMAL DISTRIBUTION: -# pnorm2d Computes bivariate Normal probability function -# dnorm2d Computes bivariate Normal density function -# rnorm2d Generates bivariate normal random deviates -################################################################################ - - -pnorm2d = -function(x, y = x, rho = 0) -{ # pnorm2d: A copy from R package "sn" - - # Description: - # Computes bivariate Normal probability function - - # Arguments: - # x, y - two numeric values or vectors of the same length at - # which the probability will be computed. - - # Value: - # returns a numeric vector of probabilities of the same length - # as the input vectors - - # FUNCTION: - - # Probaility: - X = cbind(x, y) - ans = apply(X, 1, .pnorm2d, rho = rho) - attr(ans, "control") = c(rho = rho) - - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.pnorm2d = -function(X, rho = 0) -{ # pnorm2d: A copy from R package "sn" - - # Description: - # Bivariate Normal probability function - - # Arguments: - # x, y - two numeric values at which the probability will - # be computed. - - # Value: - # returns a numeric vector of probabilities of the same length - # as the input vectors - - # FUNCTION: - - # Probability: - x = X[1] - y = X[2] - if (x == 0 & y == 0) { - return(0.25 + asin(rho)/(2 * pi)) - } - p = 0.5 * (pnorm(x) + pnorm(y)) - if (x == 0) { - p = p - 0.25 * sign(y) - } else { - if (is.finite(x)) { - Y = (y - rho * x)/(x * sqrt(1 - rho^2)) - } else { - Y = -rho/sqrt(1-rho^2) - } - p = p - .TOwen(x, Y) - } - if (y == 0) { - p = p - 0.25 * sign(x) - } else { - if (is.finite(y)) { - X = (x - rho * y)/(y * sqrt(1 - rho^2)) - } else { - X = -rho/sqrt(1-rho^2) - } - p = p - .TOwen(y, X) - } - if (is.finite(x) & is.finite(y)) { - if ((x * y < 0) | ((x * y == 0) & (x + y) < 0)) { - p = p - 0.5 - } - } - - # Return Value: - return(p) -} - - -# ------------------------------------------------------------------------------ - - -.TInt = -function(h, a, jmax, cut.point) -{ # T.int: A copy from R package "sn" - - # Note: - # Required by .pnorm2d and .TOwen - - # FUNCTION: - - .fui = function(h, i) (h^(2 * i))/((2^i) * gamma(i + 1)) - seriesL = seriesH = NULL - i = 0:jmax - low = (h <= cut.point) - hL = h[low] - hH = h[!low] - L = length(hL) - if (L > 0) { - b = outer(hL, i, .fui) - cumb = apply(b, 1, cumsum) - b1 = exp(-0.5 * hL^2) * t(cumb) - matr = matrix(1, jmax + 1, L) - t(b1) - jk = rep(c(1, -1), jmax)[1:(jmax + 1)]/(2 * i + 1) - matr = t(matr * jk) %*% a^(2 * i + 1) - seriesL = (atan(a) - as.vector(matr))/(2 * pi) - } - if (length(hH) > 0) { - seriesH = atan(a) * exp(-0.5 * (hH^2) * a/atan(a)) * - (1 + 0.00868 * (hH^4) * a^4)/(2 * pi) - } - series = c(seriesL, seriesH) - id = c((1:length(h))[low], (1:length(h))[!low]) - series[id] = series - - # Return Value: - series -} - - -# ------------------------------------------------------------------------------ - - -.TOwen = -function (h, a, jmax = 50, cut.point = 6) -{ # T.Owen: A copy from R package "sn" - - # Note: - # Required by .pnorm2d - - # FUNCTION: - - if (!is.vector(a) | length(a) > 1) - stop("a must be a vector of length 1") - if (!is.vector(h)) - stop("h must be a vector") - aa = abs(a) - ah = abs(h) - if (aa == Inf) - return(0.5 * pnorm(-ah)) - if (aa == 0) - return(rep(0, length(h))) - na = is.na(h) - inf = (ah == Inf) - ah = replace(ah, (na | inf), 0) - if (aa <= 1) { - owen = .TInt(ah, aa, jmax, cut.point) - } else { - owen = 0.5 * pnorm(ah) + pnorm(aa * ah) * (0.5 - pnorm(ah)) - - .TInt(aa * ah, (1/aa), jmax, cut.point) - } - owen = replace(owen, na, NA) - owen = replace(owen, inf, 0) - ans = return(owen * sign(a)) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -dnorm2d = -function(x, y = x, rho = 0) -{ # A function implemented by Diethelm Wuertz - - # Arguments: - # x,y - two numeric vectors - # rho - the linear correlation, a numeric value between - # minus one and one. - - # FUNCTION: - - # Argument: - xoy = (x^2 - 2*rho*x*y + y^2)/ (2*(1 - rho^2)) - - # Density: - density = exp(-xoy) / ( 2*pi*sqrt(1-rho^2)) - attr(density, "control") = c(rho = rho) - - # Return Value: - density -} - - -# ------------------------------------------------------------------------------ - - -.dnorm2d = -function(x, y = x, rho = 0) -{ # A function implemented by Diethelm Wuertz - - # Arguments: - # x,y - two numeric vectors - # rho - the linear correlation, a numeric value between - # minus one and one. - - # Note: - # Partly copied from contributed R package 'mvtnorm' - # Author Friedrich Leisch - - # FUNCTION - - # Settings: - mean = c(0,0) - sigma = diag(2) - sigma[1,2] = sigma[2,1] = rho - log = FALSE - x = cbind(x, y) - - # From mvtnorm - Check: - if (is.vector(x)) { - x = matrix(x, ncol = length(x)) - } - if (missing(mean)) { - mean = rep(0, length = ncol(x)) - } - if (missing(sigma)) { - sigma = diag(ncol(x)) - } - if (ncol(x) != ncol(sigma)) { - stop("x and sigma have non-conforming size") - } - if (nrow(sigma) != ncol(sigma)) { - stop("sigma meanst be a square matrix") - } - if (length(mean) != nrow(sigma)) { - stop("mean and sigma have non-conforming size") - } - - # From mvtnorm - Density: - distval = mahalanobis(x, center = mean, cov = sigma) - logdet = sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values)) - logretval = -(ncol(x)*log(2*pi) + logdet + distval)/2 - if(log) return(logretval) - ans = exp(logretval) - attr(ans, "control") = c(rho = rho) - - # Return value: - ans -} - - -# ------------------------------------------------------------------------------ - - -rnorm2d = -function(n, rho = 0) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Generates bivariate normal random deviates - - # Arguments: - # n - number of random deviates to be generated - # rho - the linear correlation, a numeric value between - # minus one and one. - - # Note: - # Partly copied from contributed R package 'mvtnorm' - # Author Friedrich Leisch - - # FUNCTION - - # Settings: - mean = c(0,0) - sigma = diag(2) - sigma[1,2] = sigma[2,1] = rho - - # From mvtnorm - Random Numbers: - ev = eigen(sigma, symmetric = TRUE)$values - if (!all(ev >= -sqrt(.Machine$double.eps) * abs(ev[1]))) - warning("sigma is numerically not positive definite") - sigsvd = svd(sigma) - ans = t(sigsvd$v %*% (t(sigsvd$u) * sqrt(sigsvd$d))) - ans = matrix(rnorm(n * ncol(sigma)), nrow = n) %*% ans - ans = sweep(ans, 2, mean, "+") - attr(ans, "control") = c(rho = rho) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.rnorm2d = -function(n, rho = 0) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Alternative direct algorithm from Lindskog Master Thesis - - # Arguments: - # n - number of random deviates to be generated - # rho - the linear correlation, a numeric value between - # minus one and one. - - # FUNCTION: - - # Random Deviates - x = matrix(c(1, rho, rho,1), 2) - V = NULL - U = chol(x) - siz = dim(x)[1] - for(i in 1:n) { - Z = rnorm(siz) - res = t(U)%*%Z - V = cbind(V,res) - } - rmn = t(V) - - # Return Value: - rmn -} - - -################################################################################ diff -Nru fcopulae-3042.82.1/inst/obsolete/R/bv-dt.R fcopulae-4021.84/inst/obsolete/R/bv-dt.R --- fcopulae-3042.82.1/inst/obsolete/R/bv-dt.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/bv-dt.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,148 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: BIVARIATE STUDENT-T DISTRIBUTION: -# pt2d Computes bivariate Student-t probability function -# dt2d Computes bivariate Student-t density function -# rt2d Generates bivariate Student-t random deviates -################################################################################ - - -pt2d = -function(x, y = x, rho = 0, nu = 4) -{ # pnorm2d: A copy from R package "sn" - - # Description: - # Computes bivariate Student-t probability function - - # Arguments: - # x, y - two numeric values or vectors of the same length at - # which the probability will be computed. - - # Example: - # pt2d(rnorm(5), rnorm(5), 0.5, 5) - - # Value: - # returns a numeric vector of probabilities of the same length - # as the input vectors - - # FUNCTION: - - # Normal Limit: - if (nu == Inf) return(pnorm2d(x = x, y = y, rho = rho)) - - # Settings: - sigma = diag(2) - sigma[1, 2] = sigma[2, 1] = rho - X = cbind(x, y) - - # Probaility: - ans = pmvst(X, dim = 2, mu = c(0, 0), Omega = sigma, - alpha = c(0, 0), df = nu) - attr(ans, "control") = c(rho = rho, nu = nu) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -dt2d = -function(x, y = x, rho = 0, nu = 4) -{ # A function implemented by Diethelm Wuertz - - # Arguments: - # n - number of random deviates to be generated - # rho - the linear correlation, a numeric value between - # minus one and one. - - # Description: - # Computes bivariate Student-t density function - - # Example: - # dt2d(rnorm(5), rnorm(5), 0.5, 5) - - # Note: - # Partly copied from contributed R package 'sn' - - # FUNCTION: - - # Normal Limit: - if (nu == Inf) return(dnorm2d(x = x, y = y, rho = rho)) - - # Argument: - xoy = (x^2 - 2*rho*x*y + y^2)/ (2*(1 - rho^2)) - - # Density: - density = (1 + 2*xoy/nu)^(-(nu+2)/2) / (2*pi*sqrt(1-rho^2)) - attr(density, "control") = c(rho = rho, nu = nu) - - # Return value: - density -} - - -# ------------------------------------------------------------------------------ - - -rt2d = -function(n, rho = 0, nu = 4) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Generates bivariate Student-t random deviates - - # Arguments: - # n - number of random deviates to be generated - # rho - the linear correlation, a numeric value between - # minus one and one. - - # Note: - # Partly copied from contributed R package 'mvtnorm' - # Author Friedrich Leisch - - # FUNCTION: - - # Normal Limit: - if (nu == Inf) return(rnorm2d(n = n, rho = rho)) - - # Random Deviates: - ans = rnorm2d(n, rho)/sqrt(rchisq(n, nu)/nu) - attr(ans, "control") = c(rho = rho, nu = nu) - - # Return Value: - ans -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/mv-distributions.R fcopulae-4021.84/inst/obsolete/R/mv-distributions.R --- fcopulae-3042.82.1/inst/obsolete/R/mv-distributions.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/mv-distributions.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,671 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: PARAMETER ESTIMATION: -# fMV S4 Object of class 'fMV' -# mvFit Fits a MV Normal or Student-t Distribution -# print.fMV S3: Print method for objects of class 'fMV' -# plot.fMV S3: Plot method for objects of class 'fMV' -# summary.fMV S3: Summary method for objects of class 'fMV' -# .mvnormFit Fits a Multivariate Normal Distribution -# .mvstFit Fits a Multivariate Student-t Distribution -# .mvsnormPlot Plots for Multivariate Normal Distribution -# .mvstPlot Plots for Multivariate Student-t Distribution -# REQUIREMENTS: DESCRIPTION: -# "mvtnorm" Contributed R - Package -# "sn" | "mnormt" Contributed R - Package -################################################################################ - - -################################################################################ -# PARAMETER FIT: - - -setClass("fMV", - representation( - call = "call", - method = "character", - model = "list", - data = "data.frame", - fit = "list", - title = "character", - description = "character") -) - - -# ------------------------------------------------------------------------------ - - -mvFit = -function(x, method = c("snorm", "st"), fixed.df = NA, title = NULL, -description = NULL, trace = FALSE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - - # FUNCTION: - - # Fit: - if (method[1] == "snorm") { - # Normal Fit: - fit = .mvsnormFit(x = x, trace = trace, ...) - fit$df = Inf - } - if (method[1] == "st") { - # Student-t Fit: - fit = .mvstFit(x = x, fixed.df = fixed.df, trace = trace, ...) - } - - # Add to fit: - fit$method = method[1] - class(fit) = "list" - - # Model Slot: - model = list(beta = fit$beta, Omega = fit$Omega, - alpha = fit$alpha, df = fit$df) - - # Title Slot: - if (is.null(title)) { - if (method[1] == "snorm") - title = "Multivariate Normal Distribution" - if (method[1] == "st") - title = "Multivariate Student-t Distribution" - } - - # Description Slot: - if (is.null(description)) description = description() - - # Return Value: - new("fMV", - call = as.call(match.call()), - method = as.character(method[1]), - model = model, - data = as.data.frame(x), - fit = fit, - title = as.character(title), - description = as.character(description) ) -} - - -# ------------------------------------------------------------------------------ - - - -setMethod("show", "fMV", - function(object) -{ # A function implemented by Diethelm Wuertz - - # Description: - - # Arguments: - - # FUNCTION: - - # Extract fit: - fit = object@fit - - # Print: - cat("\nCall:\n ") - print.default(fit$call) - - cat("\nParameter Sstimates:\n") - print.default(fit$dp) - - cat("\nParameter Errors:\n") - print.default(fit$se) - - # cat("\nOptimization:\n") - # print.default(fit$optim) -}) - - -# ------------------------------------------------------------------------------ - - -plot.fMV = -function(x, which = "ask", ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - - # Arguments: - - # FUNCTION: - - # Plot: - if (x@fit$method == "snorm") { - # Multivariate Skew Normal Distribution: - return(.mvsnormPlot(x = x@fit, which = which, ...)) - } - if (x@fit$method == "st") { - # Multivariate Skew Student-t Distribution: - return(.mvstPlot(x = x@fit, which = which, ...)) - } -} - - -# ------------------------------------------------------------------------------ - - -summary.fMV = -function(object, which = "ask", doplot = TRUE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - - # Arguments: - - # FUNCTION: - - # Print: - print(x = object, ...) - - # Plot: - if (doplot) plot(x = object, which = which, doplot, ...) - - # Return Value: - invisible(object) -} - - -################################################################################ -# INERNAL FUNCTIONS: - - -.mvsnormFit = -function(x, trace = FALSE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Internal Function - - # Arguments: - - # FUNCTION: - - # Settings: - y = x - y.name = deparse(substitute(y)) - y.names = dimnames(y)[[2]] - y = as.matrix(y) - colnames(y) = y.names - k = ncol(y) - freq = rep(1, nrow(y)) - n = sum(freq) - X = rep(1, nrow(y)) - X = as.matrix(X) - m = ncol(X) - dimnames(y) = list(NULL, outer("V", as.character(1:k), paste, sep = "")) - y.names = as.vector(dimnames(y)[[2]]) - qrX = qr(X) - - # Fit: - mle = msn.mle(X = X, y = y, freq = freq, trace = trace, ...) - mle$call = match.call() - mle$y = y - mle$y.names = y.names - - # Parameters: - mle$beta = beta = mle$dp$beta - mle$xi = xi = X %*% beta - mle$Omega = Omega = mle$dp$Omega - mle$alpha = alpha = mle$dp$alpha - - # Test: - # dev.norm = msn.dev(c(qr.coef(qrX, y), rep(0, k)), X, y, freq) - # test = dev.norm + 2 * mle$logL - # p.value = 1 - pchisq(test, k) - # mle$test.normality = list(LRT = test, p.value = p.value) - - # Save for Plot: - Xb = qr.fitted(qrX, y) - res = qr.resid(qrX, y) - mle$k = k - mle$n = n - mle$pp = qchisq((1:n)/(n + 1), k) - mle$rad.n = apply((y - Xb) * ((y - Xb) %*% solve(var(res))), 1, sum) - mle$rad.sn = apply((y - xi) * ((y - xi) %*% solve(Omega)), 1, sum) - - # Return Value: - class(mle) = "snFit" - mle -} - - -# ------------------------------------------------------------------------------ - - -.mvstFit = -function(x, fixed.df = NA, trace = FALSE, ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Internal Function - - # Arguments: - - # FUNCTION: - - # Settings: - y = as.matrix(x) - k = ncol(y) - y.name = deparse(substitute(y)) - dimnames(y) = list(NULL, outer("V", as.character(1:k), paste, sep = "")) - y.names = dimnames(y)[[2]] - - freq = rep(1, nrow(y)) - n = sum(freq) - - X = as.matrix(rep(1, nrow(y))) - qrX = qr(X) - m = ncol(X) - - # Fit: - mle = mst.mle(X = X, y = y, freq = freq, fixed.df = fixed.df, - trace = trace, ...) - mle$call = match.call() - mle$y = y - mle$y.names = y.names - - # Parameters: - mle$beta = beta = mle$dp$beta - mle$xi = xi = X %*% beta - mle$Omega = Omega = mle$dp$Omega - mle$alpha = alpha = mle$dp$alpha - mle$df = df = mle$dp$df - - # Save for Plot: - Xb = qr.fitted(qrX, y) - res = qr.resid(qrX, y) - mle$k = k - mle$n = n - mle$pp = k * qf((1:n)/(n + 1), k, df) - mle$rad.n = as.vector(apply(res * (res %*% solve(var(res))), 1, sum)) - mle$rad.sn = as.vector(apply((y - xi)*((y - xi) %*% solve(Omega)), 1, sum)) - - # Return Value: - class(mle) = "stFit" - mle -} - - -# ------------------------------------------------------------------------------ - - -.mvsnormPlot = -function(x, which = "ask", ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Internal Plot Function - - # Arguments: - # x - the slot @fit from an object of class "fMV" - - # FUNCTION: - - # Settings: - dim = ncol(x$y) - - # Plot Title: - plot1Title = "Scatterplots" - if (dim == 1) plot1Title = "Histogram Plot" - - # Plot: - interactivePlot( - x = x, - choices = c( - plot1Title, - "Normal QQ-Plot", - "Skew-Normal QQ-Plot", - "Normal PP-Plot", - "Skew-Normal PP-Plot"), - plotFUN = c( - ".mvsnorm.plot.1", - ".mvsnorm.plot.2", - ".mvsnorm.plot.3", - ".mvsnorm.plot.4", - ".mvsnorm.plot.5"), - which = which) - - # Return Value: - invisible(x) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.1 <- -function(x) -{ - # Plot: - dim = x$k - if(dim == 1) .mvsnorm.plot.1A(x) else .mvsnorm.plot.1B(x) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.1A <- -function(x) -{ - # Plot: - z = x - y0 <- z$y - xi0 <- apply(z$xi, 2, mean) - y0 <- as.vector(y0) - x <- seq(min(pretty(y0, 10)), max(pretty(y0, 10)), length = 100) - omega <- sqrt(diag(z$Omega)) - dp0 <- c(xi0, omega, z$alpha) - xlab <- z$y.name - hist(y0, prob = TRUE, breaks = "FD", xlab = xlab, - ylab = "density", border = "white", col = "steelblue4", - main = z$y.name) - lines(x, dsn(x, dp0[1], dp0[2], dp0[3])) - if (length(y0) < 201) - points(y0, rep(0, z$n), pch = 1) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.1B <- -function(x) -{ - # Plot: - opt = options() - options(warn = -1) - pairs( - x$y, - labels = x$y.names, - panel = function(x, y, Y, y.names, xi, Omega, alpha) { - for (i in 1:length(alpha)) { - if (all(Y[, i] == x)) - Ix = i - if (all(Y[, i] == y)) - Iy = i } - points(x, y) - marg = msn.marginal(xi, Omega, alpha, c(Ix, Iy)) - xi.marg = marg$xi - Omega.marg = marg$Omega - alpha.marg = marg$alpha - x1 = seq(min(x), max(x), length = 30) - x2 = seq(min(y), max(y), length = 30) - dsn2.plot(x1, x2, xi.marg, Omega.marg, alpha.marg, - add = TRUE, col = "steelblue4")}, - Y = x$y, - y.names = dimnames(x$y)[[2]], - xi = apply(x$xi, 2, mean), - Omega = x$Omega, - alpha = x$alpha) - options(opt) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.2 <- -function(x) -{ - # Plot: - plot(x$pp, sort(x$rad.n), pch = 1, ylim = c(0, max(x$rad.n, x$rad.sn)), - xlab = "Chi-square Percentiles", - ylab = "Mahalanobis Distances") - abline(0, 1, lty = 3) - title(main = "Normal QQ-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.3 <- -function(x) -{ - # Plot: - plot(x$pp, sort(x$rad.sn), pch = 1, ylim = c(0, max(x$rad.n, x$rad.sn)), - xlab = "Percentiles of chi-square distribution", - ylab = "Mahalanobis distances") - abline(0, 1, lty = 3) - title(main = "Skew-Normal QQ-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.4 <- -function(x) -{ - # Plot: - plot((1:x$n)/(x$n + 1), sort(pchisq(x$rad.n, x$k)), - xlab = "", ylab = "") - abline(0, 1, lty = 3) - title(main = "Normal PP-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvsnorm.plot.5 <- -function(x) -{ - # Plot: - plot((1:x$n)/(x$n + 1), sort(pchisq(x$rad.sn, x$k)), - xlab = "", ylab = "") - abline(0, 1, lty = 3) - title(main = "Skew-Normal PP-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvstPlot = -function(x, which = "ask", ...) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Internal Plot Function - - # Arguments: - # x - the slot @fit from an object of class "fMV" - - # FUNCTION: - - # Settings: - dim = ncol(x$y) - - # Plot Title: - plot1Title = "Scatterplots" - if (dim == 1) plot1Title = "Histogram Plot" - - # Plot: - plot1Title = "Scatterplots" - if (dim == 1) plot1Title = "Histogram Plot" - interactivePlot( - x = x, - choices = c( - plot1Title, - "Normal QQ-Plot", - "Skew-Normal QQ-Plot", - "Normal PP-Plot", - "Skew-Normal PP-Plot"), - plotFUN = c( - ".mvst.plot.1", - ".mvst.plot.2", - ".mvst.plot.3", - ".mvst.plot.4", - ".mvst.plot.5"), - which = which) - - # Return Value: - invisible(x) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.1 <- -function(x) -{ - # Plot: - dim = x$k - if(dim == 1) .mvst.plot.1A(x) else .mvst.plot.1B(x) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.1A <- -function(x) -{ - # Plot: - z = x - y0 <- z$y - xi0 <- apply(z$xi, 2, mean) - y0 <- as.vector(y0) - x <- seq(min(pretty(y0, 10)), max(pretty(y0, 10)), length = 100) - omega <- sqrt(diag(z$Omega)) - dp0 <- c(xi0, omega, z$alpha, z$df) - xlab <- z$y.name - hist(y0, prob = TRUE, breaks = "FD", xlab = xlab, - ylab = "density", border = "white", col = "steelblue4", - main = z$y.name) - lines(x, dst(x, dp0[1], dp0[2], dp0[3], dp0[4])) - if (length(y0) < 201) - points(y0, rep(0, z$n), pch = 1) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.1B <- -function(x) -{ - # Plot: - opt = options() - options(warn = -1) - pairs( - x$y, - labels = x$y.names, - panel = function(x, y, Y, y.names, xi, Omega, alpha, df) { - for (i in 1:length(alpha)) { - if (all(Y[, i] == x)) - Ix = i - if (all(Y[, i] == y)) - Iy = i } - points(x, y) - marg = msn.marginal(xi, Omega, alpha, c(Ix, Iy)) - xi.marg = marg$xi - Omega.marg = marg$Omega - alpha.marg = marg$alpha - x1 = seq(min(x), max(x), length = 30) - x2 = seq(min(y), max(y), length = 30) - dst2.plot(x1, x2, xi.marg, Omega.marg, alpha.marg, - df, add = TRUE, col = "steelblue4")} , - Y = x$y, - y.names = dimnames(x$y)[[2]], - xi = apply(x$xi, 2, mean), - Omega = x$Omega, - alpha = x$alpha, - df = x$df) - options(opt) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.2 <- -function(x) -{ - # Plot: - plot(x$pp, sort(x$rad.n), pch = 1, ylim = c(0, max(x$rad.n, x$rad.sn)), - xlab = "Chi-square Percentiles", - ylab = "Mahalanobis Distances") - abline(0, 1, lty = 3) - title(main = "Normal QQ-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.3 <- -function(x) -{ - # Plot: - plot(x$pp, sort(x$rad.sn), pch = 1, ylim = c(0, max(x$rad.n, x$rad.sn)), - xlab = "Percentiles of chi-square distribution", - ylab = "Mahalanobis distances") - abline(0, 1, lty = 3) - title(main = "Skew-Normal QQ-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.4 <- -function(x) -{ - # Plot: - plot((1:x$n)/(x$n + 1), sort(pchisq(x$rad.n, x$k)), - xlab = "", ylab = "") - abline(0, 1, lty = 3) - title(main = "Normal PP-Plot", sub = x$y.name) -} - - -# ------------------------------------------------------------------------------ - - -.mvst.plot.5 <- -function(x) -{ - # Plot: - plot((1:x$n)/(x$n + 1), sort(pchisq(x$rad.sn, x$k)), - xlab = "", ylab = "") - abline(0, 1, lty = 3) - title(main = "Skew-Normal PP-Plot", sub = x$y.name) -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/mv-dsnorm.R fcopulae-4021.84/inst/obsolete/R/mv-dsnorm.R --- fcopulae-3042.82.1/inst/obsolete/R/mv-dsnorm.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/mv-dsnorm.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,181 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: DESCRIPTION: -# dmvsnorm Multivariate Skew Normal Density Function -# pmvsnorm Multivariate Skew Normal Probability Function -# rmvsnorm Multivariate Skew Normal Random Deviates -# REQUIREMENTS: DESCRIPTION: -# "mvtnorm" Contributed R - Package -# "sn" | "mnormt" Contributed R - Package -################################################################################ - - -################################################################################ -# Multivariate Skew Normal Distribution - - -dmvsnorm = -function(x, dim = 2, mu = rep(0, dim), Omega = diag(dim), -alpha = rep(0, dim)) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Multivariate Skew Normal Density Function - - # Note: - # Requires dsn() and dmsn() from R package sn - - # FUNCTION: - - # Settings: - xi = mu - ans = NA - - # Univariate Case: - if (is.vector(x) & dim == 1) { - ans = dsn(x, location = xi[1], scale = as.vector(Omega)[1], - shape = alpha[1]) - } - - # Multivariate Case: - if (is.matrix(x)) { - if (dim == ncol(x)) { - ans = dmsn(x = x, xi = xi, Omega = Omega, alpha = alpha) - } - } - - # Check for conflicting Dimensions: - if (is.na(ans[1])) { - stop("conflicting x and dim") - } - - # Return Value: - as.vector(ans) -} - - -# ------------------------------------------------------------------------------ - - -pmvsnorm = -function(q, dim = 2, mu = rep(0, dim), Omega = diag(dim), -alpha = rep(0, dim)) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Multivariate Skew Normal Probability Function - - # Algorithm: - - # Note: - # Requires psn() and pmsn() from R package sn - - # FUNCTION: - - # Settings: - x = q - xi = mu - ans = NA - - # Univariate Case: - if (is.vector(x) & dim == 1) { - ans = psn(x, location = xi[1], scale = as.vector(Omega)[1], - shape = alpha[1]) - } - - # Multivariate Case: - if (is.matrix(x)) { - if (dim == ncol(x)) { - ans = NULL - for (i in 1:nrow(x) ) { - ans = c(ans, pmsn(x = x[i,], xi = xi, Omega = Omega, - alpha = alpha)) - } - } - } - - # Check for conflicting Dimensions: - if (is.na(ans[1])) { - stop("conflicting x and dim") - } - - # Return Value: - as.vector(ans) -} - - -# ------------------------------------------------------------------------------ - - -rmvsnorm = -function(n, dim = 2, mu = rep(0, dim), Omega = diag(dim), -alpha = rep(0, dim)) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Multivariate Skew Normal Random Number Generator - - # Algorithm: - - # Note: - # Requires rsn() and rmsn() from R package sn - - # FUNCTION: - - # Settings: - ans = NA - xi = mu - - # Univariate Case: - if (dim == 1) { - ans = as.matrix(rsn(n, location = xi[1], - scale = as.vector(Omega)[1], shape = alpha[1])) - } - - # Multivariate Case: - if (dim > 1) { - ans = rmsn(n, xi = xi, Omega = Omega, alpha = alpha) - } - - # Check for conflicting Dimensions: - if (is.na(ans[1])) { - stop("dim must be greater 1") - } - - # Return Value: - rownames(ans) = as.character(1:n) - colnames(ans) = as.character(1:dim) - ans -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/R/mv-dst.R fcopulae-4021.84/inst/obsolete/R/mv-dst.R --- fcopulae-3042.82.1/inst/obsolete/R/mv-dst.R 2010-05-05 07:44:38.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/R/mv-dst.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,173 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2007, Diethelm Wuertz, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ -# FUNCTION: DESCRIPTION: -# dmvst Multivariate Skew Sudent-t Density Function -# pmvst Multivariate Skew Sudent-t Probability Function -# rmvst Multivariate Skew Sudent-t Random Deviates -# REQUIREMENTS: DESCRIPTION: -# "mvtnorm" Contributed R - Package -# "sn" | "mnormt" Contributed R - Package -################################################################################ - - -################################################################################ - - -dmvst = -function(x, dim = 2, mu = rep(0, dim), Omega = diag(dim), -alpha = rep(0, dim), df = 4) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Multivariate Skew Sudent-t Density Function - - # Arguments: - - # FUNCTION: - - # Settings: - xi = mu - ans = NA - - # Univariate Case: - if (is.vector(x) & dim == 1) { - ans = dst(x, location = xi[1], scale = as.vector(Omega)[1], - shape = alpha[1], df = Inf) - } - - # Multivariate Case: - if (is.matrix(x)) { - if (dim == ncol(x)) { - ans = dmst(x = x, xi = xi, Omega = Omega, alpha = alpha, df = df) - } - } - - # Check for conflicting Dimensions: - if (is.na(ans[1])) { - stop("conflicting x and dim") - } - - # Return Value: - as.vector(ans) -} - - -# ------------------------------------------------------------------------------ - - -pmvst = -function(q, dim = 2, mu = rep(0, dim), Omega = diag(dim), -alpha = rep(0, dim), df = 4) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Multivariate Skew Sudent-t Probability Function - - # Arguments: - - # FUNCTION: - - # Settings: - x = q - xi = mu - ans = NA - - # Univariate Case: - if (is.vector(x) & dim == 1) { - ans = pst(x, location = xi[1], scale = as.vector(Omega)[1], - shape = alpha[1], df = df) - } - - # Multivariate Case: - if (is.matrix(x)) { - if (dim == ncol(x)) { - ans = NULL - for (i in 1:nrow(x) ) { - ans = c(ans, pmst(x = x[i,], xi = xi, Omega = Omega, - alpha = alpha, df = df)) - } - } - } - - # Check for conflicting Dimensions: - if (is.na(ans[1])) { - stop("conflicting x and dim") - } - - # Return Value: - as.vector(ans) -} - - -# ------------------------------------------------------------------------------ - - -rmvst = -function(n, dim = 2, mu = rep(0, dim), Omega = diag(dim), -alpha = rep(0, dim), df = 4) -{ # A function implemented by Diethelm Wuertz - - # Description: - # Multivariate Skew Sudent-t Random Number Generator - - # Arguments: - - # FUNCTION: - - # Settings: - ans = NA - xi = mu - - # Univariate Case: - if (dim == 1) { - ans = as.matrix(rst(n, location = xi[1], - scale = as.vector(Omega)[1], shape = alpha[1], df = df)) - } - - # Multivariate Case: - if (dim > 1) { - ans = rmst(n, xi = xi, Omega = Omega, alpha = alpha, df = df) - } - - # Check for conflicting Dimensions: - if (is.na(ans[1])) { - stop("dim must be greater 1") - } - - # Return Value: - rownames(ans) = as.character(1:n) - colnames(ans) = as.character(1:dim) - ans -} - - -################################################################################ - diff -Nru fcopulae-3042.82.1/inst/obsolete/src/adapt2.f fcopulae-4021.84/inst/obsolete/src/adapt2.f --- fcopulae-3042.82.1/inst/obsolete/src/adapt2.f 2013-03-18 05:17:50.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/src/adapt2.f 1970-01-01 00:00:00.000000000 +0000 @@ -1,622 +0,0 @@ -CDW The multivariate integration package adapt was added to for use in the -CDW Rmetrics package fCopula. Thanks to Prof. Alan Genz who put his code -CDW for the use in fCopulae under the GPL-2 License. -CDW Message-ID: <4AD7A74B.3020108@math.wsu.edu> -CDW Date: Thu, 15 Oct 2009 15:50:51 -0700 -CDW From: Alan Genz -CDW User-Agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.8.1.21) -CDW Gecko/20090402 SeaMonkey/1.1.16 -CDW MIME-Version: 1.0 -CDW To: Diethelm Wuertz -CDW CC: Alan C Genz -CDW Subject: Re: adapt -CDW References: <4AD3032B.4090801@itp.phys.ethz.ch> -CDW In-Reply-To: <4AD3032B.4090801@itp.phys.ethz.ch> -CDW Content-Type: text/plain; charset=ISO-8859-1; format=flowed -CDW Content-Transfer-Encoding: 7bit -CDW Status: O -CDW Dear Prof. Wuertz, -CDW Thank you for your message and your interest in my adaptive integration -CDW Fortran code. I would be pleased if you included my code in your open -CDW source R fCopulae package under the Gnu GPL2 license. You have my -CDW permission to do this. -CDW Sincerely, -CDW Alan Genz - - -cMM this is the original adapt code with one modification. -cMM instead of calling the external function "functn", a fixed -cMM external routine adphlp is always called, and passed a pointer -cMM to the external s function. -cMM Michael Meyer, October 1989. - - subroutine adapt(ndim,a,b,minpts,maxpts,eps,relerr, - * lenwrk,wrkstr,finest,ifail) -c***begin prologue adapt -c adaptive multidimensional integration subroutine -c author: A. C. Genz, Washington State University -c 19 March 1984 -c************** parameters for adapt ******************************** -c***** input parameters -c ndim number of variables, must exceed 1, but not exceed 20 -c a real array of lower limits, with dimension ndim -c b real array of upper limits, with dimension ndim -c minpts minimum number of function evaluations to be allowed. -c on the first call to adapt minpts should be set to a -c non negative value (caution... minpts is altered by adapt). -c It is possible to continue a calculation to greater accuracy -c by calling adapt again by decreasing eps (described below) -c and resetting minpts to any negative value. -c minpts must not exceed maxpts. -c maxpts maximum number of function evaluations to be allowed, -c which must be at least rulcls, where -c rulcls = 2**ndim+2*ndim**2+6*ndim+1 -c -c for ndim = 2 3 4 5 6 7 8 9 10 12 15 20 -c maxpts >= rulcls = 25 45 73 113 173 269 433 729 1285 4457 33309 1049497 -c -c a SUGGESTED value for maxpts is 100 times the above values. -c -c functn externally declared user defined function to be integrated. -c it must have parameters (ndim,z), where z is a real array -c of dimension ndim. -cTSL this function has been replaced by the fixed function adhlp -c eps required relative accuracy -c lenwrk length of array wrkstr of working storage, the routine -c needs (2*ndim+3)*(1+maxpts/rulcls)/2 for lenwrk if -c maxpts function calls are used. -c for guidance, if you set maxpts to 100*rulcls (see table -c above) then acceptable values for lenwrk are -c -c for ndim = 2 3 4 5 6 7 8 9 -c lenwrk = 357 561 1785 3417 6681 13209 26265 52377 -c -c***** OUTPUT parameters -c -c minpts actual number of function evaluations used by adapt -c wrkstr real array of working storage of dimension (lenwrk). -c relerr estimated relative accuracy of finest -c finest estimated value of integral ["FINal ESTimate"] -c ifail : return code -c -c ifail=0 for normal exit, when estimated relative accuracy -c relerr is less than eps with maxpts or less function -c calls made. -c ifail=1 if maxpts was too small for adapt to obtain the -c required relative accuracy eps. -c In this case adapt returns a value of finest -c with estimated relative accuracy relerr. -c ifail=2 if lenwrk too small for maxpts function calls. -c In this case adapt returns a value of finest with -c estimated accuracy relerr using the working storage -c available, but relerr will be greater than eps. -c ifail=3 if ndim < 2, ndim > 20, -c ifail=4 if minpts > maxpts, -c ifail=5 if maxpts < rulcls or other memory problems -c (which will only be found later) -c*********************************************************************** -c***end prologue adapt - - implicit none - -C-- Arguments: -C double precision functn -C external functn - - integer ndim, minpts,maxpts, lenwrk, ifail - double precision a(ndim), b(ndim), eps, relerr, wrkstr(lenwrk), - & finest - -C-- Local Variables: - double precision center(20), width(20) - double precision errmin, rgnerr, rgnval, half, zero,one,two - - integer divaxo, divaxn, divflg, funcls, index1, index2, - * j, k, maxcls, rgnstr, rulcls, sbrgns, sbtmpp, subrgn, subtmp - - data zero/0d0/, one/1d0/, two/2d0/ - -c Check arguments; fail w/ code '3' or '4' - relerr=one - funcls=0 - ifail=3 - if(ndim.lt.2.or.ndim.gt.20) goto 990 - ifail=4 - if(minpts.gt.maxpts) goto 990 - ifail=5 -c -c***** initialisation of subroutine -c - half=one/two - rgnstr =2*ndim+3 - errmin = zero - maxcls = 2**ndim + 2*ndim**2 + 6*ndim+1 - maxcls = min0(maxcls,maxpts) - divaxo=0 -c -c***** end subroutine initialisation - if(minpts.lt.0) then - sbrgns=wrkstr(lenwrk-1) - goto 280 - endif - - do 30 j=1,ndim - width(j)=(b(j)-a(j))*half - 30 center(j)=a(j)+width(j) - finest=zero - wrkstr(lenwrk)=zero - divflg=1 - subrgn=rgnstr - sbrgns=rgnstr - -C-- REPEAT --- (outermost loop) ------- - 40 call bsrl(ndim,center,width,maxcls,rulcls, - * errmin,rgnerr,rgnval,divaxo,divaxn) - finest=finest+rgnval - wrkstr(lenwrk)=wrkstr(lenwrk)+rgnerr - funcls = funcls + rulcls -c -c***** place results of basic rule into partially ordered list -c***** according to subregion error - if(divflg .eq. 0) then -c -c***** when divflg=0 start at top of list and move down list tree to -c find correct position for results from first half of recently -c divided subregion - 200 subtmp=2*subrgn - if(subtmp.le.sbrgns) then - if(subtmp.ne.sbrgns) then - sbtmpp=subtmp+rgnstr - if(wrkstr(subtmp).lt.wrkstr(sbtmpp)) subtmp=sbtmpp - endif - 210 if(rgnerr.lt.wrkstr(subtmp)) then - do 220 k=1,rgnstr - index1=subrgn-k+1 - index2=subtmp-k+1 - wrkstr(index1)=wrkstr(index2) - 220 continue - subrgn=subtmp - - goto 200 - endif - endif - - else -c -c*****when divflg=1 start at bottom right branch and move up list -c tree to find correct position for results from second half of -c recently divided subregion - 230 subtmp=(subrgn/(rgnstr*2))*rgnstr - if(subtmp.ge.rgnstr) then - if(rgnerr.gt.wrkstr(subtmp)) then - do 240 k=1,rgnstr - index1=subrgn-k+1 - index2=subtmp-k+1 - wrkstr(index1)=wrkstr(index2) - 240 continue - subrgn=subtmp - goto 230 - endif - endif - endif - -c***** store results of basic rule in correct position in list - 250 wrkstr(subrgn)=rgnerr - wrkstr(subrgn-1)=rgnval - wrkstr(subrgn-2)=divaxn - do 260 j=1,ndim - subtmp=subrgn-2*(j+1) - wrkstr(subtmp+1)=center(j) - wrkstr(subtmp)=width(j) - 260 continue - if(divflg .eq. 0) then -c*** when divflg=0 prepare for second application of basic rule - center(divaxo)=center(divaxo)+two*width(divaxo) - sbrgns=sbrgns+rgnstr - subrgn=sbrgns - divflg=1 -c*** loop back to apply basic rule to other half of subregion - go to 40 - endif -c -c***** end ordering and storage of basic rule results -c***** make checks for possible termination of routine -c - 270 relerr=one - if(wrkstr(lenwrk).le.zero) wrkstr(lenwrk)=zero - if(dabs(finest).ne.zero) relerr=wrkstr(lenwrk)/dabs(finest) - if(relerr.gt.one) relerr=one - - if(sbrgns+rgnstr.gt.lenwrk-2) ifail=2 - if(funcls+funcls*rgnstr/sbrgns.gt.maxpts) ifail=1 - if(relerr.lt.eps.and.funcls.ge.minpts) ifail=0 - if(ifail.lt.3) goto 990 -c -c***** prepare to use basic rule on each half of subregion with largest -c error - 280 divflg=0 - subrgn=rgnstr - subtmp = 2*sbrgns/rgnstr - maxcls = maxpts/subtmp - errmin = dabs(finest)*eps/dfloat(subtmp) - wrkstr(lenwrk)=wrkstr(lenwrk)-wrkstr(subrgn) - finest=finest-wrkstr(subrgn-1) - divaxo=wrkstr(subrgn-2) - do 290 j=1,ndim - subtmp=subrgn-2*(j+1) - center(j)=wrkstr(subtmp+1) - 290 width(j)=wrkstr(subtmp) - width(divaxo)=width(divaxo)*half - center(divaxo)=center(divaxo)-width(divaxo) -c -c***** loop back to apply basic rule -c - goto 40 -c -c***** termination point -c - 990 minpts=funcls - wrkstr(lenwrk-1)=sbrgns - return - end - subroutine bsrl(s, center,hwidth, maxvls,funcls, - * errmin,errest,basest,divaxo,divaxn) - - implicit none -C-- Arguments: - integer s - double precision center(s), hwidth(s) - integer maxvls,funcls, divaxo,divaxn - double precision errmin, errest, basest -C - EXTERNAL adphlp - double precision adphlp -C-- Local Variables: - double precision intvls(20), z(20), fulsms(200), weghts(200) - - integer intcls, i, mindeg, maxdeg, maxord, minord - integer ifail - - double precision zero, one, two, ten, dif, errorm, - * sum0, sum1, sum2, difmax, x1, x2 - - data zero/0d0/, one/1d0/, two/2d0/, ten/10d0/ - - maxdeg = 12 - mindeg = 4 - minord = 0 - do 10 maxord = mindeg,maxdeg - call symrl(s, center, hwidth, minord, maxord, intvls, - * intcls, 200, weghts, fulsms, ifail) - if (ifail.eq.2) goto 20 - errest = dabs(intvls(maxord) -intvls(maxord-1)) - errorm = dabs(intvls(maxord-1)-intvls(maxord-2)) - if (errest.ne.zero) - & errest = errest* - & dmax1(one/ten,errest/dmax1(errest/two,errorm)) - if (errorm.le. 5.*errest) goto 20 - if (2*intcls.gt.maxvls) goto 20 - if (errest.lt.errmin) goto 20 - 10 continue - 20 difmax = -1 - x1 = one/two**2 - x2 = 3.*x1 - do 30 i = 1,s - z(i) = center(i) - 30 continue -cmmm - sum0 = adphlp(s,z) - do 40 i = 1,s - z(i) = center(i) - x1*hwidth(i) -cmmm - sum1 = adphlp(s,z) - z(i) = center(i) + x1*hwidth(i) - - sum1 = sum1 + adphlp(s,z) - z(i) = center(i) - x2*hwidth(i) - - sum2 = adphlp(s,z) - z(i) = center(i) + x2*hwidth(i) - - sum2 = sum2 + adphlp(s,z) - z(i) = center(i) - - dif = dabs((sum1-two*sum0) - (x1/x2)**2*(sum2-two*sum0)) - if (dif.ge.difmax) then - difmax = dif - divaxn = i - endif - 40 continue - if (sum0.eq.sum0+difmax/two) divaxn = mod(divaxo,s) + 1 - basest = intvls(minord) - funcls = intcls + 4*s - return - end - double precision function flsm(s,center,hwidth,x,m,mp,maxord, - * g,sumcls) -c -c*** function to compute fully symmetric basic rule sum -c - integer s, m(s), mp(s), maxord, sumcls, ixchng, lxchng, i, l, - * ihalf, mpi, mpl - double precision g(maxord), x(s), intwgt, zero, one,two, intsum, - * center(s), hwidth(s) - double precision adphlp - - zero = 0 - one = 1 - two = 2 - - intwgt = one - do 10 i=1,s - mp(i) = m(i) - if (m(i).ne.0) intwgt = intwgt/two - intwgt = intwgt*hwidth(i) - 10 continue - sumcls = 0 - flsm = zero -c -c******* compute centrally symmetric sum for permutation mp - 20 intsum = zero - do 30 i=1,s - mpi = mp(i) + 1 - x(i) = center(i) + g(mpi)*hwidth(i) - 30 continue - 40 sumcls = sumcls + 1 -cmmm - intsum = intsum + adphlp(s,x) - do 50 i=1,s - mpi = mp(i) + 1 - if(g(mpi).ne.zero) hwidth(i) = -hwidth(i) - x(i) = center(i) + g(mpi)*hwidth(i) - if (x(i).lt.center(i)) go to 40 - 50 continue -c******* end integration loop for mp -c - flsm = flsm + intwgt*intsum - if (s.eq.1) return -c -c******* find next distinct permutation of m and loop back -c to compute next centrally symmetric sum - do 80 i=2,s - if (mp(i-1).le.mp(i)) go to 80 - mpi = mp(i) - ixchng = i - 1 - if (i.eq.2) go to 70 - ihalf = ixchng/2 - do 60 l=1,ihalf - mpl = mp(l) - imnusl = i - l - mp(l) = mp(imnusl) - mp(imnusl) = mpl - if (mpl.le.mpi) ixchng = ixchng - 1 - if (mp(l).gt.mpi) lxchng = l - 60 continue - if (mp(ixchng).le.mpi) ixchng = lxchng - 70 mp(i) = mp(ixchng) - mp(ixchng) = mpi - go to 20 - 80 continue -c***** end loop for permutations of m and associated sums -c - return - end - subroutine nxprt(prtcnt, s, m) -c -c*** subroutine to compute the next s partition -c - implicit none - integer s, m(s), prtcnt - - integer i,k, msum - - if (prtcnt.gt.0) go to 20 - do 10 i=1,s - m(i) = 0 - 10 continue - prtcnt = 1 - return - 20 prtcnt = prtcnt + 1 - msum = m(1) - if (s.eq.1) go to 60 - do 50 i=2,s - msum = msum + m(i) - if (m(1).le.m(i)+1) go to 40 - m(1) = msum - (i-1)*(m(i)+1) - do 30 k=2,i - m(k) = m(i) + 1 - 30 continue - return - 40 m(i) = 0 - 50 continue - 60 m(1) = msum + 1 - return - end - subroutine symrl(s, center, hwidth, minord, maxord, intvls, - * intcls, numsms, weghts, fulsms, fail) -c multidimensional fully symmetric rule integration subroutine -c -c this subroutine computes a sequence of fully symmetric rule -c approximations to a fully symmetric multiple integral. -c written by a. genz, mathematical institute, university of kent, -c canterbury, kent ct2 7nf, england -c -c************** parameters for symrl ******************************** -c*****input parameters -c s integer number of variables, must exceed 0 but not exceed 20 -c f externally declared user defined real function integrand. -c it must have parameters (s,x), where x is a real array -c with dimension s. -c minord integer minimum order parameter. on entry minord specifies -c the current highest order approximation to the integral, -c available in the array intvls. for the first call of symrl -c minord should be set to 0. otherwise a previous call is -c assumed that computed intvls(1), ... , intvls(minord). -c on exit minord is set to maxord. -c maxord integer maximum order parameter, must be greater than minord -c and not exceed 20. the subroutine computes intvls(minord+1), -c intvls(minord+2),..., intvls(maxord). -c g real array of dimension(maxord) of generators. -c all generators must be distinct and nonnegative. -c numsms integer length of array fulsms, must be at least the sum of -c the number of distinct partitions of length at most s -c of the integers 0,1,...,maxord-1. an upper bound for numsms -c when s+maxord is less than 19 is 200 -c******output parameters -c intvls real array of dimension(maxord). upon successful exit -c intvls(1), intvls(2),..., intvls(maxord) are approximations -c to the integral. intvls(d+1) will be an approximation of -c polynomial degree 2d+1. -c intcls integer total number of f values needed for intvls(maxord) -c weghts real working storage array with dimension (numsms). on exit -c weghts(j) contains the weight for fulsms(j). -c fulsms real working storage array with dimension (numsms). on exit -c fulsms(j) contains the fully symmetric basic rule sum -c indexed by the jth s-partition of the integers -c 0,1,...,maxord-1. -c fail integer failure output parameter -c fail=0 for successful termination of the subroutine -c fail=1 when numsms is too small for the subroutine to -c continue. in this case weghts(1), weghts(2), ..., -c weghts(numsms), fulsms(1), fulsms(2), ..., -c fulsms(numsms) and intvls(1), intvls(2),..., -c intvls(j) are returned, where j is maximum value of -c maxord compatible with the given value of numsms. -c fail=2 when parameters s,minord, maxord or g are out of -c range -c*********************************************************************** -cmmm external f -ctsl real f -ctsl double precision f -c*** for double precision change real to double precision -c in the next statement - integer d, i, fail, k(20), intcls, prtcnt, l, m(20), maxord, - * minord, modofm, numsms, s, sumcls - double precision intvls(maxord), center(s), hwidth(s), gisqrd, - * glsqrd, - * intmpa, intmpb, intval, one, fulsms(numsms), weghts(numsms), - * two, momtol, momnkn, momprd(20,20), moment(20), zero, g(20) - double precision flsm, wht -c patterson generators - data g(1), g(2) /0.0000000000000000,0.7745966692414833/ - data g(3), g(4) /0.9604912687080202,0.4342437493468025/ - data g(5), g(6) /0.9938319632127549,0.8884592328722569/ - data g(7), g(8) /0.6211029467372263,0.2233866864289668/ - data g(9), g(10), g(11), g(12) /0.1, 0.2, 0.3, 0.4/ -c -c*** parameter checking and initialisation - fail = 2 - maxrdm = 20 - maxs = 20 - if (s.gt.maxs .or. s.lt.1) return - if (minord.lt.0 .or. minord.ge.maxord) return - if (maxord.gt.maxrdm) return - zero = 0 - one = 1 - two = 2 - momtol = one - 10 momtol = momtol/two - if (momtol+one.gt.one) go to 10 - hundrd = 100 - momtol = hundrd*two*momtol - d = minord - if (d.eq.0) intcls = 0 -c*** calculate moments and modified moments - do 20 l=1,maxord - floatl = l + l - 1 - moment(l) = two/floatl - 20 continue - if (maxord.ne.1) then - do 40 l=2,maxord - intmpa = moment(l-1) - glsqrd = g(l-1)**2 - do 30 i=l,maxord - intmpb = moment(i) - moment(i) = moment(i) - glsqrd*intmpa - intmpa = intmpb - 30 continue - if (moment(l)**2.lt.(momtol*moment(1))**2) moment(l) = zero - 40 continue - endif - do 70 l=1,maxord - if (g(l).lt.zero) return - momnkn = one - momprd(l,1) = moment(1) - if (maxord.eq.1) go to 70 - glsqrd = g(l)**2 - do 60 i=2,maxord - if (i.le.l) gisqrd = g(i-1)**2 - if (i.gt.l) gisqrd = g(i)**2 - if (glsqrd.eq.gisqrd) return - momnkn = momnkn/(glsqrd-gisqrd) - momprd(l,i) = momnkn*moment(i) - 60 continue - 70 continue - fail = 1 -c -c*** begin LOOP -c for each d find all distinct partitions m with mod(m))=d -c - 80 prtcnt = 0 - intval = zero - modofm = 0 - call nxprt(prtcnt, s, m) - 90 if (prtcnt.gt.numsms) return -c -c*** calculate weight for partition m and fully symmetric sums -c*** when necessary -c - if (d.eq.modofm) weghts(prtcnt) = zero - if (d.eq.modofm) fulsms(prtcnt) = zero - fulwgt = wht(s,moment,m,k,modofm,d,maxrdm,momprd) - sumcls = 0 - if (weghts(prtcnt).eq.zero .and. fulwgt.ne.zero) fulsms(prtcnt) = - * flsm(s, center, hwidth, moment, m, k, maxord, g, sumcls) - intcls = intcls + sumcls - intval = intval + fulwgt*fulsms(prtcnt) - weghts(prtcnt) = weghts(prtcnt) + fulwgt - call nxprt(prtcnt, s, m) - if (m(1).gt.modofm) modofm = modofm + 1 - if (modofm.le.d) go to 90 -c -c*** end loop for each d - if (d.gt.0) intval = intvls(d) + intval - intvls(d+1) = intval - d = d + 1 - if (d.lt.maxord) go to 80 -c -c*** set failure parameter and return - fail = 0 - minord = maxord - return - end - double precision function wht(s, intrps, m, k, modofm, d, - * maxrdm, momprd) -c*** subroutine to calculate weight for partition m -c - integer s, m(s), k(s), d, maxrdm, mi, ki, m1, k1, modofm - double precision intrps(s), zero, momprd(maxrdm,maxrdm) - - zero = 0 - do 10 i=1,s - intrps(i) = zero - k(i) = 0 - 10 continue - m1 = m(1) + 1 - k1 = d - modofm + m1 - 20 intrps(1) = momprd(m1,k1) - if (s.eq.1) go to 40 - do 30 i=2,s - mi = m(i) + 1 - ki = k(i) + mi - intrps(i) = intrps(i) + momprd(mi,ki)*intrps(i-1) - intrps(i-1) = zero - k1 = k1 - 1 - k(i) = k(i) + 1 - if (k1.ge.m1) go to 20 - k1 = k1 + k(i) - k(i) = 0 - 30 continue - 40 wht = intrps(s) - return - end diff -Nru fcopulae-3042.82.1/inst/obsolete/src/adapt_callback.c fcopulae-4021.84/inst/obsolete/src/adapt_callback.c --- fcopulae-3042.82.1/inst/obsolete/src/adapt_callback.c 2013-03-18 05:17:50.000000000 +0000 +++ fcopulae-4021.84/inst/obsolete/src/adapt_callback.c 1970-01-01 00:00:00.000000000 +0000 @@ -1,63 +0,0 @@ -#include "S.h" -#include "Rinternals.h" - -/* Added declaration of FORTRAN by Yohan Chalabi */ -void F77_NAME(adapt)(int*, /* ndim */ - double*, /* lower */ - double*, /* upper */ - int*, /* minpts */ - int*, /* maxpts */ - double*, /* eps */ - double*, /* relerr */ - int*, /* lenwrk */ - double*, /* wrkstr */ - double*, /* finest */ - int*); /* ifail */ - -static SEXP rho; -static SEXP f; - -/* All this routine does is call the approriate fortran - function. We need this so as to properly pass the S function */ -/* changed to doubles for R by Thomas Lumley */ -void cadapt(int *ndim, double *lower, double *upper, - int *minpts, int *maxpts, - void *functn, void *env, - double *eps, double *relerr, - int *lenwrk, double *finest, int *ifail) -{ - double *wrkstr; - - wrkstr = (double *) S_alloc(*lenwrk, sizeof(double)); - - /* store the R function and its environment */ - rho=env; - f=functn; - - F77_CALL(adapt)(ndim,lower,upper,minpts,maxpts,eps,relerr,lenwrk, - wrkstr,finest,ifail); -} - -/* This is the fixed routine called by adapt */ -/* changed to double for R, also rewritten to use eval() */ - -double F77_NAME(adphlp)(int *ndim, double *z) -{ - SEXP args,resultsxp,callsxp; - double result; - int i; - - PROTECT(args=allocVector(REALSXP,*ndim)); - for (i=0;i<*ndim;i++){ - REAL(args)[i]=z[i]; - } - - PROTECT(callsxp=lang2( f,args)); - PROTECT(resultsxp=eval(callsxp,rho)); - - result=REAL(resultsxp)[0]; - - UNPROTECT(3); - - return(result); -} diff -Nru fcopulae-3042.82.1/MD5 fcopulae-4021.84/MD5 --- fcopulae-3042.82.1/MD5 2020-03-07 11:06:14.000000000 +0000 +++ fcopulae-4021.84/MD5 2022-07-20 08:20:02.000000000 +0000 @@ -1,12 +1,12 @@ -069404d67832a534f33ba999e3b2d142 *ChangeLog -c68167d424066d9bf28d476901194b8a *DESCRIPTION +456eba40e6c30f71057b679c81f707d1 *ChangeLog +2121663662b4e49bafd66ca0adac1678 *DESCRIPTION 9d365d9b1fff0b1a443e7a7bcdade222 *NAMESPACE b224c94b7b4cadb927e4bb7f10fb39d3 *R/ArchimedeanCopulae.R 2ecff1b4ef03dca8bd390416e70d7728 *R/ArchimedeanDependency.R cb195a3002bdeaf47dcf8213e0086f58 *R/ArchimedeanGenerator.R df46e935e2e8fe9c0b8389a11398c745 *R/ArchimedeanModelling.R d1448727aeb5c2bfb70ea051f3d29b14 *R/ArchimedeanSlider.R -d2708cbe8df88efac4bb1ef8196d369f *R/EllipticalCopulae.R +6bab46620fac6baf5493156b3c0874a9 *R/EllipticalCopulae.R d648d61dfc570f419e48ea47a03779aa *R/EllipticalDependency.R 4fa6db2ede6ff3999d42506461aa6c11 *R/EllipticalGenerator.R ff1b598b5b38c986cb9ae589348bca3c *R/EllipticalModelling.R @@ -18,19 +18,6 @@ 24b6e829324e611e40b53f598f99020e *R/aaaCopulaeClass.R cc66b4093da5eaf07f5692f598fed2a3 *R/aaaCopulaeEnv.R 03c202578e38501cae6fb4e1c599c7e2 *R/zzz.R -28904ec67ad108b7d933bed2d9602b65 *inst/obsolete/R/biv-binning.R -33ebc9f0c61e57bf206752e2ab08a84d *inst/obsolete/R/biv-density.R -bf80eda2207641e1486de2649263b1cb *inst/obsolete/R/biv-gridding.R -fc4cd13b043e75b73ab522de6f8a54c3 *inst/obsolete/R/builtin-adapt.R -fc65d1675df34e46e66a4588990ee940 *inst/obsolete/R/bv-dcauchy.R -dc89d7904bfd40ff3c1ef5c0fa96b2c0 *inst/obsolete/R/bv-delliptical.R -c0ae6dd505d87ba282fc4415ae929df6 *inst/obsolete/R/bv-dnorm.R -0b956cf8e243e6b48402b9db69f034fe *inst/obsolete/R/bv-dt.R -6728cc8658694dfd46b671e3e4c8c27b *inst/obsolete/R/mv-distributions.R -4c3f51acecd49217a4fe9d51392f3bf4 *inst/obsolete/R/mv-dsnorm.R -e6fe483c6db38c6eba55b045fa52df39 *inst/obsolete/R/mv-dst.R -0ae32d48faf982f676ee2721b1ea4e02 *inst/obsolete/src/adapt2.f -bf0cc4d5a94a7545fdb9d6dcd73306af *inst/obsolete/src/adapt_callback.c 5073775e2caf6351d1b1e15479b1c73c *inst/unitTests/Makefile b3dfb205136831d39e558c232cdb025a *inst/unitTests/runTests.R 7b2fb68b07b1bcb3031af0bbaae9401d *inst/unitTests/runit.ArchimedeanCopulae.R diff -Nru fcopulae-3042.82.1/R/EllipticalCopulae.R fcopulae-4021.84/R/EllipticalCopulae.R --- fcopulae-3042.82.1/R/EllipticalCopulae.R 2014-09-16 13:06:06.000000000 +0000 +++ fcopulae-4021.84/R/EllipticalCopulae.R 2022-07-18 13:58:20.000000000 +0000 @@ -204,9 +204,11 @@ # FUNCTION: - # Use: X = .rnorm2d(n, rho) or alternatively: - X = fMultivar:::.rnorm2d(n = n, rho = rho) - + ## Use: X = .rnorm2d(n, rho) or alternatively: + ## X = fMultivar:::.rnorm2d(n = n, rho = rho) + ## The above isn;t exported - use the exported function + X = fMultivar::rnorm2d(n = n, rho = rho) + # Generate Z <- NULL for(i in (1:n)) Z <- rbind(Z, pnorm(X [i,]))