Hifime UH1 DAC under Linux and Jack

The Hifime UH1 USB DAC and Headphone Amp combo has some great features like:
- reference-quality DAC chip: Texas Instruments PCM5102 (a flagship chip used in other >$1000 DAC's)
- Samples rates from 44.1, 48, 88.2, 96, 176.4, 192, 384kHz, and bit rates 16, 24 and 32.
- optional headphone amplifier (TPA6120, which I couldn't get much information about), as well as a direct line-out just from the DAC itself. This a great feature for someone who may have to drive high-impedance headphones, but stills wants the flexibility of using other neutral amplifiers.
- Asynchronous USB transfer with two TCXO as master clock
- under $200 USB, ships from Hong Kong (not from the USA-based HifimeDIY)

What's really nice is that Hifime explicitly mentions support for Linux (without drivers!), a major selling point for me as I don't want to play around with kernel patches, etc.

Using with GNU/Linux JACK, Pulseaudio, and ALSA.
The DAC shows up as HiFimeDIY SA9227 USB Audio. For most linux distributions, one will have to contend with Pulseaudio to get any sound out of the DAC. The easiet way is to use 'pavucontrol' (lightweight Pulseaudio GUI control) to explicitly change the 'Audio Interfaces'.

Unfortunately, Pulseaudio automatically down-samples audio. So, for bit perfect audio, one needs to use a music player with ALSA support. I use gmusicbrowser, which is highly customizable (reminiscent of my only Win7 program Foobar2000). Click on Main (menu) -> Settings->Audio (Tab)-> Output device: set to "Alsa"; Click "advanced options" and enter "hw:Audio" as the Alsa option. Viola! Enjoy audio at 192000Hz and 24bit!

- See more at: http://hifimediy.com/hifime-uh1#sthash.twSvxSOI.dpuf

Bayesian Mark-Recapture: Parsimony and Pollock's Closed Robust Design

At the Bienniel Marine Mammal Conference 2016, I presented a Hierachical Bayesian version of Pollock's Closed Robust Design (see presentation here). This was motivated by some bad behaviour of Program Mark in estimating detection probabilities and migration parameters for "sparse data" -- I say sparse, but for my sample sizes with dolphins, failure is all too common. In both real data, and through an extensive simulation exercise, I found that Mark was highly prone to "boundary-level" MLEs for survival and migration and sometimes detection probabilities, leading to erratic abundance estimates.

Interestingly, mildy weak Bayesian priors (Beta(1,1)) help smooth over this bad behaviour ("Objective" priors do not).

The image to the right shows simulations of the PCRD, comparing Maximum Likelihood estimation versus a Bayesian model... notice the infinite CI's and boundary level MLEs for Program MARK! (See the full article link below)

Hierarchical Bayesian PCRD
The world is going Bayesian, but few people understand why. Here, the true power of Bayesian Inference is revealed in a hierarchical specification for the survival, detection and migration parameters: in particular, we shrink between "time-invariant theta(.)" models and "time-varying theta(t)" models. I found that applying a particular type of shrinkage-inducing hyperprior, called the "Scaled Half Student-t distribution", lead to somewhat similar estimates as "model-averaging" over many MLE-based fixed effect models (run in Mark).

Hierarchical Bayes provides a type of automatic "Occam's Razor" (Berger 2004), in that when there is sparse/weak data, the hyprprior shrinks the model to a simple model, but when there is a lot of evidence in the data or a large sample size, the hyperprior is easily overwhelmed and the model takes on a lot of complexity. Reassuringly, the final estimates seem insensitive to exact hyper-parameter values. As shown in the following animation...



I decided to upload an R & JAGS tutorial on github: https://github.com/faraway1nspace/PCRD_JAGS_demo Users can run simple fixed effects Bayesian PCRD models, or they can run a Hierarchical Bayesian version that includes: i) Individual Heterogeneity in detection probabilities; and ii) recruitment processes; and iii) uses a hyperprior to apply a type of "multimodel" inference. All things seem to unify in HB.

JAGS code: Fixed effects Bayesian PCRD
Here JAGS code to run the simple fixed model of the PCRD: phi(.)gamma'(.)gamma''(t)p(t,s). This is just for interested readers. Serious users should consult the github demonstration
. This implements the PCRD as a Bayesian HMM. Anyone who uses this JAGS code should cite: "Rankin RW, Nicholson K, Allen SJ, Krützen M, Bejder L, Pollock KH. 2016. A full-capture Hierarchical Bayesian model of Pollock's Closed Robust Design and application to dolphins. Frontiers in Marine Science, in Press. DOI: 10.3389/fmars.2016.00025". The article is available here

model{
# priors
phi ~ dbeta(1,1) # apparent survival probability (adjust for ecological
gamma1 ~ dbeta(1.2,1,2) # temp migration: probabilty stays out conditional on being out
# loop through primary periods: parameters for detection probability pd_t, and gamma2_t
for(t_ in 1:T){ # T primary periods...
pd_a[t_] ~ dgamma(3,2) # hyperprior on detect probs
pd_b[t_] ~ dgamma(3,2) # hyperprior on detect probs
for(tt_ in 1:T2[t_]){ # loop through secondary periods...
pd[t_,tt_] ~ dbeta(pd_a[t_],pd_b[t_]) # secondard-period-level detection probs
}
gamma2[t_] ~ dbeta(1,1) # temp migration: prob migrates outside conditional on being inside
# recruitment process from 'eigenvector decomposition'
lambda[1,t_] <- (1-gamma1)/(gamma2[t_]-gamma1+1) # recruitment ratio, or long-term prob of being inside
lambda[2,t_] <- 1-lambda[1,t_] #
psi[t_] ~ dunif(0,1) # inclusion probability
# trmat: transition matrix for Markovian latent-states
# 1 =not yet in population;2=dead;3=offsite;4=onsite (only observable state)
# transition are from the column --> rows
# trmat[row,column,time] = [state at time=t_; state at time t_-1; time=t_]
trmat[1,1,t_] <- 1-psi[t_] # excluded from pop
trmat[2,1,t_] <- 0 # dead
trmat[3,1,t_] <- psi[t_]*lambda[2,t_] # inclusion into pop, outside study are
trmat[4,1,t_] <- psi[t_]*lambda[1,t_] # inclusion into pop, inside study area
trmat[1,2,t_]<- 0
trmat[2,2,t_]<- 1 # stay dead
trmat[3,2,t_]<- 0
trmat[4,2,t_]<- 0
trmat[1,3,t_]<- 0
trmat[2,3,t_]<- 1-phi # dies outside
trmat[3,3,t_]<- gamma1*phi # stays outside | outside
trmat[4,3,t_]<- (1-gamma1)*phi # reenters study area | outside
trmat[1,4,t_]<- 0 #
trmat[2,4,t_]<- 1-phi # dies inside
trmat[3,4,t_]<- gamma2[t_]*phi # leaves study area | inside
trmat[4,4,t_]<- (1-gamma2[t_])*phi # stays inside | inside
} # t_
# likelihood: loop through M individuals, both real and pseudo-individuals
for (i in 1:M){
#draw latent state at primary period 1:
# ... by definition, everyone starts in z=1 (not-in-population) at time=0
z[i,1]~ dcat(trmat[1:4,1,1]) # first z strictly excluded from pop
# likelihood for first primary period
for(tt_ in 1:T2[1]){ # loop through secondary periods
# Bernouilli process, conditional on z=4, otherwise no observation
y[i,tt_,1] ~ dbern(pd[1,tt_]*equals(z[i,1],4))
}
alive_i[i,1] <- step(z[i,1]-3) # count if i is alive or not
Nin_i[i,1] <- equals(z[i,1],4) # count if i is within study area
# loop through primary periods after 1st primary periods
for(t_ in 2:T){
# state process: draw z(t_) conditional on z(t_-1)
z[i,t_] ~ dcat(trmat[1:4, z[i,t_-1] , t_])
# likelihood: loop through secondary period observations
for(tt_ in 1:T2[t_]){
# Bernouilli observation error, conditional on z=4
y[i,tt_,t_] ~ dbern(pd[t_,tt_]*equals(z[i,t_],4))
}
# tally population size
alive_i[i,t_] <- step(z[i,t_]-3) # check alive or not
Nin_i[i,t_] <- equals(z[i,t_],4) # count if i is within study area
} # t_
} # i
# estimate population size per primary periods
for(t_ in 1:T){
alive[t_] <- sum(alive_i[,t_]) # number alive
Nin[t_] <- sum(Nin_i[,t_]) # number in study area
} # t_
}


JAGS code: Hierarchical Bayesian PCRD
Here is JAGS code for a Hierarchical Bayesian Model, phi(t)gamma'(t)gamma''(t)p(t,s). The Hierarchical Part uses 'hyperpriors' to deliberately shrink time-varying specifications, like phi(t), to a simpler time-invariant specifications, phi(.). The result is a type of "model parsimony", with a similar type of smoothing as "model averaging". The following is just a visual demo for the casual reader; serious users should refer to the github page and download the full R code. Anyone who uses this JAGS code should cite: "Rankin RW, Nicholson K, Allen SJ, Krützen M, Bejder L, Pollock KH. 2016. A full-capture Hierarchical Bayesian model of Pollock's Closed Robust Design and application to dolphins. Frontiers in Marine Science, in Press. DOI:10.3389/fmars.2016.00025". The article is available here
model{
#hyperpriors: logit(theta_mu)~Normal and theta_sd ~ half-t
phi.mu ~ dunif(pr.phiunif[1],pr.phiunif[2]) #mean survival with a Uniform prior
sigma.phi ~ dt(0,pr.tauphi[1],pr.tauphi[2]) T(0,) #mean survival dispersion, half-t hyperprior
g1.mu ~ dnorm(0,pr.g1mu) #mean gamma1, temp migration out | out
sigma.g1~dt(0,pr.taug1[1],pr.taug1[2]) T(0,) #mean gamma1 dispersion, half-t hyperprior
g2.mu ~ dnorm(0,pr.g2mu) #mean gamma2, temp migration out | in
sigma.g2~dt(0,pr.taug2[1],pr.taug2[2]) T(0,) #mean gamma2 dispersion, half-t hyperprior
pd.mu ~ dnorm(0,pr.pdmu) #mean detection prob, overall
sigma.pdmu~dt(0,pr.taupdmu[1],pr.taupdmu[2]) T(0,) #primary period detection prob dispersion
sigma.pd2nd~dt(0,pr.taupd2nd[1],pr.taupd2nd[2]) T(0,) #secondary periods detection prob dispersion
sigma.eps ~ dt(0,pr.taueps[1],pr.taueps[2]) T(0,) #individual detection prob dispersion
#time-variant parameters
for(t_ in 1:T){ #loop through primary periods
pd_mu[t_]~dnorm(pd.mu,pow(sigma.pdmu,-2)) #primary period mean detaction prob (logit)
lgamma1[t_]~dnorm(g1.mu,pow(sigma.g1,-2)) #prob migrate out|out (logit)
gamma1[t_] <- ilogit(lgamma1[t_]) #prob migrate out|out (probability)
lgamma2[t_]~dnorm(g2.mu,pow(sigma.g2,-2)) #prob migrate out|in (logit)
gamma2[t_] <- ilogit(lgamma2[t_]) #prob migrate out|in (probability)
#RECRUITMENT: psi is the 'inclusion probability' and lambda is the 'recruitment ratio'
psi[t_]~dunif(0,1) #inclusion probability
lambda[t_] <- (1-gamma1[t_])/(gamma2[t_]-gamma1[t_]+1) #long-term probability inside study area
#NOTE, lambda could also be a random variable with a beta prior
#secondary-period detection probabilities
for(tt_ in 1:T2[t_]){ #loop through secondary periods
pd[t_,tt_] ~ dnorm(pd_mu[t_],pow(sigma.pd2nd,-2))
} #tt_
} #t_
#first state transition (pure nusance; strictly from outside-pop to part of marked-pop)
trmat0[1] <- (1-psi[1]) #remains not-yet-in-pop
trmat0[2] <- 0
trmat0[3] <- psi[1]*(1-lambda[1]) #inclusion into pop, goes outside study are
trmat0[4] <- psi[1]*lambda[1] #inclusion into pop, goes inside study
#state transitions (2:T)
for(t_ in 1:(T-1)){
lphi[t_]~dnorm(log(phi.mu/(1-phi.mu)), pow(sigma.phi,-2)) #survival prob (logit)
phi[t_]<-ilogit(lphi[t_])
#state transitions
#trmat: transition matrix for Markovian latent-states
#1 =not yet in population; 2=dead;3=offsite;4=onsite (only observable state)
#transition are from the column --> rows
#trmat[row,column,time] = [state at time=t_; state at time t_-1; time=t_]
#notice that the primary periods are offset by 1 (because we already dealt with T=1)
trmat[1,1,t_]<- 1-psi[t_+1] #excluded from pop
trmat[2,1,t_] <- 0 #dead
trmat[3,1,t_] <- psi[t_+1]*(1-lambda[t_+1]) #inclusion into pop,outside study
trmat[4,1,t_] <- psi[t_+1]*lambda[t_+1] #inclusion into pop,inside study
trmat[1,2,t_]<- 0
trmat[2,2,t_]<- 1 #stay dead
trmat[3,2,t_]<- 0
trmat[4,2,t_]<- 0
trmat[1,3,t_]<- 0
trmat[2,3,t_]<- 1-phi[t_] #dies outside
trmat[3,3,t_]<- gamma1[t_+1]*phi[t_] #stays outside | outside
trmat[4,3,t_]<- (1-gamma1[t_+1])*phi[t_] #reenters study area | outside
trmat[1,4,t_]<- 0 #
trmat[2,4,t_]<- 1-phi[t_] #dies inside
trmat[3,4,t_]<- gamma2[t_+1]*phi[t_] #leaves study area | inside
trmat[4,4,t_]<- (1-gamma2[t_+1])*phi[t_] #stays inside | inside
} #t_
#loop through M individuals
for (i in 1:M){
#state transitions and likelihood for the first primary period
z[i,1]~ dcat(trmat0) #z at time 0 is strictly 'not-yet-in-pop'
alive_i[i,1] <- step(z[i,1]-3) #count if i is alive or not
Nin_i[i,1] <- equals(z[i,1],4) #count if i is within study area
eps_i[i] ~ dnorm(0,pow(sigma.eps,-2)) #random effects at individual levels
#Observation error y[i,tt_,t_] ~ Bernoulli conditional on being inside z=4
for(tt_ in 1:T2[1]){ #loop through secondary periods
y[i,tt_,1] ~ dbern(equals(z[i,1],4)/(1+exp(-pd[1,tt_]-eps_i[i]))) #inverse-logit transform
}
#state transition and likelihood for primary periods 2:T
for(t_ in 2:T){
#State process: draw z(t_) conditional on z(t_-1)
z[i,t_] ~ dcat(trmat[1:4, z[i,t_-1] , t_-1])
#Observation error y[i,tt_,t_] ~ Bernoulli condition on being inside z=4
for(tt_ in 1:T2[t_]){ #loop through secondary periods
y[i,tt_,t_] ~ dbern(equals(z[i,t_],4)/(1+exp(-pd[t_,tt_]-eps_[i]))) #inverse-logit transform
}
#check whether i individual is alive and inside
alive_i[i,t_] <- step(z[i,t_]-3) #check alive
Nin_i[i,t_] <- equals(z[i,t_],4) #count if i is within study area
} #t_
} #i
#tally population size
for(t_ in 1:T){
alive[t_] <- sum(alive_i[,t_]) #number alive
Nin[t_] <- sum(Nin_i[,t_]) #number in study area
} #t_
} #end model


The Bandcamp scavenger hunt for jobs: API, cat, and wav

I was bored and jetlagged this Xmas, so I decided to follow Bandcamps job-hunting scavenger hunt: they lead prospective hires through several Unix-savvy tasks across the internet. I'd love to work in an applied musicology setting, but because I'm an ecologist, its highly unlikely they'd actually hire me. So, as pre-rejection revenge, I'm going to spoil the scavenger hunt! Well, not really, but its a great demostration for my students WHAT skills are needed in the massive-information economy, and how unconventional (and applied) their hiring process can be.

Task 1: source code
On the Bandcamp jobsite, prospects are asked to look at the Websites source code. Tool: wget
$ wget --output-document-file /tmp/bd.html https://bandcamp.com/jobs/#systems
... Upon downloading, one must open the file and scour the annotations for the next instruction...

Task 2: Bandcamp API
"Job hunters, use the Bandcamp API (use API key 'buithagsitannagalladrodlask') to get the url of band id 1 and go there. (Hint: not all API versions are the same."
Which is an annoying detour through the particulars of their API: its a simple URL-based information retreival system. Send a URL, receive JSON information. Answer: http://api.bandcamp.com/api/band/1/info?band_id=1&key=buithagsitannagalladrodlask which leads us to http://bandcamp.bandcamp.com for our next clue.

Task 3: ruby decrypting
http://bandcamp.bandcamp.com invites the job-hunter to download a music file, like any other Bandcamp download site. However, the output is not a FLAC, it is 7 large (seemingly) junk files in some scrambled ordered, variously called: "Bandcamp - You’re Getting Warm - part51.txt". Presumeably, the mystery files are something to be concatentated, in some order, the clue being a string of ruby code: "srand(42); 7.times {puts rand(100)}". This seemingly random operation is seeded, for a reproducible output. 51,92,14,71,60,20,82

Task 4: cat a wav
The ruby code spits out the order of the files. Now what? Lets use Unix 'cat' to make a new file. The header in the first file '51' is "RIFF": some quick internet-searching reveals that the intended format is a WAV audio file. So, we output the concatenation as a wav file and it plays. A synth-voice tells us where to send resumes and what is the secret-codeword to add so they don't delete our application.

Hurray! That's a lot of work for an ecologist. But what's a more impressive prospect? An ecologist who knows Linux and APIs, or a comp-sci student who is idle enough to follow a silly scavenger hunt?

remote connect R, SSH, EMACS, ESS + persistent sessions

I'm heading to Duke University on an exchange to study statistics, and hope to maintain access to my high-performance desktop at Murdoch University in Australia (aka, 'The Legend') which runs Ubuntu. After googling around a couple of options (and thanks in particular to this helpful exchange with Dirk Eddelbuettel at stackoverflow) I've set up a means to keep R processes running persistently over ssh, even after logging out. I thank the Great Speggetti Monster that I decided to learn EMACS! (see below for an alternative way to create persistent R sessions withouts using EMACS)

The following instructions assume that the remote computer is running in a Linux Environment, is available to the outside world with an IP address, has an ssh server installed and IS running (I installed openssh-server), and of course has a version of EMACS 24 installed.

On the local machine, open a terminal and enter the command:

ssh -X user@hostipaddress

... then go through prompts for passwords and passkeys.
On the first log-in, activate the emacs-daemon and start emacs buffers

emacs --daemon
emacsclient -nw

The daemon is what allows processes to keeping running after disconnection and to resume the session later when you log back in. The second line will open EMACS while displaying the buffers in the ssh terminal. One can proceed to open R scripts and start R sessions, etc., using familar commands suggest as M-x R-mode, etc.

Now, you can simply exit EMACS (C-x C-c) and logout of the remote computer (logout), but the buffers and processes will persist! To reconnect simply restart ssh access to the remote computer and redo the emacsclient -nw.

ssh -X user@hostipaddress
emacsclient -nw

... R should still be running :) For some reason, my C-space and C-S space keybindings do not work (though all other customized keybindings do), so be prepared to C-c C-n and C-c C-j.

Reconnecting a remote R session (without EMACS)

For those of you who don't use EMACS, I found tmux to serve as an easy alternative way to disconnect and reconnect to R sessions. While in a remote ssh session, one must simply start the tmux session by typing the command "tmux", then type "R" to initialize an R instance (This assumes the remote computer has installed tmux). When you want to disconnect the ssh connection, type "CTRL-b d" to exit tmux (while the R session nonetheless persists), then you're free to logout and close the remote connection. Later, when you want to re-connect to R, reestablish the ssh connection and type "tmux attach" to re-enter tmux. Your R session will still be there, chugging along as if nothing happened.

Planned contrasts and ANCOVA in Excel / LibreOffice Calc

I just finished TA'ing Advanced Quantitative Methods in Psychology at Murdoch University, and I'd like to share an Excel/Calc file (note its an .ods) to help students visualize ANCOVA/regression. Users can tweak regression coefficients to minimize the Sum of Square Residuals and/or find the lines of best fit. For the math enthusiasts, it also demos the underlying matrix multiplication needed to solve the problem analytically. The Experiment: we are interested in whether treatments of sucrose (8%, 32%, 64%) effect rats' speed at completing a maze (DV) while controlling for the rats age (aka, the 'covariate' in ANCOVA).

One should never 'do stats' in excel (so I also provide R code below), but this is a great, quick way to visualize several statistics concepts that many non-staticians may not understand, including:
- how a mix of continuous and categorical variables can be visualized and different parallel lines, - how categorical variables are represented as indicator variables in the model matrix. - how Planned Constrasts determine the type of indicator columns (the default in the .ods file is coded using the 'first' treatment as the reference treatment, like in R, but users can do their own. Try the Helmert!)
- how to solve least squares regression coefficients via matrix multiplication

If I were a teacher in an intro stats course, I'd ask students to reproduce the above. One operation I couldn't figure out was how to calculate the S.E. of the regression coefficients: this requires more matrix operations than Excel offers.

# ANCOVA in R
# data: DV is speed to complete maze
# categorical variables: 3 sucrose treatments
# covariate: rate AGE
library(car)
d<-data.frame(speed=c(1.4,2,3.2,1.4,2.3,4,5,4.7,6.2,3.1,3.2,4,4.5,6.4,4.4,4.1,5.8,6.6,6.5,5.9,5.9,3,5.9,5.6), trt=factor(rep(c(0,1,2),each=8),labels=c("sucrose8","sucrose32","sucrose64")), age=c(12,5,8,7,11,2,3,4,7,5,12,12,8,4,3,2,5,8,12,2,7,11,4,3))
# calculate regression coefficients and statistics manually
# then compare to worksheet and R
X <- model.matrix(~age+trt,data=d) # model matrix
betas<-solve(t(X)%*%X)%*%(t(X)%*%d$speed) # regression coefficients
SSm<-sum(((X%*%betas) - mean(d$speed))**2) # explained variation
SSr<- sum(((X%*%betas) - d$speed)**2) # residual sum of squares
df_m<-ncol(X)-1; df_r<-nrow(X)-ncol(X) # degrees of freedom
MSm<-SSm/df_m # Mean Square of Model
MSr <-SSr/df_r # means square error / estimate of variance
F.stat <- MSm/MSr

# standard error of regression coefficients
betas.se<- sqrt(diag(solve(t(X)%*%X)*MSr)) # MSr is an estimate of sigma^2
betas.t <- betas/betas.se # chance of seeing this large of effect ~ t distr
betas.p <- 2*(1-pt(q=abs(betas.t),df=df_r)) # p-values of regression coefs


# redo above, but with R's linear model function
m1 <- lm(speed~age+trt,d)
F.stat1<-summary(m1)$fstatistic # compare F.stat and F.stat1
coef(m1) # compare with betas
summary(m1)$coefficients # compare with betas.se,t,and p

# try with different Planned Constrasts:
# use Helmert contrast to test the hypothesis that 8% is different from
# 32% and 64% is different from both 8 and 32%
m.helm <- lm(speed~ age+trt, data=d, contrasts = list(trt = "contr.Helmert"))


lmms: powerful Industrial music studio

Last Tuesday I saw Nine Inch Nails live in Perth (with guest appearance by How to Destroy Angels!), and they inspired me to finally compose synth / industrial music. This is surprisingly possible when you have a GNU/Linux distro and the industry-grade software lmms, Linux MultiMedia Studio, something akin to the commercial Protools or FL Studio. You can engineer your own sounds and compose your own music with a large collection of instruments, effects, and synthesizers. For many years I have played around with fragments of melodies and drumbeats, but here is my first "song". With dose of creativity and OCD, this shows what you can do in a mere 24 hours.



For the uninitiated, the lingo and layout of lmms or Ardour3 or other Digital Audio workstations can be very intimidating. I hardly know where to start, nor how to give advice. Thanks to the rise of Dubstep, there are hundreds of very helpful youtube tutorials such as: http://www.youtube.com/watch?v=naQW2RkpfYY. I like this one because it is a nice visual intro to detuning and combing saw-tooth sound waves via an 'oscillator' to make a sound (instrument?) recognizable in, say, Skrillex.

Composing music with a DAW has revolutionized how I listen to songs: I can't help but listen for the subtle signatures of familiar saw-waves or moog-filters or peculiar sounds and then try to recreate them. It feels great, like becoming an active participant rather than just a passive consumer.

Diagnosing normality in R: QQ Plots and Shapiro-Wilk

I've become a teaching assistance for a 3rd year 'Stats for Psychologists' course in Australia. I thought I might share a little visualization to help my students intuit skew, normality, QQ-plots, and the Shapiro-Wilks test versus the Kolmogorov-Smirnov test. See the R code below the graphic output: I generate random data from a Gamma distribution (left) and a Normal distribution, to show how QQ-plots compare to the sample histograms, and under what sample sizes the tests of normality have any power to catch skew.

Highlights:
- as sample size decreases (heading down the graphic), both the KS-test and the Shapiro-Wilk loose their power to 'detect normality' (i.e., p values > 0.05)
- Shapiro-Wilk is more powerful than the more general KS-test.
- for QQ plots, one should focus on systematic trends along the inner-quantiles and ignore large deviations in the tails. The tails for the samples drawn from the Normal distribution (Right) can be wonky, but nonetheless the Shapiro-Wilk tests are not significant (as they shouldn't be).
- Subplotting in R is terrible (i.e., the nice histograms in the upper corners of the QQ plots): if you look at my code, most of the nastiest that detracts from the simple qqnorm and qqline functions is to make the inset histogram (on a positive note, I learned out to finely control where plots are positioned, rather than relying on the par(mfrow) command.
- Histograms are very difficult to interpret at low sample sizes. Herein lies the beauty of QQ-plots to show subtle patterns.
- I haven't shown a follow up simulation, but the KS-test is very sensitive to bimodal distributions of 2 mixtures of Normal distributions. The Shapiro-Wilk is not!

Left, different samples from a Gamma distribution (shape=5,scale=5), and Right, draws from a Normal distribution with the same mean and variance as the draws from the Gamma distribution.

Here is some R code to run the above simulation and get a feel for the QQ-plots and tests.

# SETUP SIMULATION PARAMETERS
N <- 100# sample size
shp_ <- 5; scale_ <- 5 # shape and scale parameters for Gamma distr
dat.skew <- as.numeric(scale(rgamma(n=N,shp_,scale_))) # random skewed data
dat.norm <- rnorm(N,mean(dat.skew),sqrt(var(dat.skew)))# random normal data
d <- list(skew=dat.skew, norm=dat.norm)
pltreg <- matrix(c(0,0.5,2/3,1, # x1,x2,y1,y2
0.5,1,2/3,1,
0,0.5,1/3,2/3,
0.5,1,1/3,2/3,
0,0.5,0,1/3,
0.5,1,0,1/3
),nrow=6,ncol=4,byrow=T)# arduous way to handle subplotting

# LOOP THROUGH DIFFERENT SAMPLE SIZES, PLOT 3x2
j <- 0 # counter to track which subplot we're on
invisible(par()) # make blank plot canvas
for(subn in c(1,0.3,0.15)){ # loop through different samples sizes
for(i in c("skew","norm")){ # loop through Gamma and Normal data
j <- j+1 # tracking plot region
sub.d <- sort(sample(d[[i]],round(N*subn))) # take subsample
par(fig=pltreg[j,],mar=c(4,4,1,1),new=TRUE*(j>1),mgp=c(1.2,0.5,0),bty="l")# pain in the ass way to position plots along the 3x2 canvas layout.
qqnorm(sub.d, datax=TRUE, col="forestgreen", pch=19,main=c(skew=paste("QQ (N=",round(N*subn)," from Gamma)",sep=""),norm=paste("QQ (N=",round(N*subn)," from Normal)",sep=""))[i]) # MAIN FUNCTION
qqline(sub.d,datax=TRUE) # 1:1 line
sw_p <- round(shapiro.test(sub.d)$p.value,3) # test of normality
ks_p <- round(ks.test(x=sub.d,y="pnorm")$p.value,3) # ''
text(0.8*diff(range(sub.d))+min(sub.d),y=0,labels=paste("N:",round(N*subn),"\nShapiro-Wilk:",sw_p,"\n","KS-test:",ks_p),pos=1)
par(fig=c(pltreg[j,c(1,1,4,4)]-c(0,-0.25,0.33/1.5,0)),new=TRUE)# setup inset histogram
hist(sub.d,col="grey",breaks=15,main="",xlab="",ylab="",xaxt="n",yaxt="n") # histogram
} # inner loop through two distributions
} # loop through sub-setting sample


Made using R v3, ESS and EMACS24 on Bodhilinux 2.4