[Distr-commits] r661 - in branches/distr-2.3/pkg: distr/R distrEx distrEx/R distrEx/man distrEx/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jul 8 19:15:47 CEST 2010
Author: ruckdeschel
Date: 2010-07-08 19:15:46 +0200 (Thu, 08 Jul 2010)
New Revision: 661
Added:
branches/distr-2.3/pkg/distrEx/R/kMAD.R
branches/distr-2.3/pkg/distrEx/man/kMAD.Rd
branches/distr-2.3/pkg/distrEx/src/kMad.c
Modified:
branches/distr-2.3/pkg/distr/R/NegbinomDistribution.R
branches/distr-2.3/pkg/distrEx/NAMESPACE
branches/distr-2.3/pkg/distrEx/R/AllGeneric.R
branches/distr-2.3/pkg/distrEx/src/distrEx.dll
Log:
-Nataliya produced a C-version of kMad; integrated to distrEx now
-parameter size in NegbinomDistribution can now be real-valued
Modified: branches/distr-2.3/pkg/distr/R/NegbinomDistribution.R
===================================================================
--- branches/distr-2.3/pkg/distr/R/NegbinomDistribution.R 2010-06-15 18:51:07 UTC (rev 660)
+++ branches/distr-2.3/pkg/distr/R/NegbinomDistribution.R 2010-07-08 17:15:46 UTC (rev 661)
@@ -23,9 +23,9 @@
if(length(size(object)) != 1)
stop("size has to be a numeric of length 1")
if(size(object) < 0)
- stop("size has to be a not negative natural")
- if(!identical(floor(size(object)), size(object)))
- stop("size has to be a not negative natural")
+ stop("size has to be non-negative")
+# if(!identical(floor(size(object)), size(object)))
+# stop("size has to be a not negative natural")
else return(TRUE)
})
Modified: branches/distr-2.3/pkg/distrEx/NAMESPACE
===================================================================
--- branches/distr-2.3/pkg/distrEx/NAMESPACE 2010-06-15 18:51:07 UTC (rev 660)
+++ branches/distr-2.3/pkg/distrEx/NAMESPACE 2010-07-08 17:15:46 UTC (rev 661)
@@ -38,7 +38,7 @@
"+", "*",
"name", "name<-",
"E", "var", "IQR", "skewness", "kurtosis",
- "sd", "median", "mad",
+ "sd", "median", "mad", "kMAD",
"m1df", "m2df",
"liesInSupport")
export("EuclCondition")
Modified: branches/distr-2.3/pkg/distrEx/R/AllGeneric.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/AllGeneric.R 2010-06-15 18:51:07 UTC (rev 660)
+++ branches/distr-2.3/pkg/distrEx/R/AllGeneric.R 2010-07-08 17:15:46 UTC (rev 661)
@@ -75,6 +75,11 @@
setGeneric("mad", function(x, ...) standardGeneric("mad"))
}
+if(!isGeneric("kMAD")){
+ setGeneric("kMAD", function(x, k, ...) standardGeneric("kMAD"))
+}
+
+
if(!isGeneric("skewness")){
setGeneric("skewness", function(x, ...) standardGeneric("skewness"))
}
Added: branches/distr-2.3/pkg/distrEx/R/kMAD.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/kMAD.R (rev 0)
+++ branches/distr-2.3/pkg/distrEx/R/kMAD.R 2010-07-08 17:15:46 UTC (rev 661)
@@ -0,0 +1,33 @@
+############################################################
+## ###
+## Median of absolute deviations for skewed distribution ###
+## ###
+############################################################
+
+
+setMethod("kMAD", signature(x = "numeric", k = "numeric"),
+ function(x, k=1, na.rm = TRUE, eps = .Machine$double.eps, ...){
+ if(na.rm) x <- x[!is.na(x)]
+ if(! length(k)==1) stop ("k has to be a numeric of length 1")
+ if(k<=0) stop ("k has to be strictly positive")
+ eps1 <- min(diff(unique(sort(x))))
+ erg <- .C("kMad", as.double(x),
+ as.integer(length(x)),
+ as.integer(k),
+ d = double(1),
+ eps = as.double(min(eps1,eps)),
+ PACKAGE="distrEx")
+ return(erg$d)
+ })
+
+setMethod("kMAD", signature(x = "UnivariateDistribution", k = "numeric"),
+ function(x, k=1, up = NULL, ...){
+ if(! length(k)==1) stop ("k has to be a numeric of length 1")
+ if(k<=0) stop ("k has to be strictly positive")
+ m <- median(x)
+ if(is.null(up)) up <- min(3*IQR(x),q(x)(1)-m,m-q(x)(0))
+ fun <- function(t)
+ {p(x)(m+k*t)-p(x)(m-t)-.5}
+ return(uniroot(fun,lower=0,upper=up)$root)
+ })
+
Added: branches/distr-2.3/pkg/distrEx/man/kMAD.Rd
===================================================================
--- branches/distr-2.3/pkg/distrEx/man/kMAD.Rd (rev 0)
+++ branches/distr-2.3/pkg/distrEx/man/kMAD.Rd 2010-07-08 17:15:46 UTC (rev 661)
@@ -0,0 +1,43 @@
+\name{kMAD}
+\alias{kMAD}
+\alias{kMAD-methods}
+\alias{kMAD,UnivariateDistribution,numeric-method}
+\alias{kMAD,numeric,numeric-method}
+
+\title{Asymmetric Median of Absolute Deviations for Skewed Distributions}
+\description{
+ Function for the computation of asymmetric median absolute deviation (kMAD)
+ It coincides with ordinary median absolute deviation (MAD) for \eqn{k=1}{k=1}.
+}
+\usage{
+kMAD(x,k,...)
+\S4method{kMAD}{numeric,numeric}(x, k = 1, na.rm = TRUE,
+ eps = .Machine$double.eps, ... )
+\S4method{kMAD}{UnivariateDistribution,numeric}(x, k = 1, up = NULL, ... )
+}
+\arguments{
+ \item{x}{a numeric vector or a distribution. }
+ \item{k}{numeric; tunning parameter for asymmetrical MAD; has to be of length 1 and larger than 1.}
+ \item{na.rm}{logical; if \code{TRUE} then \code{NA} values are stripped from \code{x} before computation takes place. }
+ \item{eps}{numeric; accuracy up to which to state equality of two numeric values }
+ \item{up}{numeric; upper bound for search interval; important in distributions without left/right endpoint.}
+ \item{\dots}{additional arguments for other functions; not used so far;}
+ }
+\details{
+ For kMAD (asymmetrial MAD) is a root of the equation:
+ \deqn{\mathop{\rm kMAD}(F,k) = \inf\{t>0\;\mid \;F(m+kt)-F(m-t)\ge 1/2 \}}{kMAD(F,k) = inf{t>0|F(m+kt)-F(m-t)>=1/2}},
+ where \code{F} is the cumulative distribution function, \code{m} is the median of \code{F}.
+ }
+\references{
+ Ruckdeschel, P., Horbenko, N. (2010): Robustness Properties for Generalized Pareto Distributions. ITWM Report 182.
+}
+\author{Peter Ruckdeschel \email{Peter.Ruckdeschel at itwm.fraunhofer.de},
+Nataliya Horbenko \email{Nataliya.Horbenko at itwm.fraunhofer.de}}
+\seealso{\code{\link{mad}}}
+\examples{
+x <- rnorm(100)
+kMAD(x,k=10)
+kMAD(Norm(),k=10)
+}
+\concept{estimator}
+\keyword{scale estimator}
Modified: branches/distr-2.3/pkg/distrEx/src/distrEx.dll
===================================================================
(Binary files differ)
Added: branches/distr-2.3/pkg/distrEx/src/kMad.c
===================================================================
--- branches/distr-2.3/pkg/distrEx/src/kMad.c (rev 0)
+++ branches/distr-2.3/pkg/distrEx/src/kMad.c 2010-07-08 17:15:46 UTC (rev 661)
@@ -0,0 +1,110 @@
+#include <stdlib.h>
+#include <stdio.h>
+
+int compare_doubles(const void *a,const void *b)
+{
+ double *da = (double*)a;
+ double *db = (double*)b;
+ if (*da > *db)
+ return 1;
+ else if (*da < *db)
+ return -1;
+ else
+ return 0;
+}
+
+void kMad(double *x, int *lx, int *kp, double *d, double *eps)
+{
+double m;
+int n2, r1, r2, i1, i2,j, k = (*kp);
+
+// special cases
+ if ((lx[0]) < 1) d[0] = 0;
+ if ((lx[0]) == 1) d[0] = x[0];
+
+ if(lx[0] > 1){
+ qsort(x, *lx,sizeof(double), compare_doubles);
+
+ if ((lx[0]) % 2 == 0){
+ r1 = (lx[0])/2-1;
+ r2 = (lx[0])/2;
+ m = (x[r1]+x[r2])/2;
+ n2 = lx[0]/2;
+ }else{
+ r1 = ((lx[0])+1)/2-1;
+ r2 = r1;
+ m = x[r1];
+ n2 = (lx[0]+1)/2;
+ }
+ i1 = r1;
+ while( ( x[i1] > m-eps[0] ) && (i1 > 0) ) i1--;
+ i2 = r2;
+ while( ( x[i2] < m+eps[0] ) && ( i2 < lx[0]-2) ) i2++;
+
+
+ j = i2 - i1 + 1;
+ d[0] = (x[i2]-x[i1])/2;
+
+ while(j < n2)
+ {int i1l=0, i2l=0, jl=0, i1r=0, i2r=0, jr=0;
+ double xll=0, xrl=0, ll=0, xlr=0, xrr=0, lr=0;
+ // l = left check values, r = right check values
+ // check left
+ if(i1>0){
+ i1l = i1 - 1;
+ xll = x[i1l];
+ xrl = m + k * (m-xll);
+ ll = xrl - xll;
+ i2l = i2;
+ while((x[i2l] <= xrl +eps[0]) &&(i2l < lx[0]-2)) i2l ++;
+ jl = i2l - i1l + 1;
+ }
+ // check right
+ if(jl< n2){
+ i1 = i1l;
+ i2 = i2l;
+ j = jl;
+ d[0] = m-xll;
+ }else{
+ if(i2< lx[0]-2){
+ i2r = i2 + 1;
+ xrr = x[i2r];
+ xlr = m - (xrr-m) / k;
+ lr = xrr - xlr;
+ i1r = i1;
+ while((x[i1r] >= xlr -eps[0]) &&(i1r > 0)) i1r --;
+ jr = i2r - i1r + 1;
+ }
+ // decide which is shorter
+ if(jr > jl){
+ i1 = i1l;
+ i2 = i2l;
+ j = jl;
+ d[0] = m-xll;
+ }else{
+ i1 = i1r;
+ i2 = i2r;
+ j = jr;
+ d[0] = m-xlr;
+ }
+ }
+ }
+ }
+}
+
+// int main (){
+// int k = 10;
+// int n = 10;
+// int i;
+// double a[n];
+// double d;
+//
+// for (i=0;i<n;i++){
+// a[i] = i*8;
+// }
+//
+// kMad(a,&n,&k,&d);
+// printf("hello: %f \n", d);
+// return 0;
+// }
+//
More information about the Distr-commits
mailing list