R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("sac-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('sac')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "BootsChapt"
> 
> ### * BootsChapt
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: BootsChapt
> ### Title: Bootstrap (Permutation) Test of Change-Point(s) with One-Change
> ###   or Epidemic Alternative
> ### Aliases: BootsChapt
> ### Keywords: htest
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> 
> # one-change alternative
> k<-10
> n<-20
> x<-rnorm(n,0,1)
> x[(k+1):n]<-x[(k+1):n]+1.5
> T<-SemiparChangePoint(x, alternative = "one.change")$Sn
> BootsChapt(x, T, B = 5)
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
bootstrap p-value of Tn 
                    0.2 
>     #Choose larger B to get better approximate p-value.
> 
> 
> 
> cleanEx(); ..nameEx <- "BootsModelTest"
> 
> ### * BootsModelTest
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: BootsModelTest
> ### Title: Bootstrap Test of the Validity of the Semiparametric
> ###   Change-Point Model
> ### Aliases: BootsModelTest
> ### Keywords: htest
> 
> ### ** Examples
> 
> ## Nile data with one change-point: the annual flows drop in 1898.
> ## It is believed to be caused by the building of the first Aswan dam.
> if(! "package:stats" %in% search()) library(stats)
> data(Nile)
> require(sac) #load the package
[1] TRUE
> Nile.res<-SemiparChangePoint(Nile, alternative = "one.change")
> BootsModelTest(Nile, Nile.res$k.hat, length(Nile), B=5, Nile.res$alpha.hat, 
+     Nile.res$beta.hat)
$Delta
[1] 0.3138604

$Pvalue
[1] 0.8

>     # Choose larger B to get better approximate p-value.
>     # It takes longer to do bootstrap model test for large B.
> 
> 
> 
> cleanEx(); ..nameEx <- "CriticalValues"
> 
> ### * CriticalValues
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: CriticalValues
> ### Title: Critical Values of Tests of Change-Point(s) with One-Change or
> ###   Epidemic Alternative
> ### Aliases: Sn.alfa CV.Epidemic.Vn CV.Epidemic.Wn
> ### Keywords: htest
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> alpha<-0.05
> n<-20
> d<-1
> Sn.alfa(alpha, n, d, model="semiparametric")
[1] 9.700366
> CV.Epidemic.Vn(alpha, d)
[1] 0.3737601
> CV.Epidemic.Wn(alpha)
[1] 3.052917
> 
> 
> 
> cleanEx(); ..nameEx <- "Pvalues"
> 
> ### * Pvalues
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Pvalues
> ### Title: The p-values of Test Statistics Based on Asymptotic Distribution
> ### Aliases: p.OneChange p.Epidemic.Vn p.Epidemic.Wn
> ### Keywords: htest
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> # one-change alternative
> k<-10
> n<-30
> x<-rnorm(n,0,1)
> x[(k+1):n]<-x[(k+1):n]+1.5
> T<-SemiparChangePoint(x, alternative = "one.change")$Sn
> p.OneChange(n, d=1, T)
[1] 0.02249479
> 
> # epidemic alternative
> k<-5
> m<-10
> n<-20
> x<-rnorm(n,0,1)
> x[(k+1):m]<-x[(k+1):m]+1.5
> res<-SemiparChangePoint(x, alternative = "e")
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
> V<-res$Vn; W<-res$Wn
> p.Epidemic.Vn(V, d=1)
[1] 0.1103208
> p.Epidemic.Wn(W)
[1] 0.04370378
> 
> 
> 
> cleanEx(); ..nameEx <- "SemiparChangePoint"
> 
> ### * SemiparChangePoint
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: SemiparChangePoint
> ### Title: Semiparametric Test of Change-point(s) with One-change or
> ###   Epidemic Alternative
> ### Aliases: SemiparChangePoint
> ### Keywords: htest models
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> # one-change alternative
> k<-10
> n<-30
> x<-rnorm(n,0,1)
> x[(k+1):n]<-x[(k+1):n]+1.5
> SemiparChangePoint(x, alternative = "one.change")
$k.hat
[1] 10

$ll
 [1] -102.03592 -100.55701 -100.50180  -98.58183 -100.15195  -99.82235
 [7]  -98.11321  -97.88209  -97.88650  -97.62050  -96.11530  -98.62219
[13]  -99.43404  -99.32114  -97.61652  -99.13115  -99.41321  -99.68080
[19] -100.56188 -101.17521 -101.55127 -101.89401 -102.02980 -102.03330
[25] -101.87696 -102.02469 -102.03589 -102.02474 -101.90455 -101.71220
[31] -102.03592

$Sn
[1] 11.84124

$alpha.hat
(Intercept) 
  -1.222546 

$beta.hat
       x 
1.444638 

> 
> # epidemic alternative
> k<-5
> m<-10
> n<-20
> x<-rnorm(n,0,1)
> x[(k+1):m]<-x[(k+1):m]+1.5
> SemiparChangePoint(x,  alternative = "epidemic")
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
Warning in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
	 fitted probabilities numerically 0 or 1 occurred
$k.hat
[1] 7

$m.hat
[1] 10

$ll
           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
 [1,] -59.91465 -59.69005 -59.75287 -59.55097 -57.92767 -58.87672 -59.39492
 [2,] -59.69005 -59.91465 -59.90710 -59.77465 -58.28901 -59.21250 -59.63870
 [3,] -59.75287 -59.90710 -59.91465 -59.72648 -57.56715 -59.11478 -59.64176
 [4,] -59.55097 -59.77465 -59.72648 -59.91465 -55.94434 -59.36118 -59.81055
 [5,] -57.92767 -58.28901 -57.56715 -55.94434 -59.91465 -59.72899 -59.50395
 [6,] -58.87672 -59.21250 -59.11478 -59.36118 -59.72899 -59.91465 -59.71564
 [7,] -59.39492 -59.63870 -59.64176 -59.81055 -59.50395 -59.71564 -59.91465
 [8,] -59.78718 -59.89124 -59.90065 -59.90463 -58.96414 -59.21703 -59.43264
 [9,] -59.78710 -59.58285 -59.50487 -59.12045 -56.57428 -56.86137 -56.88768
[10,] -59.01380 -58.52410 -58.34722 -57.58415 -52.45139 -52.99176 -51.46046
[11,] -59.40056 -59.05427 -58.95352 -58.44996 -55.60251 -56.24271 -56.74592
[12,] -59.69763 -59.47122 -59.41493 -59.08672 -57.18474 -57.74002 -58.21074
[13,] -59.62234 -59.37075 -59.31996 -58.98160 -57.05876 -57.68161 -58.19816
[14,] -59.58288 -59.32128 -59.28122 -58.94859 -57.05462 -57.72896 -58.27431
[15,] -59.88720 -59.77826 -59.75624 -59.58010 -58.40163 -58.86961 -59.23748
[16,] -59.83477 -59.91202 -59.91455 -59.88331 -59.26195 -59.55025 -59.75167
[17,] -59.78387 -59.90335 -59.91111 -59.89660 -59.31500 -59.59972 -59.78756
[18,] -59.83450 -59.91461 -59.91321 -59.85902 -59.11074 -59.47456 -59.71425
[19,] -59.46910 -59.83202 -59.87039 -59.91440 -59.42138 -59.69880 -59.85375
[20,] -59.51397 -59.89600 -59.91074 -59.88471 -59.09847 -59.52332 -59.76588
           [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
 [1,] -59.78718 -59.78710 -59.01380 -59.40056 -59.69763 -59.62234 -59.58288
 [2,] -59.89124 -59.58285 -58.52410 -59.05427 -59.47122 -59.37075 -59.32128
 [3,] -59.90065 -59.50487 -58.34722 -58.95352 -59.41493 -59.31996 -59.28122
 [4,] -59.90463 -59.12045 -57.58415 -58.44996 -59.08672 -58.98160 -58.94859
 [5,] -58.96414 -56.57428 -52.45139 -55.60251 -57.18474 -57.05876 -57.05462
 [6,] -59.21703 -56.86137 -52.99176 -56.24271 -57.74002 -57.68161 -57.72896
 [7,] -59.43264 -56.88768 -51.46046 -56.74592 -58.21074 -58.19816 -58.27431
 [8,] -59.91465 -55.94434 -53.41299 -57.65592 -58.90867 -58.90600 -58.97245
 [9,] -55.94434 -59.91465 -58.02774 -59.56538 -59.89073 -59.86733 -59.86552
[10,] -53.41299 -58.02774 -59.91465 -59.63892 -59.23908 -59.57882 -59.67467
[11,] -57.65592 -59.56538 -59.63892 -59.91465 -59.55465 -59.81746 -59.85947
[12,] -58.90867 -59.89073 -59.23908 -59.55465 -59.91465 -59.89356 -59.89633
[13,] -58.90600 -59.86733 -59.57882 -59.81746 -59.89356 -59.91465 -59.91299
[14,] -58.97245 -59.86552 -59.67467 -59.85947 -59.89633 -59.91299 -59.91465
[15,] -59.65032 -59.86747 -59.04809 -59.39979 -59.72012 -59.49626 -58.92160
[16,] -59.90566 -59.51302 -58.05138 -58.56007 -59.06482 -58.56987 -57.65257
[17,] -59.91235 -59.47492 -58.08504 -58.60747 -59.10543 -58.75490 -58.30590
[18,] -59.89641 -59.60001 -58.46542 -58.93844 -59.36478 -59.15805 -58.95919
[19,] -59.91015 -59.32920 -57.92441 -58.51256 -59.04024 -58.79686 -58.58098
[20,] -59.91120 -59.50695 -58.34154 -58.86932 -59.32005 -59.16355 -59.05228
          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]
 [1,] -59.88720 -59.83477 -59.78387 -59.83450 -59.46910 -59.51397
 [2,] -59.77826 -59.91202 -59.90335 -59.91461 -59.83202 -59.89600
 [3,] -59.75624 -59.91455 -59.91111 -59.91321 -59.87039 -59.91074
 [4,] -59.58010 -59.88331 -59.89660 -59.85902 -59.91440 -59.88471
 [5,] -58.40163 -59.26195 -59.31500 -59.11074 -59.42138 -59.09847
 [6,] -58.86961 -59.55025 -59.59972 -59.47456 -59.69880 -59.52332
 [7,] -59.23748 -59.75167 -59.78756 -59.71425 -59.85375 -59.76588
 [8,] -59.65032 -59.90566 -59.91235 -59.89641 -59.91015 -59.91120
 [9,] -59.86747 -59.51302 -59.47492 -59.60001 -59.32920 -59.50695
[10,] -59.04809 -58.05138 -58.08504 -58.46542 -57.92441 -58.34154
[11,] -59.39979 -58.56007 -58.60747 -58.93844 -58.51256 -58.86932
[12,] -59.72012 -59.06482 -59.10543 -59.36478 -59.04024 -59.32005
[13,] -59.49626 -58.56987 -58.75490 -59.15805 -58.79686 -59.16355
[14,] -58.92160 -57.65257 -58.30590 -58.95919 -58.58098 -59.05228
[15,] -59.91465 -58.88425 -59.31357 -59.65643 -59.40091 -59.65774
[16,] -58.88425 -59.91465 -59.90373 -59.90976 -59.86388 -59.91269
[17,] -59.31357 -59.90373 -59.91465 -59.87538 -59.87765 -59.91446
[18,] -59.65643 -59.90976 -59.87538 -59.91465 -59.68252 -59.89779
[19,] -59.40091 -59.86388 -59.87765 -59.68252 -59.91465 -59.83563
[20,] -59.65774 -59.91269 -59.91446 -59.89779 -59.83563 -59.91465

$Vn
[1] 0.2935622

$Wn
[1] 3.134569

$alpha.hat
(Intercept) 
  -719.6614 

$beta.hat
       x 
515.4679 

> 
> 
> 
> cleanEx(); ..nameEx <- "cumsum.test"
> 
> ### * cumsum.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cumsum.test
> ### Title: Nonparametric Test for Change-Point with One-change or Epidemic
> ###   Alternative
> ### Aliases: cumsum.test
> ### Keywords: nonparametric
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> # one-change alternative
> k<-10
> n<-30
> x<-rnorm(n,0,1)
> x[(k+1):n]<-x[(k+1):n]+1.5
> cumsum.test(x, alternative = "one.change")
$Sn
         [,1]
[1,] 2.660218

$k.hat
[1] 10

> # epidemic alternative
> k<-10
> m<-20
> n<-30
> x<-rnorm(n,0,1)
> x[(k+1):m]<-x[(k+1):m]+1.5
> cumsum.test(x,  alternative = "epidemic")
$Sn
          [,1]
[1,] 0.3078001

$k.hat
[1] 8

$m.hat
[1] 20

> 
> 
> 
> 
> cleanEx(); ..nameEx <- "plots"
> 
> ### * plots
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plots
> ### Title: Visualized Model Diagnostic and Loglikelihood Plot
> ### Aliases: Graf.Diagnostic Plot.ll
> ### Keywords: device
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> k<-30
> n<-80
> x<-rnorm(n,0,1)
> x[(k+1):n]<-x[(k+1):n]+1.5
> res<-SemiparChangePoint(x, alternative = "one.change")
> Plot.ll(x, res$ll, col="blue")
change-point =  30 
> 
> ## Nile data with one change-point: the annual flows drop in 1898 which corresponds 
> ## to k=28. It is believed to be caused by the building of the first Aswan dam.
> if(! "package:sac" %in% search()) library(sac) 
>     #if package sac has not been loaded, load it.
> if(! "package:stats" %in% search()) library(stats)
> data(Nile)
> plot(Nile, type="p")
> Nile.res<-SemiparChangePoint(Nile, alternative = "one.change")
> Color<-c(1,2,3,4); LTY<-c(1,2,3,4)
> 
> ## Plots of estimated distribution functions
> Graf.Diagnostic(Nile, Nile.res$k.hat, length(Nile), Nile.res$alpha.hat, 
+     Nile.res$beta.hat, Color, LTY, xlab = "x", ylab = "Estimated DF's", 
+     main="Model Diagnostic for Nile Data", OneLegend = FALSE, lgnd1 = 
+     c(1100, 0.15), lgnd2 = c(600, .99), arw1=c(780, .93, 1010, .9),
+     arw2 = c(1165, .15, 1015, .24))
> 
> ## Plot of loglikelihood function
> Plot.ll(Nile, Nile.res$ll, col = "blue")
change-point =  28 
> Plot.ll(Nile, Nile.res$ll, col = "blue", xaxis.lab = seq(1871,1970, length = 100), 
+     xlab = "Year")
change-point =  28 
> 
> 
> 
> cleanEx(); ..nameEx <- "schapt"
> 
> ### * schapt
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: schapt
> ### Title: Semiparametric Analysis of Changepoint
> ### Aliases: schapt
> ### Keywords: htest robust
> 
> ### ** Examples
> 
> require(sac) #load the package
[1] TRUE
> # one-change alternative
> ## Nile data with one change-point: the annual flows drop in 1898.
> ## It is believed to be caused by the building of the first Aswan dam.
> if(! "package:sac" %in% search()) library(sac) 
>     #if package sac has not been loaded, load it.
> if(! "package:stats" %in% search()) library(stats)
> data(Nile)
> plot(Nile, type="p")
> schapt(Nile, alternative = "one.change")
$data.name
[1] "Nile"

$parameter
sample size          df 
        100           1 

$alternative
[1] "one.change"

$statistic
      Sn 
51.96226 

$change.point
 k 
28 

$estimate
      alpha        beta 
13.12708268 -0.01350496 

$p.value
[1] 0

> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  42.52 0.18 42.97 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')