#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'README' <<'END_OF_FILE' X2-12-93 X XNOTE: This is a shareware file for modal regression, X which generates the penny thickness example in X Scott, D.W. (1992). Multivariate Density Estimation: X Theory, Practice, and Visualization, John Wiley & Sons. X pp. 233-235. X X1. create a new directory and a ".Data" directory there X2. unpack the file modalreg.shar by typing "sh modalreg.shar" X3. the function "modalreg" is in "ash.s" X4. source("modalreg.eg") handles all compiling, sourcing, and X plotting X5. the use of modalreg() should be clear afterwards X XDavid W. Scott END_OF_FILE if test 542 -ne `wc -c <'README'`; then echo shar: \"'README'\" unpacked with wrong size! fi # end of 'README' fi if test -f 'modalreg.eg' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'modalreg.eg'\" else echo shar: Extracting \"'modalreg.eg'\" \(140 characters\) sed "s/^X//" >'modalreg.eg' <<'END_OF_FILE' X# 2-12-93 X X# example of modal regression for penny data X Xsource("ash.s") X Xunix("f77 -c ash.f") X Xsource("penny.plot") X Xsource("penny5.plot") END_OF_FILE if test 140 -ne `wc -c <'modalreg.eg'`; then echo shar: \"'modalreg.eg'\" unpacked with wrong size! fi # end of 'modalreg.eg' fi if test -f 'ash.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ash.s'\" else echo shar: Extracting \"'ash.s'\" \(9217 characters\) sed "s/^X//" >'ash.s' <<'END_OF_FILE' X## ash.s---------------------------------------------------------------------- X X##### Form univariate bin counts X##### Copyright 1986,1990 David W. Scott X X# x data vector X# ab half-open bin range [a,b) X# nbin number of bins X# opt optional message control T=suppress F=default X X# nc bin counts X# nskip number of points outside interval [a,b) X Xbin1 <- function(x,ab=nicerange(x),nbin=100,opt=FALSE){ X X if(ab[1]>=ab[2]) fatal("Interval vector has negative orientation") X if(nbin<=0) fatal("Number of bin intervals nonpositive") X a <- .Fortran("bin1", X as.single(x), X as.integer(length(x)), X as.single(ab), X as.integer(nbin), X nc=integer(nbin), X nskip=integer(1)) X if( !opt & a$nskip>0) warning( X paste(as.character(a$nskip),"points not in open interval ab")) X return( list(nc=a$nc,ab=ab,nskip=a$nskip) ) X} X X#### S function to compute one-dimensional averaged shifted histogram X#### Copyright 1986,1990 David W. Scott X X# nc vector of bin counts X# n assumed to be total of "nc" X# nbin assumed to be length of "nc" X# ab half-open interval containing bins [a,b) X# opt optional message control T=suppress F=default X X# t vector of bin centers X# f vector of ash density estimates X# work work vector of ash weights X Xash1 <- function(bins,m=5,opt=FALSE){ X if(m<=0) stop("Smoothing parameter m nonpositive") X nc <- bins$nc X ab <- bins$ab X nbin <- length(nc) X a <- .Fortran("ash1", X as.integer(m), X as.integer(nc), X as.integer(nbin), X as.single(ab), X t=single(nbin), X f=single(nbin), X work=single(m), X ier=single(1)) X X# optional warning if density nonzero outside ab X if( !opt & a$ier!=0) warning("Density nonzero outside interval ab") X X return(list(x=a$t,y=a$f,ab=ab,m=m,n=sum(nc),nskip=bins$nskip)) X} X X#### S function to evaluate univariate ASH X#### Copyright 1990 David W. Scott X Xfash1 <- function(x,fhat) { X tk <- fhat$x # bin midpoints X fk <- fhat$y # ash values X nt <- length(tk) # same as nbin X nx <- length(x) X fx <- rep(0,nx) X h <- tk[2]-tk[1] X ptr <- (x-tk[1])/h + 1 X ip <- floor(ptr) X rp <- ptr-ip X i <- seq(nx)[ ptr>0 & ptr<=1 ] X if(length(i)>0) fx[i] <- rp[i]*fk[1] X i <- seq(nx)[ ptr>nt & ptr<=nt+1 ] X if(length(i)>0) fx[i] <- (1-rp[i])*fk[nt] X i <- seq(nx)[ ptr>1 & ptr<=nt ] X if(length(i)>0) fx[i] <- (1-rp[i])*fk[ip[i]] + rp[i]*fk[ip[i]+1] X return(list(x=x,y=fx)) X} X X#### S function to add kernel bump to plot X#### Copyright 1990 David W. Scott X Xaddbump <- function(fhat,x=fhat$x[length(fhat$x)/2],base=FALSE) { X a <- ash1( bin1(x,fhat$ab,length(fhat$x)), fhat$m) X i <- seq(length(a$x))[a$y>0] X i <- c( i[1]-1, i, i[length(i)]+1 ) # padd each side X lines(a$x[i],a$y[i]/fhat$n) X if(base) lines(a$x[i[c(1,length(i))]],c(0,0)) # little baseline X} X X##### S function to form bivariate bin counts X##### Copyright 1986,1990 David W. Scott X X# x data matrix - n by 2 X# ab bin support rectangle (half-open) X# x-axis [ ab[1,1] , ab[1,2] ) X# y-axis [ ab[2,1] , ab[2,2] ) X# nbin number of x and y axis bins X# opt optional message control T=suppress F=default X X# nc bin count matrix (see "persp" function) X# nskip number of points outside binning rectange X Xbin2 <- function(x,ab,nbin=c(40,40),opt=FALSE){ X X if( !is.matrix(x) || ncol(x)!=2 ) stop("x must be a matrice with 2 columns") X if(missing(ab)) ab <- matrix(c(nicerange(x[,1]),nicerange(x[,2])),2,2,byrow=T) X else { if( !is.matrix(ab) || any(dim(ab)!=c(2,2)) ) X stop("ab must be a 2x2 matrix") X if ( ab[1,1]>=ab[1,2] | ab[2,1]>=ab[2,2] ) X stop("Array ab contains intervals with negative orientation") } X X if ( !opt ) { X print(paste("x axis binned over interval",as.character(signif(ab[1,1],4)), X "to",as.character(signif(ab[1,2],4)))) X print(paste("y axis binned over interval",as.character(signif(ab[2,1],4)), X "to",as.character(signif(ab[2,2],4)))) } X X if ( nbin[1] <= 0 | nbin[2] <= 0 ) X stop("Number of bins along x or y axis nonpositive") X X a <- .Fortran("bin2", X as.single(x), X as.integer(nrow(x)), X as.single(ab), X as.integer(nbin), X nc=matrix(as.integer(0),nbin[1],nbin[2]), X nskip=integer(1)) X X if( !opt & a$nskip>0 ) X print(paste(as.character(a$nskip),"points not in open rectangle ab")) X return( list(nc=a$nc,ab=ab,nskip=a$nskip) ) X} X X##### S function to compute bivariate ash estimates X##### Copyright 1986,1990 David W. Scott X X# nc bivariate bin count matrix X# n assumed to be total of "nc" and X# ## bins == ## rows & cols of "nc" X# ab bin support rectangle X# m ash x and y axis smoothing parameters X# opt optional message control T=suppress F=default X X# x length nbin[1]; bin centers z=z(x,y) X# y length nbin[2] X# z bivariate ash estimates (see "persp/contour" functions) X# work bivariate ash smoothing weights (first quadrant) X# ier indicates whether f nonzero outside bin rectangle X Xash2 <- function(bins,m=c(5,5),opt=FALSE){ X if(min(m)<=0) stop("Smoothing parameters m[1] or m[2] nonpositive") X nbin <- dim(bins$nc) X ab <- bins$ab X a <- .Fortran("ash2", X as.integer(m), X as.integer(bins$nc), X as.integer(nbin[1]), X as.integer(nbin[2]), X as.single(ab), X f=matrix(as.single(0),nbin[1],nbin[2]), X work=matrix(as.single(0),m[1],m[2]), X ier=single(1)) X X if( !opt & a$ier!=0 ) warning("Density nonzero outside rectangle ab") X return( list(x=bincen(ab[1,],nbin[1]),y=bincen(ab[2,],nbin[2]), X z=a$f,v=seq(0,max(a$f),l=11)[-1],m=m) ) X} X X#### S function to evaluate bivariate ASH X#### Copyright 1990 David W. Scott X Xfash2 <- function(x,fhat) { X if(!is.matrix(x) || ncol(x)!=2 ) stop("x in fash2 not matrix with 2 columns") X txk <- fhat$x # bin midpoints X tyk <- fhat$y # bin midpoints X fk <- c(fhat$z) # ash values -- string out as vector for indexing X ntx <- length(txk) # same as nbin[1] X nty <- length(tyk) # same as nbin[2] X nx <- nrow(x) X fxy <- rep(0,nx) # zero outside bin center support X hx <- txk[2]-txk[1] X hy <- tyk[2]-tyk[1] X ptrx <- (x[,1]-txk[1])/hx + 1; ipx <- floor(ptrx); rpx <- ptrx-ipx X ptry <- (x[,2]-tyk[1])/hy + 1; ipy <- floor(ptry); rpy <- ptry-ipy X X#### extrapolate ash2 out half a bin width to original "ab" support rectangle X k_seq(ipx)[ipx== 0&rpx>.5]; if(length(k)>0) {ipx[k]_ 1;rpx[k]_rpx[k]-1} X k_seq(ipx)[ipx==ntx&rpx<.5]; if(length(k)>0) {ipx[k]_ntx-1;rpx[k]_rpx[k]+1} X k_seq(ipy)[ipy== 0&rpy>.5]; if(length(k)>0) {ipy[k]_ 1;rpy[k]_rpy[k]-1} X k_seq(ipy)[ipy==nty&rpy<.5]; if(length(k)>0) {ipy[k]_nty-1;rpy[k]_rpy[k]+1} X X i <- seq(nx)[ ipx>=1 & ipx=1 & ipy=a[i-1] & a[i]>=a[i+1] ] X X if(length(j)>0 & !quadopt) { # linear blend maximum X k <- j[a[j]>up[2]] X if(length(k)>0) points(rep(xx[x],length(k)),yy[k],pch=ch) X k <- j[ a[j]<=up[2] & a[j]>=up[1]] X if(length(k)>0) points(rep(xx[x],length(k)),yy[k],pch="x") X} X if(length(j)>0 & quadopt) { # local quadratic fit X yp[j] <- a[j-1]-2*a[j]+a[j+1] # second difference at maximum X jj <- j[ yp[j]>=0 ] # singular case X if(length(jj)>0) yp[jj] <- yy[jj] # use middle point X jj <- j[ yp[j]<0 ] # find local maximum [-.5,.5] X if(length(jj)>0) yp[jj] <- yy[jj] + dy*0.5*(a[jj-1]-a[jj+1])/yp[jj] X k <- j[a[j]>up[2]] X if(length(k)>0) points(rep(xx[x],length(k)),yp[k],pch=ch) X k <- j[ a[j]<=up[2] & a[j]>=up[1]] X if(length(k)>0) points(rep(xx[x],length(k)),yp[k],pch="x") X} } } X X X#### S function to compute bivariate ASH regression line(s) X#### Copyright 1990 David W. Scott X Xash2reg <- function(fhat) { X x <- fhat$x; y <- fhat$y X nx <- length(x) X ny <- length(y) X rx <- rep(NA,length(x)) X for( i in seq(x) ){ X fyx <- fhat$z[i,] X sfyx <- sum(fyx) X if(sfyx>0) rx[i] <- sum(y*fyx)/sfyx X lines( x[1:i],rx[1:i]) X } X Xbrowser() X X x <- fhat$x; y <- fhat$y X dy <- y[2]-y[1]; ym <- c(y[1]-dy,y)+dy/2 # midpoints, average adjacent y's X rx2 <- rep(NA,length(x)) X X for( i in seq(x) ){ X fyx <- fhat$z[i,] X fmyx <- ( c(0,fyx) + c(fyx,0) ) /2 # average adjacent fyx's X sfmyx <- sum(fmyx) X if(sfmyx>0) rx2[i] <- sum(ym*fmyx)/sfmyx X lines( x[1:i],rx2[1:i]) X } Xbrowser() X return(list(x=x,y=rx2,m=fhat$m)) X} X Xnicerange <- function(x, beta = 0.2) { X ab <- range(x) # interval surrounding data X del <- ((ab[2] - ab[1]) * beta)/2 X return(c(ab + c( - del, del))) } X Xbincen <- function(ab,nbin) { # sequence of bin centers X del <- diff(ab)/nbin X x <- seq(ab[1]+del/2,by=del,length=nbin) X x } END_OF_FILE if test 9217 -ne `wc -c <'ash.s'`; then echo shar: \"'ash.s'\" unpacked with wrong size! fi # end of 'ash.s' fi if test -f 'ash.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ash.f'\" else echo shar: Extracting \"'ash.f'\" \(4852 characters\) sed "s/^X//" >'ash.f' <<'END_OF_FILE' Xc April 8, 1986 Xc Xc Find bin counts of data array "x(n)" for ASH estimator Xc Xc "nbin" bins are formed over the interval [a,b) Xc Xc bin counts returned in array "nc" - # pts outside [a,b) = "nskip" X Xc ##### Copyright 1986 David W. Scott X X subroutine bin1 ( x , n , ab , nbin , nc , nskip ) X dimension x(n) , nc(nbin) , ab(2) X X nskip = 0 X a = ab(1) X b = ab(2) X X do 5 i = 1,nbin X nc(i) = 0 X5 continue X X d = (b-a) / nbin X X do 10 i = 1,n X k = (x(i)-a) / d + 1.0 X if (k.ge.1 .and. k.le.nbin) then X nc(k) = nc(k) + 1 X else X nskip = nskip + 1 X end if X10 continue X X return X end Xc April 8, 1986 Xc Xc Computer ASH density estimate; Quartic (biweight) kernel Xc Xc Average of "m" shifted histograms Xc Xc Bin counts in array "nc(nbin)" - from routine "nbin1" Xc Xc "nbin" bins are formed over the interval [a,b) Xc Xc ASH estimates returned in array "f(nbin)" Xc Xc FP-ASH plotted at a+d/2 ... b-d/2 where d = (b-a)/nbin Xc Xc Note: If "nskip" was nonzero, ASH estimates incorrect near boundary Xc Note: Should leave "m" empty bins on each end of array "nc" so f OK X Xc ##### Copyright 1986 David W. Scott X X subroutine ash1 ( m , nc , nbin , ab , t , f , w , ier ) X dimension nc(nbin) , ab(2) , t(nbin) , f(nbin) , w(m) X X ier = 0 X a = ab(1) X b = ab(2) X n = 0 X Xc-compute weights Xc w-array shifted by 1 X X mm1 = m-1 X xm = m X xm4 = xm**4 X cons = (15.0*xm4) / (16.0*xm4-1.0) X w(1) = cons X do 5 i = 1,mm1 X w(i+1) = cons * ( 1.0 - (i/xm)**2 )**2 X5 continue X Xc-check if estimate extends beyond mesh X X do 7 i = 1,mm1 X if( nc(i)+nc(nbin+1-i) .gt. 0) ier = 1 X7 continue X Xc-compute ash(m) estimate X X delta = (b-a) / nbin X h = m*delta X X do 10 i = 1,nbin X t(i) = a + (i-0.5)*delta X f(i) = 0.0 X n = n + nc(i) X10 continue X X do 20 i = 1,nbin X if (nc(i).eq.0) goto 20 X X c = nc(i) / (n*h) X do 15 k = max0(1,i-mm1) , min0(nbin,i+mm1) X f(k) = f(k) + c * w( iabs(k-i)+1 ) X15 continue X20 continue X X return X end Xc April 12, 1986 bin2.f Xc Xc Find bin counts of data array "x(n,2)" for ASH estimator Xc Xc "nbin(1)" by "nbin(2)" bins are formed Xc Xc x:axis [ ab(1,1) , ab(1,2) ) Xc y:axis [ ab(2,1) , ab(2,2) ) half-open Xc Xc bin counts returned in array "nc" - # pts outside [a,b) = "nskip" X Xc ##### Copyright 1986 David W. Scott X X subroutine bin2 ( x , n , ab , nbin , nc , nskip ) X dimension x(n,2) , nbin(2) , nc(nbin(1),nbin(2)) , ab(2,2) X X nskip = 0 X ax = ab(1,1) X bx = ab(1,2) X ay = ab(2,1) X by = ab(2,2) X X do 5 j = 1,nbin(2) X do 4 i = 1,nbin(1) X nc(i,j) = 0 X4 continue X5 continue X X dx = (bx-ax) / nbin(1) X dy = (by-ay) / nbin(2) X X do 10 i = 1,n X kx = (x(i,1)-ax) / dx + 1.0 X ky = (x(i,2)-ay) / dy + 1.0 X if (kx.ge.1 .and. kx.le.nbin(1) .and. X 1 ky.ge.1 .and. ky.le.nbin(2)) then X nc(kx,ky) = nc(kx,ky) + 1 X else X nskip = nskip + 1 X end if X10 continue X X return X end Xc April 12, 1986 ash2.f Xc Xc Computer ASH density estimate; Product Quartic (biweight) kernel Xc Xc Average of "m[1] by m[2]" shifted histograms Xc Xc Bin counts in matrix "nc" - from routine "nbin2" Xc Xc ASH estimates returned in matrix "f" Xc Xc FP-ASH plotted at a+d/2 ... b-d/2 where d = (b-a)/nbin Xc Xc Note: If "nskip" was nonzero, ASH estimates incorrect near boundary Xc Note: Should leave "m" empty bins on each end of array "nc" so f OK X Xc #### Copyright 1986 David W. Scott X X subroutine ash2 ( m , nc , nbinx , nbiny , ab , f , w , ier ) X dimension m(2) , nc(nbinx,nbiny) , ab(2,2) , f(nbinx,nbiny) X dimension w( m(1) , m(2) ) X X ier = 0 X ax = ab(1,1) X bx = ab(1,2) X ay = ab(2,1) X by = ab(2,2) X Xc-compute weights Xc w-array shifted by 1 X X mx = m(1) X my = m(2) X mxm1 = mx-1 X mym1 = my-1 X xm = mx X ym = my X xm4 = xm**4 X ym4 = ym**4 X consx = (15.0*xm4) / (16.0*xm4-1.0) X consy = (15.0*ym4) / (16.0*ym4-1.0) X Xc-----ASSUMES f dimensioned larger than m(1) or m(2) X Xc-put marginal weights in f array as work array X X f(1,1) = consx X f(2,1) = consy X do 5 i = 1,mxm1 X f(1,i+1) = consx * ( 1.0 - (i/xm)**2 )**2 X5 continue X do 6 i = 1,mym1 X f(2,i+1) = consy * ( 1.0 - (i/ym)**2 )**2 X X6 continue X Xc-computer product weight array (avoids later multipications) X X do 3 j = 1,my X do 2 i = 1,mx X w(i,j) = f(1,i) * f(2,j) X2 continue X3 continue X Xc-compute ash(m) estimate X X n = 0 X do 10 j = 1,nbiny X do 9 i = 1,nbinx X f(i,j) = 0.0 X n = n + nc(i,j) X9 continue X10 continue Xc-check if estimate extends beyond mesh X X ncheck = 0 X do 12 j = my , nbiny+1-my X do 11 i = mx , nbinx+1-mx X ncheck = ncheck + nc(i,j) X11 continue X12 continue X if (ncheck .ne. n) ier = 1 X X dx = (bx-ax) / nbinx X dy = (by-ay) / nbiny X X hx = mx*dx X hy = my*dy X X do 20 j = 1,nbiny X do 19 i = 1,nbinx X X if (nc(i,j).eq.0) goto 19 X X c = nc(i,j) / (n*hx*hy) X X do 18 ky = max0(1,j-mym1) , min0(nbiny,j+mym1) X do 17 kx = max0(1,i-mxm1) , min0(nbinx,i+mxm1) X f(kx,ky) = f(kx,ky) + c * X 1 w( iabs(kx-i)+1 , iabs(ky-j)+1 ) X17 continue X18 continue X X19 continue X20 continue X X return X end END_OF_FILE if test 4852 -ne `wc -c <'ash.f'`; then echo shar: \"'ash.f'\" unpacked with wrong size! fi # end of 'ash.f' fi if test -f 'penny' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'penny'\" else echo shar: Extracting \"'penny'\" \(900 characters\) sed "s/^X//" >'penny' <<'END_OF_FILE' X1964 55.0 X1961 56.6 X1956 54.0 X1982 54.4 X1985 55.6 X1987 55.0 X1985 56.0 X1970 57.2 X1954 52.2 X1983 52.8 X1971 58.0 X1955 52.6 X1949 53.2 X1979 54.0 X1986 55.0 X1956 53.6 X1972 57.0 X1976 53.2 X1970 55.6 X1955 54.2 X1961 57.0 X1957 53.6 X1960 57.0 X1952 53.4 X1965 56.6 X1977 54.0 X1988 54.0 X1947 53.0 X1959 56.8 X1980 53.6 X1964 56.0 X1979 53.0 X1971 55.8 X1984 55.2 X1962 54.6 X1953 50.6 X1951 53.4 X1989 56.2 X1978 53.0 X1983 52.6 X1948 53.6 X1965 58.2 X1974 59.0 X1949 51.6 X1966 57.2 X1963 58.0 X1980 53.0 X1967 57.6 X1989 54.0 X1974 54.0 X1951 54.0 X1975 54.0 X1969 58.0 X1969 56.0 X1988 55.6 X1946 53.2 X1958 53.2 X1968 57.2 X1967 55.8 X1963 56.2 X1966 56.0 X1950 53.4 X1975 54.2 X1977 53.2 X1950 51.0 X1982 54.6 X1986 53.2 X1957 55.0 X1953 54.8 X1981 54.8 X1945 51.8 X1948 52.2 X1947 53.0 X1946 54.2 X1984 54.0 X1972 57.0 X1976 53.2 X1978 53.6 X1952 54.6 X1987 54.0 X1968 56.0 X1958 53.2 X1960 55.2 X1973 59.0 X1959 57.2 X1962 56.0 X1981 54.0 X1954 54.6 X1945 54.6 X1973 57.6 END_OF_FILE if test 900 -ne `wc -c <'penny'`; then echo shar: \"'penny'\" unpacked with wrong size! fi # end of 'penny' fi if test -f 'penny.plot' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'penny.plot'\" else echo shar: Extracting \"'penny.plot'\" \(860 characters\) sed "s/^X//" >'penny.plot' <<'END_OF_FILE' X#~New session: Time: 647244866; Version: "S Sat Apr 7 17:28:40 CDT 1990" X Xpenny <- matrix(scan("penny"), 90, 2, byrow = T) Xpo <- order(penny[, 1]) Xpenny[po, ] X Xpar(mfrow=c(2,2)) X Xhist(penny[, 2],xlab="penny thickness (mils)") Xplot(penny[, 1], penny[, 2], pch = "+", X xlab="year",ylab="thickness") X Xdyn.load("ash.o") Xab <- c(48, 61) Xbx <- bin1(penny[, 2], ab, 260) Xfx <- ash1(bx, 28) Xplot(fx, type = "l",xlab="penny thickness (mils)",ylab="", X main="ash m=28") Xaddbump(ash1(bx, 28)) X Xab2 <- matrix(c(1941, 48, 1993, 62), 2, 2) Xbx2 <- bin2(penny, ab2, c(40, 40)) Xfx2 <- ash2(bx2) X Xcontour(fx2, v = fx2$v, labex = F,xlab="year",ylab="thickness") Xtitle("ash2 m=(5,5) & lowess f=.2") X Xpoints(penny[,1],penny[,2],pch="o") Xlines(lowess(penny[, 1], penny[, 2], f = 0.2), lty = 3) Xstamp() X# persp(fx2$z/max(fx2$z)) X# q() X#~End session: Time: 647246282; Process: 4268 END_OF_FILE if test 860 -ne `wc -c <'penny.plot'`; then echo shar: \"'penny.plot'\" unpacked with wrong size! fi # end of 'penny.plot' fi if test -f 'penny5.plot' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'penny5.plot'\" else echo shar: Extracting \"'penny5.plot'\" \(616 characters\) sed "s/^X//" >'penny5.plot' <<'END_OF_FILE' X# 6-3-91 --- 7-3-91 X Xprint("postscript?") Xif(scan(n=1)==1) postscript("1.ps",horiz=F) X Xpar(mfrow=c(3,2),pty="s") Xpar(cex=.8) Xpar(mar=c(6,6,1,1)) X Xprint("be sure to run penny.plot first----") X Xplot(penny[, 1], penny[, 2], pch = "o", X xlab="year",ylab="penny thickness (mils)", X xlim=c(1942,1992),ylim=c(49.5,60)) X Xlines(lowess(penny[, 1], penny[, 2], f = 0.2), lwd=2) Xtext(1970,51,"lowess f=0.20") X Xplot(penny[, 1], penny[, 2], type="n", X xlab="year",ylab="", X xlim=c(1942,1992),ylim=c(49.5,60)) X Xcontour(fx2,v=fx2$v,labex=F,lty=2,add=T) X Xmodalreg(fx2,nintx=1) X Xlines(lowess(penny[, 1], penny[, 2], f = 0.2), lwd=2) X END_OF_FILE if test 616 -ne `wc -c <'penny5.plot'`; then echo shar: \"'penny5.plot'\" unpacked with wrong size! fi # end of 'penny5.plot' fi echo shar: End of shell archive. exit 0