Examples for ddst package

Data Driven Smooth Tests with ddst package

This document presents tests avaliable in the ddst package.

library("ddst")
library("polynom")

II. Data Driven Smooth Tests for Selected Goodness-of-Fit Problems

1. Data Driven Smooth Test for Uniformity

set.seed(7)
# H0 is true
z <- runif(80)
t <- ddst.uniform.test(z, compute.p = TRUE, 
                       compute.cv = TRUE, 
                       d.n = 10)
t
## 
##  Data Driven Smooth Test for Uniformity
## 
## data:  z, base: ddst.base.legendre  c: 2.4  d.n: 10  cv(0.05) : 0.0045668)
## WT = 0.31054, T = 1, p-value = 0.5989
plot(t)

# H0 is false
z <- rbeta(80,4,2)
t <- ddst.uniform.test(z, compute.p = TRUE,compute.cv = TRUE,  d.n = 10)
t
## 
##  Data Driven Smooth Test for Uniformity
## 
## data:  z, base: ddst.base.legendre  c: 2.4  d.n: 10  cv(0.05) : 0.0044916)
## WT = 40.553, T = 3, p-value = 8e-05
t$p.value
## [1] 8e-05
plot(t)

2. Data Driven Smooth Test for Exponentiality

set.seed(7)
# H0 is true
z <- rexp(80,4)
t <- ddst.exp.test(z, compute.p = TRUE, d.n = 10)
t
## 
##  Data Driven Smooth Test for Expotentiality
## 
## data:  z, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 2.1593, T* = 1, p-value = 0.1748
plot(t)

# H0 is false
z = rchisq(80,4)
t = ddst.exp.test (z, compute.p = TRUE, d.n = 10)
t
## 
##  Data Driven Smooth Test for Expotentiality
## 
## data:  z, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 15.543, T* = 1, p-value = 0.00484
t$p.value
## [1] 0.00484
plot(t)

3. Data Driven Smooth Tests for Normality

3.1. Bounded Basis Functions

set.seed(7)
# H0 is true
z <- rnorm(100)
# let's look on first 10 coordinates
d.n <- 10
t <- ddst.normbounded.test(z, compute.p = TRUE, d.n = d.n)
t
## 
##  Data Driven Smooth Test for Normality - Bounded Basis Functions
## 
## data:  z, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 2.2232, T* = 1, p-value = 0.1581
plot(t)

# H0 is false
z <- rexp(100, 1)
t <- ddst.normbounded.test(z, compute.p = TRUE, d.n = d.n)
t
## 
##  Data Driven Smooth Test for Normality - Bounded Basis Functions
## 
## data:  z, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 299.45, T* = 10, p-value < 2.2e-16
plot(t)

# for Tephra data
z <- c(-1.748789, -1.75753, -1.740102, -1.740102, -1.731467, -1.765523,
     -1.761521, -1.72522, -1.80371, -1.745624, -1.872957, -1.729121,
     -1.81529, -1.888637, -1.887761, -1.881645, -1.91518, -1.849769,
     -1.755141, -1.665687, -1.764721, -1.736171, -1.736956, -1.737742,
     -1.687537, -1.804534, -1.790593, -1.808661, -1.784081, -1.729903,
     -1.711263, -1.748789, -1.772755, -1.72756, -1.71358, -1.821116,
     -1.839588, -1.839588, -1.830321, -1.807835, -1.747206, -1.788147,
     -1.759923, -1.786519, -1.726779, -1.738528, -1.754345, -1.781646,
     -1.641949, -1.755936, -1.775175, -1.736956, -1.705103, -1.743255,
     -1.82613, -1.826967, -1.780025, -1.684504, -1.751168)
t <- ddst.normbounded.test(z, compute.p = TRUE, Dmax = d.n)
t
## 
##  Data Driven Smooth Test for Normality - Bounded Basis Functions
## 
## data:  z, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 3.7615, T* = 1, p-value = 0.08891
plot(t)

3.2. Unbounded Basis Functions

set.seed(7)
# H0 is true
z <- rnorm(100)
# let's look on first 10 coordinates
d.n <- 10

# calculate finite sample corrections
# see 6.2. Composite null hypothesis H in the appendix materials
e.v <- ddst.normunbounded.bias(n = length(z), d.n = d.n)
e.v
## $e.0
##              [,1]       [,2]         [,3]       [,4]        [,5]      [,6]
## [1,] -0.003118967 -0.6935202 -0.002956126 -0.6265916 0.005252232 0.2493225
##              [,7]      [,8]        [,9]      [,10]
## [1,] -0.005572803 0.1643331 0.005080106 -0.3043065
## 
## $v.0
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,] 0.9478394 0.8147186 0.6556153 0.5867975 0.6271109 0.5654675 0.5199736
##           [,8]      [,9]    [,10]
## [1,] 0.5518196 0.4929518 0.489381
# simulated 1-alpha qunatiles, s(n, alpha) 
# see Table 1 in the JSCS article
s.n.alpha <- 4.4
r.alpha <- 2.708

t <- ddst.normubounded.test(z, d.n, e.v$e.0, e.v$v.0, r.alpha, s.n.alpha)
t
## 
##  Data Driven Smooth Test for Normality - Unbounded Basis Functions
## 
## data:  z, d.n: 10, r.alpha: 2.708, s.n.alpha: 4.4
## NAs = 2.7435, As = 1
plot(t)

# H0 is false, same lenght n = 100
z <- rexp(100, 1)

t <- ddst.normubounded.test(z, d.n, e.v$e.0, e.v$v.0, r.alpha, s.n.alpha)
t
## 
##  Data Driven Smooth Test for Normality - Unbounded Basis Functions
## 
## data:  z, d.n: 10, r.alpha: 2.708, s.n.alpha: 4.4
## NAs = 249.08, As = 9
plot(t)

# for Tephra data
z <- c(-1.748789, -1.75753, -1.740102, -1.740102, -1.731467, -1.765523,
       -1.761521, -1.72522, -1.80371, -1.745624, -1.872957, -1.729121,
       -1.81529, -1.888637, -1.887761, -1.881645, -1.91518, -1.849769,
       -1.755141, -1.665687, -1.764721, -1.736171, -1.736956, -1.737742,
       -1.687537, -1.804534, -1.790593, -1.808661, -1.784081, -1.729903,
       -1.711263, -1.748789, -1.772755, -1.72756, -1.71358, -1.821116,
       -1.839588, -1.839588, -1.830321, -1.807835, -1.747206, -1.788147,
       -1.759923, -1.786519, -1.726779, -1.738528, -1.754345, -1.781646,
       -1.641949, -1.755936, -1.775175, -1.736956, -1.705103, -1.743255,
       -1.82613, -1.826967, -1.780025, -1.684504, -1.751168)

# calculate finite sample corrections
e.v <- ddst.normunbounded.bias(n = length(z))
e.v
## $e.0
##              [,1]       [,2]          [,3]       [,4]        [,5]      [,6]
## [1,] -0.001745287 -0.7812591 -0.0003146211 -0.5289889 0.003045283 0.3671568
##              [,7]         [,8]       [,9]      [,10]        [,11]     [,12]
## [1,] -0.005166266 -0.006556786 0.00614837 -0.2153776 -0.005973216 0.2596683
##           [,13]      [,14]        [,15]      [,16]       [,17]   [,18]
## [1,] 0.00490113 -0.1825443 -0.003298743 0.05528399 0.001527895 0.06612
##            [,19]      [,20]
## [1,] 0.000115756 -0.1487876
## 
## $v.0
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,] 0.9422175 0.7691471 0.6028542 0.5825068 0.6072724 0.5109567 0.5263274
##          [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
## [1,] 0.524183 0.4622242 0.4965074 0.4581659 0.4443366 0.4631899 0.4147061
##          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]
## [1,] 0.4448983 0.4167827 0.4114816 0.4252005 0.3854767 0.4197615
# simulated 1-alpha qunatiles, s(n, alpha) and  s.o(n, alpha)
# see Table 1 in the JSCS article
s.n.alpha <- 3.3
r.alpha <- 2.142

t <- ddst.normubounded.test(z, d.n, e.v$e.0, e.v$v.0, r.alpha, s.n.alpha)
t
## 
##  Data Driven Smooth Test for Normality - Unbounded Basis Functions
## 
## data:  z, d.n: 10, r.alpha: 2.142, s.n.alpha: 3.3
## NAs = 66.704, As = 17
plot(t)

4. Data Driven Smooth Test for Extreme Value Distribution

library(evd)
set.seed(7)


# H0 is true
x <- -qgumbel(runif(100), -1, 1)
t <- ddst.evd.test (x, compute.p = TRUE, d.n = 10)
t
## 
##  Data Driven Smooth Test for Extreme Values
## 
## data:  x, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 1.4531, T* = 1, p-value = 0.9941
plot(t)

# H0 is false
x <- rexp(80,4)
t <- ddst.evd.test (x, compute.p = TRUE, d.n = 10)
t
## 
##  Data Driven Smooth Test for Extreme Values
## 
## data:  x, base: ddst.base.legendre, c: 100, d.n: 10
## W*T* = 11731652, T* = 10, p-value < 2.2e-16
plot(t)

III. Nonparametric Data Driven Smooth Tests for Comparing Distributions

5. Data Driven Smooth Test for k-Sample Problem

set.seed(7)
# H0 is false
x <- runif(80)
y <- rexp(80, 1)
t <- ddst.twosample.test(x, y, compute.p = TRUE)
t
## 
##  Data Driven Smooth Test for Two-Sample Problem
## 
## data:  listxy, c: 2, d.N: 12
## WT = 36.108, T = (6, 6)
plot(t)

# H0 is false
x <- runif(80)
y <- rexp(80, 1)
z <- runif(80)
t <- ddst.ksample.test(list(x, y, z))
t
## 
##  Data Driven Smooth Test for Two-Sample Problem
## 
## data:  listxyz, c: 2.3, d.N: 12
## WT = 73.211, T = (3, 12, 1)
plot(t)

# H0 is true
x <- rnorm(80)
y <- rnorm(80)
z <- rnorm(80)
t <- ddst.ksample.test(list(x, y, z))
t
## 
##  Data Driven Smooth Test for Two-Sample Problem
## 
## data:  listxyz, c: 2.3, d.N: 12
## WT = 0.063224, T = (1, 1, 1)
plot(t)

6. Data Driven Test for Stochastic Ordering

set.seed(7)
library("rmutil", warn.conflicts = FALSE)
# H0 is false
# 1. Pareto(1)/Pareto(1.5)
x <- rpareto(50, 2, 2)
y <- rpareto(50, 1.5, 1.5)
t <- ddst.forstochdom.test(x, y, t = 2.2, K.N = 4)
t
## 
##  Data Driven Stochastic Ordering Test
## 
## data:  x y, t: 2.2, K.N: 4
## QT = 124.46, T = 5
plot(t)

# H0 is false
# 2. Laplace(0,1)/Laplace(1,25)
x <- rlaplace(50, 0, 1)
y <- rlaplace(50, 1, 25)
t <- ddst.forstochdom.test(x, y, t = 2.2, K.N = 4)
t
## 
##  Data Driven Stochastic Ordering Test
## 
## data:  x y, t: 2.2, K.N: 4
## QT = 129.93, T = 5
plot(t)

# H0 is true
# 3. LN(0.85,0.6)/LN(1.2,0.2)
x <- rlnorm(50, 0.85, 0.6)
y <- rlnorm(50, 1.2, 0.2)
t <- ddst.forstochdom.test(x, y, t = 2.2, K.N = 4)
t
## 
##  Data Driven Stochastic Ordering Test
## 
## data:  x y, t: 2.2, K.N: 4
## QT = 0, T = 1
plot(t)

7. Two-Sample Test Against Stochastic Dominance

set.seed(7)
# H0 is true
x <- runif(80)
y <- runif(80)
t <- ddst.againststochdom.test(x, y, alpha = 0.05, t = 2.2, k.N = 4)
t
## 
##  Data Driven Test Against Stochastic Dominance
## 
## data:  x, alpha: 0.05, t: 2.2, k.N: 4
## VT = 0.9, T = 1
plot(t)

# H0 is false
x <- runif(80)
y <- rbeta(80, 4, 2)
t <- ddst.againststochdom.test(x, y, alpha = 0.05, t = 2.2, k.N = 4)
t
## 
##  Data Driven Test Against Stochastic Dominance
## 
## data:  x, alpha: 0.05, t: 2.2, k.N: 4
## VT = 41.6, T = 5
plot(t)

8. Data Driven Smooth Test for Upward Trend Alternative in k samples

set.seed(7)
# H0 is true
x = runif(80)
y = runif(80) + 0.2
z = runif(80) + 0.4
t <- ddst.upwardtrend.test(list(x, y, z), t.p = 2.2, t.n = 2.2)
t
## 
##  Data Driven k-Sample Upward Trend Test
## 
## data:  listxyz, t.p: 2.2, t.n: 2.2
## U.T = 1105, p = 3
plot(t)

# H0 is false
x1 = rnorm(80)
x2 = rnorm(80) + 2
x3 = rnorm(80) + 4
x4 = rnorm(80) + 3
t <- ddst.upwardtrend.test(list(x1, x2, x3, x4), t.p = 2.2, t.n = 2.2)
t
## 
##  Data Driven k-Sample Upward Trend Test
## 
## data:  listx1x2x3x4, t.p: 2.2, t.n: 2.2
## U.T = -3401.6, p = 4
plot(t)

9. Data Driven Smooth Test for Umbrella Alternatives in k samples

set.seed(7)
# H0 is true
x = runif(80)
y = runif(80) + 0.2
z = runif(80)
t <- ddst.umbrellaknownp.test(list(x, y, z), p = 2, t.p = 2.2, t.n = 2.2)
t
## 
##  Data Driven k-Sample Umbrella Test
## 
## data:  listxyz, t.p: 2.2, t.n: 2.2, p: 2
## U.T = 412.82, p = 2
plot(t)

# H0 is true
x1 = rnorm(80)
x2 = rnorm(80) + 2
x3 = rnorm(80) + 4
x4 = rnorm(80) + 3
x5 = rnorm(80) + 2
x6 = rnorm(80) + 1
x7 = rnorm(80)
t <- ddst.umbrellaknownp.test(list(x1, x2, x3, x4, x5, x6, x7), p = 3, t.p = 2.2, t.n = 2.2)
t
## 
##  Data Driven k-Sample Umbrella Test
## 
## data:  listx1x2x3x4x5x6x7, t.p: 2.2, t.n: 2.2, p: 3
## U.T = 4662.7, p = 3
plot(t)

t <- ddst.umbrellaknownp.test(list(x1, x2, x3, x4, x5, x6, x7), p = 5, t.p = 2.2, t.n = 2.2)
t
## 
##  Data Driven k-Sample Umbrella Test
## 
## data:  listx1x2x3x4x5x6x7, t.p: 2.2, t.n: 2.2, p: 5
## U.T = -3101, p = 5
plot(t)

Unknown peak

# true umbrella pattern
x1 = rnorm(80)
x2 = rnorm(80) + 2
x3 = rnorm(80) + 4
x4 = rnorm(80) + 3
x5 = rnorm(80) + 2
x6 = rnorm(80) + 1
x7 = rnorm(80)

##  peak is unknown, so we test for each possible position
statistics <- sapply(1:7, function(i)
ddst.umbrellaknownp.test(list(x1, x2, x3, x4, x5, x6, x7), p = i, t.p = 2.2, t.n = 2.2)$statistic)

statistics
##       U.T       U.T       U.T       U.T       U.T       U.T       U.T 
## -3887.927 -4291.024  4627.885 -3936.036 -3162.127 -2773.183 -2589.604
max(statistics)
## [1] 4627.885