#!/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 1997-09-18 14:19 CDT by . # Source directory was `/home/scottdw/research/sagae/DIST'. # # Existing files will *not* be overwritten unless `-c' is specified. # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 4738 -rw-r--r-- h.c2.m1.s2.est.q # 1572 -rw-r--r-- h.est.q # 141 -rw-r--r-- nicerange.q # 602 -rw-r--r-- readme # 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 _sh05760; then $echo 'x -' 'creating lock directory' else $echo 'failed to create lock directory' exit 1 fi # ============= h.c2.m1.s2.est.q ============== if test -f 'h.c2.m1.s2.est.q' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'h.c2.m1.s2.est.q' '(file already exists)' else $echo 'x -' extracting 'h.c2.m1.s2.est.q' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'h.c2.m1.s2.est.q' && # 10-31-96 X # count/mean preserving (h^2 bias) # f'' continuous # minimize R(f'') X # mk0 = xk0/n etc X # mathematica commands X # fkx = ak + bk(x-mk) + ck/2(x-mk)^2 + dk/6(x-mk)^3 + ek/24(x-mk)^4 + fk/120(x-mk)^5 # fk1x = ak1 + bk1(x-mk1) + ck1/2(x-mk1)^2 + dk1/6(x-mk1)^3 + ek1/24(x-mk1)^4 + fk1/120(x-mk1)^5 X # Ak = Flatten[ Solve[ Expand[ Integrate[ fkx, {x,mk-hk/2,mk+hk/2} ] ] == mk0,ak] ] (* ak *) # Bk = Flatten[ Solve[ Expand[ Integrate[ (x-mk)fkx, {x,mk-hk/2,mk+hk/2} ] ] == mk1 ,bk] ] (* bk *) X # Ak1 = Flatten[ Solve[ Expand[ Integrate[ fk1x, {x,mk1-hk1/2,mk1+hk1/2} ] ] == mk10,ak1] ] (* ak1 *) # Bk1 = Flatten[ Solve[ Expand[ Integrate[ (x-mk1)fk1x, {x,mk1-hk1/2,mk1+hk1/2} ] ] == mk11 ,bk1] ] (* bk1 *) X # C0 = ( fkx /. x -> mk + hk/2 ) - ( fk1x /. x -> mk1 - hk1/2 ) (* continuity f ==0 *) # C1 = ( D[fkx,x] /. x -> mk + hk/2 ) - ( D[fk1x,x] /. x -> mk1 - hk1/2 ) (* continuity f'==0 *) # C2 = ( D[fkx,x,x] /. x -> mk + hk/2 ) - ( D[fk1x,x,x] /. x -> mk1 - hk1/2 ) (* continuity f''==0 *) # # Expand[ ( ( ( C0 /. Ak ) /. Ak1 ) /. Bk ) /. Bk1 ] (* C0 condition with substitutions for ak,ak1,bk,bk1 *) # Expand[ ( ( ( C1 /. Ak ) /. Ak1 ) /. Bk ) /. Bk1 ] (* C1 condition with substitutions for ak,ak1,bk,bk1 *) # Expand[ ( ( ( C2 /. Ak ) /. Ak1 ) /. Bk ) /. Bk1 ] (* C2 condition with substitutions for ak,ak1,bk,bk1 *) # # Expand[ Integrate [ (D[fkx,x,x])^2, {x,mk-hk/2,mk+hk/2} ] ] (* R(f'') criterion *) # # # # X X X # Oct 31, 1996 X # New constrained histogram estimator Sagae/Scott # Bin-Centered polynomial implementation # 0,1-order moment constraints; f,f'f'' continuous # ssmoothing minimizer # All constraints exactly on each bin interval X cbh.c2.m1.s2 <- function( x, m=10, tk ) { # m bins m+1 mesh nodes in tk X X if(missing(tk)){ab<-nicerange(x); tk <- seq(ab[1],ab[2],length=m+1) } X else{ m <- length(tk)-1 } X n <- length(x) X kbin <- integer(n) X for(k in seq(m)) { kbin[ seq(n)[ x>=tk[k] & x0) { xk <- xk - mk[k] } # centered on bin center X if(s0>0) { s0 <- s0/n; s1 <- sum(xk)/n } X mk0s[k] <- s0; mk1s[k] <- s1 X } X X # y <- c( c1...f1,c2...f2,...,cm...fm) X A <- matrix(0,4*m,4*m) # quadratic criterion X B <- matrix(0,3*(m-1),4*m) # continuity constraint matrix X C <- rep(0,3*(m-1)) # right hand side constraints X for(k in seq(1,length=m-1)) { # set up equations X io <- 0 # row offset: f continuous conditions X jo <- 4*(k-1) # column offset: variables c,d,e,f X B[io+k,jo+1] <- hk[k]^2/12; B[io+k,jo+5] <- -hk[k+1]^2/12 X B[io+k,jo+2] <- hk[k]^3/120; B[io+k,jo+6] <- hk[k+1]^3/120 X B[io+k,jo+3] <- hk[k]^4/480; B[io+k,jo+7] <- -hk[k+1]^4/480 X B[io+k,jo+4] <- hk[k]^5/6720; B[io+k,jo+8] <- hk[k+1]^5/6720 X C[io+k] <- -mk0s[k]/hk[k] + mk0s[k+1]/hk[k+1] - 6*mk1s[k]/hk[k]^2 - 6*mk1s[k+1]/hk[k+1]^2 X X io <- m-1 # row offset: f' continuous conditions X B[io+k,jo+1] <- hk[k]/2; B[io+k,jo+5] <- hk[k+1]/2 X B[io+k,jo+2] <- hk[k]^2/10; B[io+k,jo+6] <- -hk[k+1]^2/10 X B[io+k,jo+3] <- hk[k]^3/48; B[io+k,jo+7] <- hk[k+1]^3/48 X B[io+k,jo+4] <- hk[k]^4/420; B[io+k,jo+8] <- -hk[k+1]^4/420 X C[io+k] <- -12*mk1s[k]/hk[k]^3 + 12*mk1s[k+1]/hk[k+1]^3 X X io <- 2*(m-1) # row offset: f'' continuous conditions X B[io+k,jo+1] <- 1; B[io+k,jo+5] <- -1 X B[io+k,jo+2] <- hk[k]/2; B[io+k,jo+6] <- hk[k+1]/2 X B[io+k,jo+3] <- hk[k]^2/8; B[io+k,jo+7] <- -hk[k+1]^2/8 X B[io+k,jo+4] <- hk[k]^3/48; B[io+k,jo+8] <- hk[k+1]^3/48 X } X X for(k in seq(1,length=m)) { # set up equations X jo <- 4*(k-1) # column offset: variables c,d,e,f X X A[jo+1,jo+1] <- hk[k] # ck^2 X A[jo+2,jo+2] <- hk[k]^3/12 # dk^2 X A[jo+1,jo+3] <- A[jo+3,jo+1] <- hk[k]^3/12/2 # ck ek (split ck*ek, ek*ck) X A[jo+3,jo+3] <- hk[k]^5/320 # ek^2 X A[jo+2,jo+4] <- A[jo+4,jo+2] <- hk[k]^5/240/2 # dk*fk fk*dk X A[jo+4,jo+4] <- hk[k]^7/16128 # fk^2 X } # browser() X y <- solve(A) %*% t(B) %*% solve( B %*% solve(A) %*% t(B) ) %*% C X i <- seq(0,by=4,length=m) X ck <- y[1+i]; dk <- y[2+i]; ek <- y[3+i]; fk <- y[4+i] X ak <- -ck*hk^2/24 - ek*hk^4/1920 + mk0s/hk X bk <- -dk*hk^2/40 - fk*hk^4/4480 + 12*mk1s/hk^3 X X for(k in seq(m)) { X tt <- seq(tk[k],tk[k+1],l=251) X tm <- tt - mk[k] # centered X yy <- ak[k] + bk[k]*tm + ck[k]*tm^2/2 + dk[k]*tm^3/6 + ek[k]*tm^4/24 + fk[k]*tm^5/120 X lines( tt , yy ) # browser() X } # browser() # invisible(list(tk=tk,n=n,i=i,mk0s=mk0s,mk1s=mk1s,mks=mks,h=h)) } SHAR_EOF $shar_touch -am 0918120097 'h.c2.m1.s2.est.q' && chmod 0644 'h.c2.m1.s2.est.q' || $echo 'restore of' 'h.c2.m1.s2.est.q' '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 'h.c2.m1.s2.est.q:' 'MD5 check failed' 27835b6322a7b1aabc7c3b33cd7b1fa9 h.c2.m1.s2.est.q SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'h.c2.m1.s2.est.q'`" test 4738 -eq "$shar_count" || $echo 'h.c2.m1.s2.est.q:' 'original size' '4738,' 'current size' "$shar_count!" fi fi # ============= h.est.q ============== if test -f 'h.est.q' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'h.est.q' '(file already exists)' else $echo 'x -' extracting 'h.est.q' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'h.est.q' && # Oct 12, 1996 X # New constrained histogram estimator Sagae/Scott # Orthogonal polynomial implementation # 0 <= p <= 5 # All constraints exactly on each bin interval X cbh <- function( x, h, p=0, ab ) { X if(h<=0|p<0|p>5) stop("h must be positive, and p must be 0-5") X if(missing(ab)){ab<-nicerange(x)}; tk <- seq(ab[1],ab[2]+.99*h,by=h) X n <- length(x) X kbin <- floor( 1 + (x-tk[1])/h ) X i <- max(table(kbin)); plot(range(tk),c(0,1.2*i/(n*h)),type="n") # plot space X X for(k in seq(length(tk)-1)) { X mk <- tk[k]+h/2 # bin midpoint X xk <- ( x[kbin==k] - mk ) / h # scaled/centered (-.5,.5) X s0 <- length(xk) X s1 <- s2 <- s3 <- s4 <- s5 <- 0 X if(s0>0) { s0 <- s0/n; s1 <- sum(xk)/n; s2 <- sum(xk^2)/n; s3 <- sum(xk^3)/n; s4 <- sum(xk^4)/n; s5 <- sum(xk^5)/n } X X a0 <- s0; a1 <- a2 <- a3 <- a4 <- a5 <- 0 # orthogonal polynomial estimated weights X if(p>=1) { a1 <- 2*sqrt(3)*s1 X if(p>=2) { a2 <- sqrt(5) * (s0/2 - 6 *s2) X if(p>=3) { a3 <- sqrt(7) * (-3*s1 + 20*s3) X if(p>=4) { a4 <- 9*s0/8 - 45*s2 + 210*s4 X if(p>=5) { a5 <- sqrt(11) * (15*s1/4 - 70*s3 + 252*s5) X } } } } } X X tt <- c(-.5,.5); if(p>1){ tt <- seq(tt[1],tt[2],length=5*p) } X y <- rep(a0,length(tt)) X if(p>=1) { y <- y + a1 * 2*sqrt(3)*tt X if(p>=2) { y <- y + a2 * sqrt(5) * (0.5 - 6*tt^2) X if(p>=3) { y <- y + a3 * sqrt(7) * (3*tt - 20*tt^3) X if(p>=4) { y <- y + a4 * (9/8 - 45*tt^2 + 210*tt^4) X if(p>=5) { y <- y + a5 * sqrt(11) * (15*tt - 70*tt^3 + 252*tt^5) X } } } } } # browser() ##??? lines( mk + h*tt , h*y ) X lines( mk + h*tt , y/h ) X } } SHAR_EOF $shar_touch -am 0918115797 'h.est.q' && chmod 0644 'h.est.q' || $echo 'restore of' 'h.est.q' '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 'h.est.q:' 'MD5 check failed' fc735782427a5519bfb49d2f46a71429 h.est.q SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'h.est.q'`" test 1572 -eq "$shar_count" || $echo 'h.est.q:' 'original size' '1572,' 'current size' "$shar_count!" fi fi # ============= nicerange.q ============== if test -f 'nicerange.q' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'nicerange.q' '(file already exists)' else $echo 'x -' extracting 'nicerange.q' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'nicerange.q' && nicerange <- function(x,beta=.1) { ab <- range(x) # interval surrounding data X del <- (ab[2]-ab[1])*beta/2; return( c( ab+c(-del,del) ) ) } X SHAR_EOF $shar_touch -am 0918141897 'nicerange.q' && chmod 0644 'nicerange.q' || $echo 'restore of' 'nicerange.q' '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 'nicerange.q:' 'MD5 check failed' 5bbad212f525f1019735389e38867363 nicerange.q SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nicerange.q'`" test 141 -eq "$shar_count" || $echo 'nicerange.q:' 'original size' '141,' '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' && 9-17-97 X Couple of S programs for the Scott/Sagae and Sagae/Scott polynomial histogram programs. X cbh == function (once called constrained binned histogram) X X x == data vector X h == bin width X p == order of polynomial histogram (0,1,2,3,4,5 only) X ab== optional interval of support X X cbh.c2.m1.s2 == smoothing spline version of cbh X X c2 means smoothness like cubic spline X m1 means bin count and bin means matched X s2 means 2 orders of roughness smoothing applied X X x == data vector X m == number of bins (only one of m or tk should be supplied) X tk== vector of mesh points (need not be equally spaced) SHAR_EOF $shar_touch -am 0918120397 'readme' && chmod 0644 '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' 0f26a67f8c196c1ccad441a8926dbeb0 readme SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'readme'`" test 602 -eq "$shar_count" || $echo 'readme:' 'original size' '602,' 'current size' "$shar_count!" fi fi rm -fr _sh05760 exit 0