panel.dist_summary <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 3) )
    #hist(x, probability=T, add=T)
    
    #default color and line style for density plot
    density.col = 'gray'
    density.lty = 3
    
    # is this a normally distributed dataset?
    # if so, change the color of the density plot
    # The test rejects the null hypothesis if W is too small.
    s.W <- shapiro.test(x)$statistic
    if( (s.W > 0.91) == TRUE)
    	{
    	density.col = 'gray'
    	density.lty = 1
    	}
    
    # compute and plot density
    d <- density(x)
    dy <- d$y / max(d$y) * .5
    lines(d$x, dy, col=density.col, lty=density.lty)
    
    # get a small increment to use in the next tests:
    delta <- abs(min(x) - max(x)) / 100
    
    y_mean <- dy[d$x < mean(x) + delta & d$x > mean(x) - delta][1]
    y_median <- dy[d$x < median(x) + delta & d$x > median(x) - delta][1]
    
    debug
    #print(y_median)
    
    #add points on the density plot for the mean and median
    points( c(mean(x), median(x)), c(y_mean, y_median), col=c('red', 'orange'), pch=16)
    
    #add a boxplot
    boxplot(x, horizontal=TRUE, boxwex=0.3, add=T)
    
    #debugging
    #print(s.W)
}


panel.cor <- function(x, y, digits=2, prefix="", cex.cor, cor.method="pearson")
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y, method=cor.method))
    txt <- format(c(r, 0.123456789), digits=digits)[1]
    txt <- paste(prefix, txt, sep="")
    if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex * r, col='gray')
    
    # might be interesting to use ks.test
    # http://www.physics.csbsju.edu/stats/KS-test.html
}

# using built-in data sets

par(pty='s')

pairs(USJudgeRatings[1:5], upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)

pairs(iris[1:4], upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)

# use spearman correlation calculation instead of default person:
pairs(iris[1:4], upper.panel=function(...) panel.cor(..., cor.method="spearman"), lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)

# color iris specis: note location of 'col=' argument
pairs(iris[1:4], upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1, pch=16, col=c("red4", "green3", "blue4")[unclass(iris$Species)]), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2 )


pairs(trees, upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)

pairs(swiss, upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., 
col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)

# using formula notation:
pairs( ~ Fertility + Education + Catholic, data=swiss, upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)

pairs(longley, upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., 
col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)



# jocelyn's soil data
pairs(soil, upper.panel=panel.cor, lower.panel=function(...) panel.smooth(..., col.smooth=gray(.5), lty=1), diag.panel=panel.dist_summary, cex.labels = 2, font.labels=2)