## FIXME: could shorten with percentile info only?
babydat0 <- read.table(header=TRUE,text="
bucket percentile
15	0.2	
16	1.1	
17	3.9	
18	10.2	
19	21.5	
20	36.8	
21	54.1	
22	69.9	
23	82.2	
24	90.3	
25	95.3	
26	97.9	
27	99.1	
28	99.7	
29	99.9	
30	100	
")

## differences of cumulative curve
babydat <- transform(babydat0, data=diff(c(0,percentile)))
## mean, std dev
m <- with(babydat, sum(bucket*data)/sum(data))
s <- with(babydat, sqrt(sum((bucket-m)^2*data)/sum(data)))

pdf("baby_bar.pdf",width=6.1,height=4.2)
par(las=1,bty="l",mgp=c(2,1,0))
b <- with(babydat,barplot(data,space=0,
                          ylab="Percentage of babies",
                          xlab="Weight (pounds)"))
mtext(text=babydat$bucket,at=b,side=1,cex=0.8)
off <- babydat$bucket[1]-b[1]
curve(100*dnorm(x+off,mean=m,sd=s),add=TRUE,lwd=2)
dev.off()

pdf("norm_curves.pdf",width=6.1,height=4.2)
par(las=1,bty="l",mgp=c(2,1,0),yaxs="i")
sdvec <- c(4.8,2.4,1.2)
curve(100*dnorm(x,mean=m,sd=sdvec[1]),from=5,to=37,ylim=c(0,35),
      xlab="",ylab="",lwd=2)
curve(100*dnorm(x,mean=m,sd=sdvec[2]),add=TRUE,lty=2,lwd=2)
curve(100*dnorm(x,mean=m,sd=sdvec[3]),n=201,add=TRUE,lty=3,lwd=2)
xpos <- c(22.7, 24.3, 28.2)-1.2
ypos <- c(29.4,11.2, 4.5)
for (i in 1:3) {
    text(xpos[i],ypos[i], pos=4,
         substitute(~sigma==s,
                    ## expression(sigma==s),
                    list(s=sdvec[i]))
         )
}
dev.off()