#!/bin/sh
# This is a shell archive (produced by GNU sharutils 4.2).
# To extract the files from this archive, save it to some FILE, remove
# everything before the `!/bin/sh' line above, then type `sh FILE'.
#
# Made on 1998-04-18 09:25 MDT by <jah@tundra>.
# Source directory was `/home/fac/a/jah/RA/gary/BMA.lib'.
#
# Existing files will *not* be overwritten unless `-c' is specified.
#
# This shar contains:
# length mode       name
# ------ ---------- ------------------------------------------
#   2889 -rw------- For.MC3.REG
#   5993 -rw------- MC3.REG
#   1542 -rw------- MC3.REG.choose
#   1747 -rw------- MC3.REG.logpost
#   6292 -rw------- README
#   1083 -rw------- out.lmsreg
#
save_IFS="${IFS}"
IFS="${IFS}:"
gettext_dir=FAILED
locale_dir=FAILED
first_param="$1"
for dir in $PATH
do
  if test "$gettext_dir" = FAILED && test -f $dir/gettext \
     && ($dir/gettext --version >/dev/null 2>&1)
  then
    set `$dir/gettext --version 2>&1`
    if test "$3" = GNU
    then
      gettext_dir=$dir
    fi
  fi
  if test "$locale_dir" = FAILED && test -f $dir/shar \
     && ($dir/shar --print-text-domain-dir >/dev/null 2>&1)
  then
    locale_dir=`$dir/shar --print-text-domain-dir`
  fi
done
IFS="$save_IFS"
if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED
then
  echo=echo
else
  TEXTDOMAINDIR=$locale_dir
  export TEXTDOMAINDIR
  TEXTDOMAIN=sharutils
  export TEXTDOMAIN
  echo="$gettext_dir/gettext -s"
fi
touch -am 1231235999 $$.touch >/dev/null 2>&1
if test ! -f 1231235999 && test -f $$.touch; then
  shar_touch=touch
else
  shar_touch=:
  echo
  $echo 'WARNING: not restoring timestamps.  Consider getting and'
  $echo "installing GNU \`touch', distributed in GNU File Utilities..."
  echo
fi
rm -f 1231235999 $$.touch
#
if mkdir _sh15493; then
  $echo 'x -' 'creating lock directory'
else
  $echo 'failed to create lock directory'
  exit 1
fi
# ============= For.MC3.REG ==============
if test -f 'For.MC3.REG' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'For.MC3.REG' '(file already exists)'
else
  $echo 'x -' extracting 'For.MC3.REG' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'For.MC3.REG' &&
For.MC3.REG<-function(i)
{
X
#Version 1.2: April 18, 1998
#S-plus 'For' loop program to compute posterior probabilities and 
#do MCMC for variable selection and outlier identification for the program "MC3.REG".
#Refer to the program "MC3.REG" for further description. 
X
# Developer: Jennifer Hoeting (jah@lamar.ColoState.EDU)
# with the assistance of Gary Gadbury (gadbury@statsun.stat.ColoState.EDU).
# This software is not formally maintained, but we will be happy to hear
# from people who have problems with it.
X
#Permission is hereby granted to StatLib to redistribute this software.
#The software can be freely used for non-commercial purposes, and can
#be freely distributed for non-commercial purposes.
#The copyright is retained by the developer.
#Copyright 1995 Jennifer A. Hoeting
X
#Inputs:
#   i -  a count parameter supplied by the program MC3.REG
X
#Outputs:
#   the updated matrix big.list is returned to the program MC3.REG
X
#Related programs:
#   refer to this section in the program "MC3.REG".
X
X	if (flag==1) 
X  	{
X           if (sum(M0.var)!=0) M0.1<<-sum(2^((0:(length(M0.var)-1))[M0.var]))+1
X           else M0.1<<-1
X      	   if (sum(M0.out)!=0) M0.2<<-sum(2^((0:(length(M0.out)-1))[M0.out]))+1
X      	   else M0.2<<-1
X        }
X
X
#Select a new candidate model 
X
X	M1 <- MC3.REG.choose(M0.var,M0.out)
X
X        if (sum(M1$var)!=0) M1.1<<-sum(2^((0:(length(M0.var)-1))[M1$var]))+1
X        else M1.1<-1
X        if (sum(M1$out)!=0) M1.2<<-sum(2^((0:(length(M0.out)-1))[M1$out]))+1
X        else M1.2<-1
X
#If the posterior has not been computed, it is computed in this section
X
X	if (sum(big.list[,1]==M1.1 & big.list[,2]==M1.2)==0)
X	{    
X           if (M1.1==1)
X              {
X              if (outcnt!=0) a<-(dim(Ys)[1]-sum(M1$out))*log(1-PI)+
X                  sum(M1$out)*log(PI)+ MC3.REG.logpost(Ys,Xs,0,0,outs.list[M1$out],
X                  K,nu,lambda,phi)
X              else a<-MC3.REG.logpost(Ys,Xs,0,0,outs.list[M1$out],K,nu,lambda,phi)
X              }   
X           else 
X              { 
X              if (outcnt!=0) a<-(dim(Ys)[1]-sum(M1$out))*log(1-PI)+
X              sum(M1$out)*log(PI)+MC3.REG.logpost(Ys,Xs,M1$var,
X              sum(M1$var),outs.list[M1$out],K,nu,lambda,phi)
X	      else a<-MC3.REG.logpost(Ys,Xs,M1$var,sum(M1$var),outs.list[M1$out],
X              K,nu,lambda,phi)
X              }   
X         big.list<<-rbind(big.list,c(M1.1,M1.2,a,0))        
X	 }
X
#BF, the Bayes factor is computed
X
X   	BF <- exp(big.list[big.list[,1]==M1.1 & big.list[,2]==M1.2,3]-
X		  big.list[big.list[,1]==M0.1 & big.list[,2]==M0.2,3])
X
X	if (BF >= 1) flag <<- 1     
X	else flag <<- rbinom(1,1,BF)
X
X        if (flag == 1) 
X	   {
X	      M0.var <<- M1$var 
X	      M0.out <<- M1$out
X   	      M0.1 <<- M1.1
X  	      M0.2 <<- M1.2
X	    }
X    	big.list[big.list[,1]==M0.1 & big.list[,2]==M0.2,4] <<-
X	 	big.list[big.list[,1]==M0.1 & big.list[,2]==M0.2,4]+1
}
X
X
X
X
X
X
X
X
X
X
X
X
SHAR_EOF
  $shar_touch -am 0418092198 'For.MC3.REG' &&
  chmod 0600 'For.MC3.REG' ||
  $echo 'restore of' 'For.MC3.REG' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'For.MC3.REG:' 'MD5 check failed'
b47766366bcb34172bc43c2cee45a6d4  For.MC3.REG
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'For.MC3.REG'`"
    test 2889 -eq "$shar_count" ||
    $echo 'For.MC3.REG:' 'original size' '2889,' 'current size' "$shar_count!"
  fi
fi
# ============= MC3.REG ==============
if test -f 'MC3.REG' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'MC3.REG' '(file already exists)'
else
  $echo 'x -' extracting 'MC3.REG' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'MC3.REG' &&
MC3.REG<-function(all.y,all.x,num.its,M0.var,M0.out,outs.list,PI,K,nu,lambda,phi)
{ 
X
#Version 1.2: April 18, 1998
#S-plus program to perform Bayesian simultaneous variable selection and
#outlier identification via Markov chain Monte Carlo model composition
X
#This program performs Bayesian simultaneous variable selection and
#outlier identification(SVO) via Markov chain Monte Carlo model
#composition(MC^3).  Papers on this subject are available from:
#http://www.colostate.edu/Depts/STAT/Documents/hoeting.papers.html
X
# Developer: Jennifer Hoeting (jah@lamar.ColoState.EDU)
# with the assistance of Gary Gadbury (gadbury@statsun.stat.ColoState.EDU).
# This software is not formally maintained, but we will be happy to hear
# from people who have problems with it.  
# The copyright is retained by the developer.
# Copyright 1995 Jennifer A. Hoeting
X
X
#Permission is hereby granted to StatLib to redistribute this software.
#The software can be freely used for non-commercial purposes, and can
#be freely distributed for non-commercial purposes.
X
#Inputs:
#   all.y -  the response matrix (1 column)
#   all.x -  the matrix of all possible covariates
#   num.its -  the number of iterations in the Markov chain
#   M0.var -  the starting model variable set in T/F format that is the
#  	same length as the number of predictors. For example, if you
#	have 3 predictors and the starting model is X_1 and X_3, then
#	M0.var would be c(T,F,T)
#	NOTE:  starting predictor model cannot be the null model.
#   M0.out -  the starting model outlier set in T/F format (see above)
#             this can be NULL only if outs.list is NULL, otherwise must be 
#             the same length as outs.list (but can be a vector of all "F's")
#   outs.list -  the list of all potential outlier locations 
#                (e.g. c(10,12) means the 10th and 12th points are
#                 potential outliers) - can be NULL  
#   PI -  a hyperparameter indicating the probability of an 
#         outlier (for Hastings ratio)
#   K -  a hyperparameter indicating the outlier inflation factor
#   nu, lambda, phi - regression hyperparameters
X
#Outputs:
#   model.matrix -  matrix of selected models (described below)
X
#Related programs:		
#   For.MC3.REG -  a For loop called by MC3.REG to compute posteriors and perform 
#                 MCMC for variable and outlier selection
#   MC3.REG.choose -  a function called by For.MC3.REG that randomly selects a
#                    neighboring model for posterior calculation
#   MC3.REG.logpost -  a function used by For.MC3.REG to compute the log of 
#                     the posterior for more than one outlier
#   out.lmsreg -  a function to select potential outliers from all.y.  The
#                 output of this function can be used as input to outs.list
X
#In this section the input is checked for consistency
X
X   if((sum(M0.var)==0)) 
X        stop("\nInput error:  M0.var cannot be null model")
X   if((length(M0.out)!=length(outs.list)) || (length(M0.var)!=dim(all.x)[2]))
X        stop("\nInput error: M0.*** is not the right length")
X
X
# The following permanent assignments are for the S-plus 'For' loop
# All computations are performed using standardized data
X   Ys <<- scale(all.y)
X   Xs <<- scale(all.x)
X 
X   M0.var<<-M0.var
X   M0.out<<-M0.out
X   outs.list <<- outs.list
X   K<<-K
X   nu<<-nu
X   PI<<-PI
X   lambda<<-lambda
X   phi<<-phi
X   flag <<- 1
X   outcnt<<-sum(outs.list)
X  
#   big.list is a matrix with:
#     column 1 = binary number for variables (corresponding to a binary
#                conversion of the M0.var, a T/F vector, to a number)
#     column 2 = binary number for outliers  
#     column 3 = posterior model probability
#     column 4 = # of times model has been visited
X
X   big.list<<-matrix(0,1,4)
X   big.list[1,1]<<-sum(2^((0:(length(M0.var)-1))[M0.var]))+1
X  if (sum(M0.out)!=0) big.list[1,2]<<-sum(2^((0:(length(M0.out)-1))[M0.out]))+1
X   else big.list[1,2]<<-1
X   if (outcnt!=0) big.list[1,3]<<-(dim(Ys)[1]-sum(M0.out))*log(1-PI)+
X                sum(M0.out)*log(PI)+MC3.REG.logpost(Ys,Xs,M0.var,
X                sum(M0.var),outs.list[M0.out],K,nu,lambda,phi)
X   else big.list[1,3]<<-MC3.REG.logpost(Ys,Xs,M0.var,
X                sum(M0.var),outs.list[M0.out],K,nu,lambda,phi)
X
X   For(i=1:num.its,For.MC3.REG(i))  
X
#The following section constructs the output matrix called model.matrix.
#  The rows show models for which posteriors were computed.
#  Columns 1 through length(M0.var) denote in T/F format the variables contained in a model.
#  Columns length(M0.var)+1 to length(M0.var)+length(outs.list) denote in T/F format which
#          outliers from out.list are contained in a model.
#  The final two columns are the posterior model probability(Post.Mod.Pr.)
#  and the number of times a model was visited(#visits), respectively. 
#  The rows are ordered by descending Post.Mod.Pr.
X   
X   var.vect<<-matrix(as.logical(rep(big.list[,1]-1,rep(length(M0.var),
X		length(big.list[,1])))
X              %/% 2^(0:(length(M0.var) - 1)) %% 2),ncol=length(M0.var),byrow=T)
X
X      n.var<-length(M0.var)
X      ndx<-1:n.var
X      Xn<-rep("X",n.var)
X      labs<-paste(Xn,ndx,sep="")
X
X      dimnames(var.vect)<-list(c(1:length(var.vect[,1])),labs)
X
X  postprob<<-matrix((exp(big.list[,3]))/(sum(exp(big.list[,3]))),ncol=1)
X
X      dimnames(postprob)[2]_list(c("Post.Mod.Pr."))
X
X
X   visits<<-matrix(big.list[,4],ncol=1)
X
X      dimnames(visits)[2]<-list(c("#visits"))
X
X
X   if (length(outs.list)!=0) {out.vect<<-matrix(as.logical(rep(big.list[,2]-1,
X                        rep(length(outs.list),length(big.list[,2])))
X              %/% 2^(0:(length(outs.list) - 1)) %% 2),ncol=length(outs.list),byrow=T)
X      dimnames(out.vect)<-list(c(1:length(out.vect[,1])),c(outs.list))
X   model.matrix<<-cbind(var.vect,out.vect,postprob,visits)}
X   else model.matrix<<-cbind(var.vect,postprob,visits)
X
X   colno<-length(M0.var)+length(M0.out)+1
X   model.matrix<<-model.matrix[order(-model.matrix[,colno]),]
X
X   return(model.matrix)
X
}
X
X
SHAR_EOF
  $shar_touch -am 0418092298 'MC3.REG' &&
  chmod 0600 'MC3.REG' ||
  $echo 'restore of' 'MC3.REG' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'MC3.REG:' 'MD5 check failed'
61fdd3d41d2a7d9c6c61743317f4d509  MC3.REG
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'MC3.REG'`"
    test 5993 -eq "$shar_count" ||
    $echo 'MC3.REG:' 'original size' '5993,' 'current size' "$shar_count!"
  fi
fi
# ============= MC3.REG.choose ==============
if test -f 'MC3.REG.choose' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'MC3.REG.choose' '(file already exists)'
else
  $echo 'x -' extracting 'MC3.REG.choose' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'MC3.REG.choose' &&
MC3.REG.choose<-function(M0.var,M0.out) {
X
X
#Version 1.2: April 18, 1998
#Splus program to select a new model for the program "MC3.REG".
#Refer to the program "MC3.REG" for further description.
X
# Developer: Jennifer Hoeting (jah@lamar.ColoState.EDU)
# with the assistance of Gary Gadbury (gadbury@statsun.stat.ColoState.EDU).
# This software is not formally maintained, but we will be happy to hear
# from people who have problems with it.
X
#Permission is hereby granted to StatLib to redistribute this software.
#The software can be freely used for non-commercial purposes, and can
#be freely distributed for non-commercial purposes.
#The copyright is retained by the developer.
#Copyright 1995 Jennifer A. Hoeting
X
#Inputs:
#   M0.var - current predictor model in T/F format supplied by the program "For.MC3.REG".
#   M0.out - current outlier model in T/F format supplied by the program "For.MC3.REG".
X
#Outputs:
#   The chosen model is returned to the program "For.MC3.REG".
X
#Related programs:
#   refer to this section in the program "MC3.REG".
X
X
X   var <- M0.var
X 
X   in.or.out <- sample(c(1:length(M0.var),rep(0,length(M0.out))),1)
#       0 means you stay at same predictors (and jump to new outlier)
X
X   if (in.or.out == 0){	#stay at same predictors, choose an outlier	
X      	out<-M0.out
X	in.or.out2 <- sample(1:length(M0.out),1)
X	out[in.or.out2]<-!M0.out[in.or.out2]
X	}
X   else {			
X        var[in.or.out]<- !M0.var[in.or.out]
X        out <- M0.out  #stay at same outliers, choose a predictor
X	}             
return(var,out) }
SHAR_EOF
  $shar_touch -am 0418092398 'MC3.REG.choose' &&
  chmod 0600 'MC3.REG.choose' ||
  $echo 'restore of' 'MC3.REG.choose' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'MC3.REG.choose:' 'MD5 check failed'
489a646da1dbb560c62c4ee834b247c0  MC3.REG.choose
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'MC3.REG.choose'`"
    test 1542 -eq "$shar_count" ||
    $echo 'MC3.REG.choose:' 'original size' '1542,' 'current size' "$shar_count!"
  fi
fi
# ============= MC3.REG.logpost ==============
if test -f 'MC3.REG.logpost' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'MC3.REG.logpost' '(file already exists)'
else
  $echo 'x -' extracting 'MC3.REG.logpost' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'MC3.REG.logpost' &&
MC3.REG.logpost<-function(Y, X, model.vect, p, i, K, nu, lambda, phi)
{
X
#Version 1.2: April 18, 1998
#Splus program to compute the log of the posterior for more than 1 outlier
#Refer to the program "MC3.REG" for further description.
X
# Developer: Jennifer Hoeting (jah@lamar.ColoState.EDU)
# with the assistance of Gary Gadbury (gadbury@statsun.stat.ColoState.EDU).
# This software is not formally maintained, but we will be happy to hear
# from people who have problems with it.
X
#Permission is hereby granted to StatLib to redistribute this software.
#The software can be freely used for non-commercial purposes, and can
#be freely distributed for non-commercial purposes.
#The copyright is retained by the developer.
#Copyright 1995 Jennifer A. Hoeting
X
#Inputs:
#   Y -  the column vector of responses
#   X -  the design matrix WITHOUT the column of ones
#   model.vect -  the vector of columns of X that should be included
#                 for the posterior calcuation
#   p -  the number of predictors
#   i -  the list of the outlier locations (i.e., 10th and 12th points
#        are outliers), in the form c(10,12)
#   K -  the value of the inflation factor for outliers
#   nu, lambda, phi - regression hyperparameters
X
#Outputs:
#   The log of the posterior for the current model
X
#Related programs:
#   refer to this section in the program "MC3.REG".
X
X	n <- dim(Y)[1]
X	ones <- rep(1, n)
X	A <- cbind(ones, X[, model.vect])
X	V <- diag(c(1, rep(phi^2, p)))
X	ones[i] <- K^2
X	det <- diag(ones) + A %*% V %*% t(A)
X	divs <- prod(eigen(det, T, T)$values)^0.5
X	denom <- (t(Y) %*% solve(det, Y)) + (nu * lambda)
X	lgamma((n + nu)/2) + log(nu * lambda) * (nu/2) - log(pi) * (n/2) - 
X		lgamma(nu/2) - log(divs) - ((nu + n)/2) * log(denom)
}
X
SHAR_EOF
  $shar_touch -am 0418092198 'MC3.REG.logpost' &&
  chmod 0600 'MC3.REG.logpost' ||
  $echo 'restore of' 'MC3.REG.logpost' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'MC3.REG.logpost:' 'MD5 check failed'
feb6bca8c67db8132dd0c5c6b684055e  MC3.REG.logpost
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'MC3.REG.logpost'`"
    test 1747 -eq "$shar_count" ||
    $echo 'MC3.REG.logpost:' 'original size' '1747,' 'current size' "$shar_count!"
  fi
fi
# ============= README ==============
if test -f 'README' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'README' '(file already exists)'
else
  $echo 'x -' extracting 'README' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'README' &&
README text file for the unix shar file BMA.shar.
X
X
1. All programs in BMA.shar have been developed by
X   Jennifer Hoeting (jah@lamar.ColoState.EDU), with the assistance of
X   Gary Gadbury (gadbury@statsun.stat.ColoState.EDU).
X
X   This software is not formally maintained, but we will be happy to hear
X   from people who have problems with it.
X
X   Permission is hereby granted to StatLib to redistribute this software.
X   The software can be freely used for non-commercial purposes, and can
X   be freely distributed for non-commercial purposes only.
X   The copyright is retained by the developer.
X   Copyright 1995 Jennifer A. Hoeting
X
2. References:
X
X  Raftery, A.E., Madigan, D., and Hoeting,
X  J.A. (1997) ``Bayesian Model Averaging for Linear Regression
X  Models'',   Journal of the American Statistical
X  Association, 92, 179-191.
X
X  Hoeting, J.A., Madigan, D., and Raftery, A.E. (1996)
X  ``A Method for Simultaneous Variable Selection and Outlier
X  Identification in Linear Regression'',
X  Journal of Computational Statistics and Data Analysis, 22, 251-270.
X
X  Hoeting, J. (1994), "Accounting for Model Uncertainty
X  in linear regression", Ph.D. dissertation, University of Washington.
X
X  These papers are available via the World Wide Web:
X   http://www.stat.colostate.edu/~jah
X
3. BMA(Bayesian Model Averaging).shar is a collection of S-plus programs
X   that perform Bayesian simultaneous variable selection and
X   outlier identification (SVO) via Markov chain Monte Carlo model
X   composition (MC^3).  There are five modules contained in BMA.shar in
X   addition to the README file.  They are:
X   MC3.REG - this is the parent program that receives input, structures
X             the data for calculations, and formats the output.  
X             See examples below.
X   For.MC3.REG - an S-plus For loop called by MC3.REG that computes
X                 posterior probabilities and performs MCMC for variable
X                 selection and outlier identification.
X   MC3.REG.choose - called by the For loop to randomly select a new predictor
X                    model or a new outlier model.
X   MC3.REG.logpost - called by both the parent program and the For loop to
X                     compute the log of the posterior model probability.
X   out.lmsreg - a stand alone program to select potential outliers from a data set.
X                The output of out.lmsreg can be used as input to the outs.list
X                parameter of MC3.REG.
X
4. The following examples show input commands for both out.lmsreg and
X   MC3.REG for selected data sets, and the resulting output (truncated).
X
X
Example 1:   Scottish hill racing data (see Hoeting 1994).
X
Sample input -
X
>b<-out.lmsreg(races[,-1],races[,1],2)
X
>RACES.run1<-MC3.REG(y.races,x.races,num.its=20000,c(F,T),rep(T,length(b)),
X            b,.1,7,.2,.1684,9.2)
X
>RACES.run1[1:20,]
X    X1 X2 6 7 10 11 14 15 17 18 19 26 33 35 Post.Mod.Pr. #visits
X 40  1  1 0 1  0  0  0  0  0  1  0  0  1  0  0.563009192   11234
X 34  1  1 0 1  0  0  0  0  0  1  1  0  1  0  0.108741574    2303
X 41  1  1 0 1  0  0  1  0  0  1  0  0  1  0  0.055648357    1302
X 57  1  1 1 1  0  0  0  0  0  1  0  0  1  0  0.051300315     947
X 68  1  1 0 1  0  0  0  0  0  1  0  0  0  0  0.042710082     697
X 47  1  1 0 1  0  1  0  0  0  1  0  0  1  0  0.018255390     372
X 42  1  1 0 1  0  0  0  0  0  1  0  1  1  0  0.017550885     246
X 70  1  1 0 1  0  0  0  1  0  1  0  0  1  0  0.014377033     388
X 66  1  1 0 1  1  0  0  0  0  1  0  0  1  0  0.012362058     204
X 33  1  1 1 1  0  0  0  0  0  1  1  0  1  0  0.011938033     174
X 67  1  1 0 1  0  0  0  0  0  1  0  0  1  1  0.011764571     224
X 36  1  1 0 1  0  0  1  0  0  1  1  0  1  0  0.011598439     243
X 69  1  1 0 1  0  0  0  0  1  1  0  0  1  0  0.009540701     232
106  1  1 1 1  0  0  1  0  0  1  0  0  1  0  0.005918905     125
X 25  1  1 0 1  0  0  0  0  0  1  1  1  1  0  0.004149772      69
X 39  1  1 0 1  0  1  0  0  0  1  1  0  1  0  0.003716219      76
109  1  1 0 1  0  0  0  1  0  1  1  0  1  0  0.003273965      52
X 37  1  1 0 1  0  0  0  0  0  1  1  0  0  0  0.003219818      64
X 77  1  1 0 1  0  0  1  0  0  1  0  0  0  0  0.003137438      56
X 35  1  1 0 1  0  0  0  0  0  1  1  0  1  1  0.002376023      55
X
---------------------------------------------------------------------------
X
Example 2:  Crime data (see Hoeting 1994).
X
In this example  we did not consider outliers, so the outlier model and
the parameter outs.list are set to NULL. Since outliers are not
considered,  hyperparameters PI and K are also set to zero.
X
Note:   The program accepts a starting model for outliers set equal to
X        NULL only if outs.list (the set of potential outliers) is also
X        NULL.
X
X
Sample input/output:
X
>crime.run1<-MC3.REG(y.crime.log,x.crime.log,30000,rep(T,15),NULL,NULL,0,0,
X        nu= 2.58, lambda = 0.28,phi = 2.85)
X
>crime.run1[1:20,]
X     X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 Post.Mod.Pr. #visits
X 167  1  0  1  1  0  0  0  0  1   0   1   0   1   1   0  0.024538963     760
X 458  1  0  1  1  0  0  0  0  0   0   1   0   1   1   0  0.017479142     495
X 144  1  0  1  1  0  0  0  0  1   0   0   0   1   1   0  0.016261682     462
X 340  1  0  1  0  1  0  0  0  1   0   1   0   1   1   0  0.015532250     586
X 168  1  0  1  1  0  0  0  0  1   0   1   0   1   1   1  0.015034825     432
X 142  1  0  1  1  0  0  0  0  1   0   0   0   1   1   1  0.015019941     318
X  37  0  0  1  1  0  0  0  1  1   0   0   0   1   1   0  0.014777818     460
X  74  1  0  1  1  0  0  0  0  0   0   0   0   1   1   0  0.012303096     408
X 846  1  0  1  1  0  0  0  0  0   0   1   0   1   0   0  0.011384699     293
X 295  1  0  1  0  1  0  0  0  0   0   1   0   1   1   0  0.011137065     304
X 283  1  0  1  1  1  0  0  0  1   0   1   0   1   1   0  0.010916401     302
X 457  1  0  1  1  0  0  0  0  0   0   0   0   1   0   0  0.009489289     332
X 194  1  0  1  1  0  0  0  0  1   0   1   1   1   1   0  0.009409087     307
X 339  1  0  1  0  1  0  0  0  1   0   0   0   1   1   0  0.009305826     326
X 175  1  0  1  1  0  0  0  1  1   0   0   0   1   1   0  0.009001765     240
1062  1  0  1  1  1  0  0  0  0   0   1   0   1   1   0  0.008725018     273
X 396  0  0  1  0  1  0  0  1  1   0   0   0   1   1   0  0.008660988     219
X 136  1  0  1  1  0  0  0  0  1   0   0   1   1   1   1  0.008474762     171
X 116  1  0  1  0  1  0  0  0  1   0   1   1   1   1   0  0.008177776     307
X 237  0  0  1  1  0  0  0  0  1   0   0   0   1   1   0  0.007964352     302
X
X
SHAR_EOF
  $shar_touch -am 0418092498 'README' &&
  chmod 0600 'README' ||
  $echo 'restore of' 'README' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'README:' 'MD5 check failed'
1b7ca8d21b0a078f4c0ed19b166e420e  README
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'README'`"
    test 6292 -eq "$shar_count" ||
    $echo 'README:' 'original size' '6292,' 'current size' "$shar_count!"
  fi
fi
# ============= out.lmsreg ==============
if test -f 'out.lmsreg' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'out.lmsreg' '(file already exists)'
else
  $echo 'x -' extracting 'out.lmsreg' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'out.lmsreg' &&
out.lmsreg <- function(x,y,delta)  {
X
#Version 1.2: April 18, 1998
#Splus function to identify potential outliers
X
# Developer: Jennifer Hoeting (jah@lamar.ColoState.EDU)
# with the assistance of Gary Gadbury (gadbury@statsun.stat.ColoState.EDU).
# This software is not formally maintained, but we will be happy to hear
# from people who have problems with it.
X 
#Permission is hereby granted to StatLib to redistribute this software.
#The software can be freely used for non-commercial purposes, and can
#be freely distributed for non-commercial purposes.
#The copyright is retained by the developer.
#Copyright 1995 Jennifer A. Hoeting
X 
#Inputs:
#   x - the design matrix
#   y - the observations
#   delta -  the threshold set by the user.  Standardized residuals
#            from LMS regression that are larger than delta are identified as
#            potential outliers
X
#Output:
#   The result of this procedure can be used to supply the outs.list parameter
#   of the program "MC3.REG".
X
X
X   abc<-lmsreg(x,y)$resid
X   (1:length(y))[abs(as.vector(abc)/mad(abc))>=delta]}
SHAR_EOF
  $shar_touch -am 0418092398 'out.lmsreg' &&
  chmod 0600 'out.lmsreg' ||
  $echo 'restore of' 'out.lmsreg' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'out.lmsreg:' 'MD5 check failed'
56faa1a8b68757f4152256232b87989a  out.lmsreg
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'out.lmsreg'`"
    test 1083 -eq "$shar_count" ||
    $echo 'out.lmsreg:' 'original size' '1083,' 'current size' "$shar_count!"
  fi
fi
rm -fr _sh15493
exit 0
