[Robast-commits] r513 - branches/robast-0.9/pkg/RobExtremes/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 12 17:27:01 CEST 2012
Author: horbenko
Date: 2012-09-12 17:27:01 +0200 (Wed, 12 Sep 2012)
New Revision: 513
Added:
branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R
branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R
branches/robast-0.9/pkg/RobExtremes/R/plotScaledIC.R
Log:
plotting Functions for RobExtremes and test functions
Added: branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R 2012-09-12 15:27:01 UTC (rev 513)
@@ -0,0 +1,52 @@
+##########################################
+## ##
+## Function for drawing IC ##
+## on observations ##
+## ##
+##########################################
+
+##@x - observations
+##@IC - 1-dim influence function (object)
+##@col - color of points
+##@bg - background for semitranparence
+##@cex - scaling of points
+
+
+drawPointsIC = function(x,IC,col=NULL,bg=NULL,cex=NULL,xlab=NULL,ylab=NULL){
+
+##evaluation of IC on observations
+IC0 = evalIC(IC,as.matrix(x))
+
+if (is.null(col)) col = rgb(192,192,192,maxColorValue=255,alpha=0.8)
+if (is.null(bg)) bg = c("#70707050")
+if (is.null(bg)) xlab = "x"
+if (is.null(ylab)) {ylab = numeric()
+for (i in 1:dim(IC0)[1]) ylab[i] = paste("IC of",names(L2fam at param@main[i]),"parameter")}
+
+
+par(mfrow=c(dim(IC0)[1],1),mar=c(2,4,3,2))
+
+for (i in 1:dim(IC0)[1]){
+
+##plotting IC
+plot(x,IC0[i,],ylab=ylab[i],xlab=xlab,'l',main = ifelse(i==1,paste(IC at name)," "),col="gray")
+
+##plotting observations onto IC
+points(x,IC0[i,]
+,pch = 21
+,cex = abs(IC0[i,]*sqrt(2))
+,col = col
+,bg = bg
+)
+
+grid()
+ }
+}
+
+##Example
+L2fam = GParetoFamily(shape = 0.7)
+IC = optIC(L2fam,risk=asCov())
+x = q(L2fam at distribution)(seq(0,1,length=50))
+
+##with probit transformation
+drawPointsIC(x,IC)
Added: branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R 2012-09-12 15:27:01 UTC (rev 513)
@@ -0,0 +1,61 @@
+##########################################
+## ##
+## Wrapper for outlyingnessPlot.R ##
+## ##
+## ##
+##########################################
+
+##projection distance
+ qfun = function(x){p0 = p(X)(x); q0 = q(X)(p0)}
+ QProj <- function(){new("NormType", name="Quantiles", fct=qfun)}
+
+##@x - dataset
+##@X - random variable
+##@fam - parameter family
+##@alpha - confidence level for quantile
+#
+plotOutlyingness = function(x,X=GPareto(),fam=GParetoFamily(),alpha=0.95){
+
+ ##logarithmic representation (for distributions with positive support)
+ fam at distribution = log(fam at distribution)
+
+ ##classical IC
+ ICmle <- optIC(model=fam,risk=asCov())
+
+ ##parameter for plotting
+ par(cex=1,bty="n")
+
+##call of routine from RobAStBase
+outlyingPlotIC(x
+ ,IC.x = ICmle
+ ,IC.y = ICmle
+ ,dist.x = QProj()
+ #NormType() - Euclidean norm, default - Mahalanobis norm
+ #,dist.y = NormType()
+ ,adj = 0.1
+ ,pch = 16
+ ,col = rgb(152,152,152,maxColorValue=255)
+ ,col.idn = rgb(102,102,102,maxColorValue=255)
+ ,cex.idn = 1.7
+ ,col.cutoff = rgb(202,202,202,maxColorValue=255)
+ ,offset = 0
+ ,cutoff.quantile.y = alpha
+ ,cutoff.quantile.x = alpha
+ ,cutoff.x = cutoff()
+ ,cutoff.y = cutoff.sememp()
+ ,robCov.x = TRUE
+ ,robCov.y = TRUE
+ ,tf.x = function(x)log(x)
+ ,cex.main = 1.5
+ ,cex.lab = 1.5
+ ,cex = 1.5
+ #,col.lab=FhGred
+ ,lwd.cutoff = 3
+ #,jitt.fac = 300
+ ,col.abline = rgb(102,102,102,maxColorValue=255)
+ ,cex.abline = 1.5
+ ,main = ""#"Outlyingness Plot"
+ ,xlab="Theoretical log-quantiles"
+ ,ylab="Mahalanobis distance"
+)
+}
Added: branches/robast-0.9/pkg/RobExtremes/R/plotScaledIC.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/plotScaledIC.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/plotScaledIC.R 2012-09-12 15:27:01 UTC (rev 513)
@@ -0,0 +1,157 @@
+##########################################
+## ##
+## IC Scaling of axes ##
+## ##
+##########################################
+
+
+##@IC - influence function
+##@L2fam - L2 parametric family
+##@probit - if probit transformation should be done, otherwise is a log-transformation of x-axis
+##@xlim, ylim - plot limits
+##@xticks, yticks - ticks for plot
+
+
+
+plotScaledIC = function(IC,L2fam,probit=TRUE,xlim=NULL,ylim=NULL,xticks=NULL,yticks=NULL){
+
+## get distribution from L2-diff family
+X = L2fam at distribution
+
+## get L2 derivative from L2-diff family
+L2derivative = L2fam at L2deriv.fct(L2fam at param)
+
+## dimension of derivative
+dim = length(L2derivative)
+
+xgrid = seq(1e-9,1-1e-9,length=1000)
+
+## quantile transformation of x-scale
+x = q(X)(xgrid)
+
+## evaluation of IC on quantiles
+ICeval = evalIC(IC,as.matrix(x))
+
+ICtr = ICeval
+
+## set plot limits in (0,1)-cub for probit-scaled IC
+if(probit){
+xlim = c(0,1)
+ylim = c(0,1)}
+
+minIC = min(ICtr)
+imin = which(ICtr == minIC)
+
+maxIC = max(ICtr)
+imax = which(ICtr == maxIC)
+
+## calculate reasonable plot limits for a non-transformed IC
+if ((is.null(ylim))&&(is.null(xlim))){
+
+if (!probit){
+
+if (dim == 1){
+ xmin = x[imin]
+ xmax = x[imax]
+}else{
+ dist = numeric(length(xgrid))
+for (i in 1:length(xgrid)){
+ sum = 0
+for (j in 1:dim){sum = ICtr[j,i]^2+sum}
+ dist[i] = sqrt(sum)
+ }
+}
+
+idist.max = which(dist == max(dist))
+idist.min = which(dist == min(dist))
+
+xmax.IC = x[idist.max]
+xmin.IC = x[idist.min]
+
+cat("\n IC takes its maximal values at",xmax.IC, "and minimal values at",xmin.IC,"\n")
+
+xmin = min(xmin.IC,xmax.IC)
+xmax = max(xmin.IC,xmax.IC)
+
+spanx = (xmax-xmin)/length(x)
+spany = (maxIC-minIC)/length(x)
+
+xlim = c(min(x),max(x))
+ylim = c(minIC,maxIC)
+ }
+}
+
+if (probit){
+## optimal min bias IC
+ robModel <- InfRobModel(center = L2fam, neighbor = ContNeighborhood(radius = 0.5))
+ ICmbr = optIC(model = robModel, risk = asBias())
+ sd = 2*ICmbr at clip
+
+## probit transformation of IC
+ for (i in 1:dim) ICtr[i,] = pnorm(ICeval[i,],sd=sd)
+}
+
+if(is.null(xticks)&&is.null(yticks)){
+## ticks for x,y-axis
+if (is.null(xticks)) xticks = c(-Inf,-1000,-100,-10,-5,-1,-0.5,-0.1,0.0,0.1,0.5,1,5,10,100,1000,Inf)
+if (is.null(yticks)) yticks = c(-Inf,-1000,-100,-10,-5,-1,-0.5,-0.1,0.0,0.1,0.5,1,5,10,100,1000,Inf)
+}
+
+if (probit){
+## transformation of ticks for x-axis
+ xax = round(p(X)(xticks),1)
+## probit transformation of ticks for y-axis
+ yax = round(pnorm(yticks,sd=sd),1)
+}else{
+ xax = round(q(X)(p(X)(xticks)),1)
+ yax = round(yticks,1)
+}
+
+## do ticks lie to near to each other?
+indy = !duplicated(yax)
+yax = yax[indy]
+yticks = yticks[indy]
+indx = !duplicated(xax)
+xax = xax[indx]
+xticks = xticks[indx]
+
+## function creating a grid on plot
+fun = function(){for(i in 1:length(yax)){
+ abline(v=xax[i],col="grey",lty="dashed",lwd=0.1)
+ abline(h=yax[i],col="grey",lty="dashed",lwd=0.1)
+ }
+}
+
+## plotting
+if (!probit){
+ plot(x, ICtr[1,], 'l', ylim=ylim, xlim=xlim, main=" ", xlab="x", ylab="IC", axes=FALSE, col=1,log="x")}
+else{
+ plot(xgrid, ICtr[1,], 'l', ylim=ylim, xlim=xlim, main=" ", xlab="x", ylab="probit(IC)", axes=FALSE, col=1)}
+
+##grid generation
+fun()
+
+## plotting all of the IC components
+if (dim > 1){
+ if (!probit){
+ for (i in 2:dim) lines(x, ICtr[i,], col = i)}
+ else{
+ for (i in 2:dim) lines(xgrid, ICtr[i,], col = i)
+ }
+}
+
+## drawing the ticks and labels on the plot
+axis(side=1, at= xax, labels = c(expression(paste(-infinity)),xticks[2:(length(xticks)-1)],expression(paste(infinity))), pos=0)
+axis(side=2, at= yax, labels = c(expression(paste(-infinity)),yticks[2:(length(yticks)-1)],expression(paste(infinity))))
+}
+
+##Example
+
+L2fam = GParetoFamily(shape = 0.7)
+IC = optIC(L2fam,risk=asCov())
+
+##with probit transformation
+plotScaledIC(IC,L2fam,probit = TRUE)
+
+##without probit transformation
+plotScaledIC(IC,L2fam,probit = FALSE,xlim=c(0.01,100),ylim=c(-10,10),yticks=c(-Inf,-10,-5,-1,0,1,5,10,Inf),xticks=c(0.01,1,5,10,100,Inf))
\ No newline at end of file
More information about the Robast-commits
mailing list