added plot position calculation

build.bash

@@ -0,0 +1,32 @@
+#!/bin/bash
+
+R=R
+if [ ${#} -gt 0 ]; then
+    if [ "$1" = "-dev" ]; then
+        R=Rdev  
+    fi
+fi
+
+echo "######## Removing building information..."
+rm -rf output
+
+echo "######## Copying DESCRIPTION and NAMESPACE to pkg directory..."
+cp build/DESCRIPTION pkg
+cp build/NAMESPACE pkg
+
+echo "######## Generate documentation..."
+$R -q -f roxygen.R
+
+echo "######## Building package in output..."
+mkdir output
+cd output
+$R CMD build ../pkg
+echo "######## Testing package..."
+for x in *.tar.gz 
+do 
+    $R CMD check --as-cran $x
+done
+
+echo "**BUILT USING $R"
+$R --version
+

build/DESCRIPTION

@@ -0,0 +1,14 @@
+Package: extremevalues
+Maintainer: Mark van der Loo <mark.vanderloo@gmail.com>
+License: GPL-3
+Title: Univariate outlier detection using model distributions.
+LazyData: no
+Type: Package
+LazyLoad: yes
+Author: Mark van der Loo
+Description: hallo
+Version: 3.0-0
+Depends: R (>= 2.13.0)
+URL: https://github.com/markvanderloo/rspa
+Date: 2012-07-15
+Suggests: knitr, testthat

document.bash

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+R -f roxygen.R
+R CMD Rd2pdf --force --no-preview -o manual.pdf ./pkg
+
+

examples/plotpositions.R

@@ -0,0 +1,10 @@
+
+# simple plot positions
+plotpositions(c(1,3,5))
+
+# plot positions: duplicates get mean position
+plotpositions(c(1,3,5,3))
+
+# plot positions: duplicates get original position
+plotpositions(c(1,3,5,3),duplicates='asis')
+

pkg/NEWS

@@ -0,0 +1,33 @@
+version 3.0-0
+- BREAKING CHANGE: Complete rewrite of the package.
+- plot positions of duplicates now get mean plot position by default.
+
+version 2.2
+- Fixed minor bug in namespace file (thanks to Duncan Murdoch for pointing that out)
+- Changed dependancy on gWidgets and gWidgetstcltk to "suggest" for faster loading.
+- Changed documentation of returned R2-value (thanks to Charles Lobos for useful discussion and remarks)
+- Changed calculation of variance of residuals for method II (slightly underestimated previously, thanks to Javier Orche Galindo for noting this)
+- Fixed a bug causing extimated parameters in method II/pareto not to be passed to output.
+
+version 2.1
+- Killed a bug in normal distribution / method II
+- Updated references, some typo's in the man pages
+- Added version 0.1 of GUI interface
+- Small stuff in helpfiles
+
+version 2.0
+- getOutliers.r - Changed rho default, added left limit, added some input
+  checks
+- getLognormalLimit.r - added left limit
+- getExponentialLimit.r - added left limit
+- getWeibullLimit.r - added left limit
+- getNormalLimit.r - added left limit
+- getParetoLimit.r - added left limit
+- Added Method II 
+- Renaming main functions
+- Revised and extended plotting capabilities
+- Major overhaul of the help files
+- Fixed a bug in the exponential fit
+- Changed to Makkonen's equation for plot positions
+
+

pkg/R/.R

@@ -0,0 +1,6 @@
+
+
+source('plotpositions.R')
+
+
+

pkg/R/extremevalues.R

@@ -0,0 +1,6 @@
+#' A package for univariate outlier detection
+#' 
+#' @name extremevalues
+#' @docType package
+#' @useDynLib extremevalues
+{}

pkg/R/plotpositions.R

@@ -0,0 +1,44 @@
+
+all_finite <- function(x){
+   stopifnot(is.numeric(x))
+   storage.mode(x) <- "double"
+   .Call("all_finite_double",x)
+}
+
+#' Compute QQ-plotpositions
+#'
+#' @section Details:
+#' The QQ-plot indices are computed as as follows. For each \eqn{x_i} the rank \eqn{j_i} is determined.
+#' The plot index for \eqn{x_i} is then given by \eqn{i/(N+1)}. When \eqn{\boldsymbol{x}} has duplicate
+#' entries, plot indices for the duplicates are computed for each duplicate and replaced with the mean (by default).
+#' Passing \code{duplicates='asis'}, the original rank is used.
+#'
+#'
+#' @param x numeric vector. An error is thrown when \code{x} contains nonnumerical or nonfinite values.
+#' @param duplicates A character constant determining how duplicate entries are indexed (see details). 
+#' @return plot positions \eqn{0< p < 1} for the values in \code{x}.
+#' @export
+#' @example ../examples/plotpositions.R
+#' 
+plotpositions <- function(x,duplicates=c('mean','asis')){
+   stopifnot(
+      is.numeric(x),
+      all_finite(x)
+   )
+   duplicates <- match.arg(duplicates)
+   N <- length(x)
+   i <- order(x)
+   if (identical(duplicates,'asis')){
+      pp <- (1:N)/(N+1)
+      return(pp[i])
+   } else {
+      storage.mode(x) <- "double"
+      pp <- .Call("R_plotpositions",x[i])
+   }
+   pp[i]
+}
+
+
+
+
+

pkg/inst/testPlotpositions.R

@@ -0,0 +1,22 @@
+
+library(testthat)
+
+context('QQ-plot positions')
+test_that('plot postion results', {
+   expect_equal( plotpositions(1), 0.5 )
+   expect_equal( plotpositions(c(1,1), c(0.5,0.5)))
+   expect_equal(
+      round(plotpositions(c(1,2)),8), 
+      round(c(1/3,2/3),8) 
+   )
+   expect_equal( 
+      round(plotpositions(1,1,duplicates='asis'),8),
+      round(c(1/3,2/3),8)
+   )
+})
+
+test_that('plot position errors',{
+   expect_error( plotposition(NA) )
+   expect_error( plotposition(NaN) )
+})
+

pkg/src/R_all_finite.c

@@ -0,0 +1,25 @@
+
+#include <R.h>
+#include <Rdefines.h>
+
+SEXP all_finite_double(SEXP x){
+   PROTECT(x);
+   double *xx = REAL(x);
+
+   SEXP y;
+   PROTECT(y = allocVector(LGLSXP,1));
+
+   int i, b = 1;
+   for (i=0; i<length(x); i++) {
+      b = R_finite(xx[i]);
+      if (!b) break;
+   }
+
+
+   LOGICAL(y)[0] = b;
+   UNPROTECT(2);
+
+   return y; 
+}
+
+

pkg/src/R_plotpositions.c

@@ -0,0 +1,21 @@
+
+#include <R.h>
+#include <Rdefines.h>
+#include "plotpositions.h"
+
+SEXP R_plotpositions(SEXP x){
+   PROTECT(x);
+
+   SEXP I;
+   PROTECT(I = allocVector(REALSXP, length(x)));
+
+   plotpositions(REAL(x), length(x), REAL(I)); 
+
+   UNPROTECT(2);
+
+   return I; 
+}
+
+
+
+

pkg/src/plotpositions.c

@@ -0,0 +1,55 @@
+
+// Sum of sequence i, i+1,...,j.
+// - "double" arithmetic is used to avoid integer overflow.
+double seq_sum(int i, int j){
+   int k = ( (j - i) + 1 );
+   int m = k / 2;
+   double u, di = (double) i, dj = (double) j;
+   if ( k % 2 == 0 ){
+      u = m * (dj + di);
+   } else {
+      u = m * (dj + di ) + dj - m;
+   }
+   return u;
+}
+
+
+// QQ-plot positions, replacing double entries with the mean position.
+// - Positions determined  as described in Makkonen (2008).
+// - It is assumed that x is sortend in ascending order!
+void plotpositions(double *x, int n, double *I){
+   int i = 0, j, k = 1;
+   double m, pp;
+   int N = n + 1;
+   while( i < n ){
+      // find number (m)  and index (i,i+1,...,j) of duplicates
+      j = i+1;
+      while ( x[j] == x[i] ){ 
+         j++;
+      };
+      m = (double) j-i;
+      // compute mean QQ-plot position and assign
+      pp = seq_sum(i+1,j)/(m*N);
+      for ( k = i; k < j; k++ ){ 
+         I[k] = pp;
+      }
+      // shift index
+      i = j;
+   }
+
+}
+
+/* example: uncomment to run.
+#include <stdio.h>
+void main(){
+   double x[4] = {1.0,3.0,3.0,4.0};
+   double I[4];
+
+   plotpositions(x, 4, I);
+   for (int i=0; i<4; i++){
+      printf("%4.2f, %4.2f", x[i], I[i]);
+      printf("\n");
+   }
+
+}
+*/

pkg/src/plotpositions.h

@@ -0,0 +1,10 @@
+
+#ifndef ev_plotpositions_h
+#define ev_plotpositions_h
+
+void plotpositions(double *, int, double *);
+
+#endif
+
+
+

pkg/tests/testpackage.R

@@ -0,0 +1,5 @@
+
+if ( require(testthat) )
+   test_package("rspa")
+
+

roxygen.R

@@ -0,0 +1,14 @@
+library(roxygen2)
+options(error=traceback)
+unlink( 'pkg/man', TRUE)
+
+setwd('pkg')
+roxygenize( '.'
+          , roxygen.dir='.'
+          , copy.package=FALSE
+          , unlink.target=TRUE
+		    )
+
+if (length(list.files('inst/doc')) == 0){
+   unlink( 'inst/doc', TRUE)   
+}

test/autotest.R

@@ -0,0 +1,7 @@
+library(testthat)
+
+
+
+dyn.load("../pkg/src/plotpositions.so")
+auto_test("../pkg/R", "../pkg/inst/tests")
+