postscript("l01.ps");options(width=64)
#*******************************************************/
# Mark.A: Densities standardized to the same IQR.      */
#*******************************************************/
source("common.R")
fun.comparedensityplot()
#*****************************************************/
# Mark.B: Simulate data from various                 */
# distributions and show size for normal theory test.*/
#*****************************************************/
library("VGAM")#vgam gives laplace distribution
library("BSDA")#Library gives z test
nsamp<-10000
nobsv<-c(10,40)
for(j in 1:2){
  size<-fun.comparepower(m=nobsv[j],nsamp=nsamp,
     dist=list("rnorm","rcauchy","rlaplace","runif"),
     hypoth=c(0,0,0,-.5),altd=c("two.sided","greater"))
 cat("\nSize for T and Sign Tests, sample size",nobsv[j],"\n")
 print(size)
}
#**********************************************/
# Mark.C: Why is Cauchy conservative? Examine */
# density of studentized Cauchy.              */
#**********************************************/
fun.studentizedcaucyplot(10,10000)
#Now calculate power for t test and sign test.
for(j in 1:2){
  power<-fun.comparepower(m=nobsv[j],nsamp=nsamp,
     dist=list("rnorm","rcauchy","rlaplace"),hypoth=1)
 cat("\nPower for T and Sign Tests, sample size",nobsv[j],"\n")
 print(power)
}
#***************************/;
# Mark.D: Arsenic Example. */;
#***************************/;
#*************************************************************/
# Data from http://lib.stat.cmu.edu/datasets/Arsenic         */
# reformatted into ASCII on the course home page.  Data re-  */
# flect arsenic levels in toenail clippings; covariates in-  */
# clude age, sex (1=M), categorical measures of quantities   */
# used for drinking and cooking, arsenic in the water, and   */
# arsenic in the nails.  To make arsenic.dat from Arsenic, do*/
#antiword Arsenic|awk '((NR>39)&&(NR<61)){print}'>arsenic.dat*/
# Potential threshold for ill health effects for toenails is */
# .26 http://www.efsa.europa.eu/de/scdocs/doc/1351.pdf       */
#*************************************************************/
arsenic<-as.data.frame(scan('arsenic.dat',
  what=list(age=0,sex=0,drink=0,cook=0,arswater=0,arsnails=0)))
library("BSDA")#Gives sign test.
cat("\n gn test and intervals\n")
#SIGN.test below operates strangely.  Most R functions return as
#a list all of their results, and then rely on a print function
#to format and print them.  SIGN.test prints them directly, and 
#saving the results saves confidence intervals, but #nothing else.
SIGN.test(arsenic$arsnails,md=.26,alternative="greater")
SIGN.test(arsenic$arsnails,md=.26)
attach(arsenic)
fun.invertsigntest<-function(dataset,alpha,maint){
  rd<-range(dataset)
  mu<-rd[1]-.2*diff(rd)+ (0:10000)/10000*1.4*diff(rd)
  statval<-apply(outer(dataset,mu,"-")<=0,2,sum)
  plot(mu,statval,type="l",
     main=paste("Construction of CI for",maint,"Location"))
  abline(h=a<-qbinom(alpha/2,length(dataset),.5),lty=2)
  abline(h=b<-qbinom(1-alpha/2,length(dataset),.5),lty=2)
  ci<-sort(dataset)[c(a,b+1)]
  for(i in 1:2) abline(v=ci[i],lty=3)
  return(ci)
}
fun.invertsigntest(log(arsnails),.05,"Log Nail Arsenic")