ktp.Rout.save
R : Copyright 2006, The R Foundation for Statistical Computing
Version 2.3.1 (2006-06-01)
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 an HTML browser interface to help.
Type 'q()' to quit R.
>
> library(aster)
Loading required package: trust
>
> do.chisq.test <- function(x, k, mu, max.bin) {
+ stopifnot(all(x > k))
+ stopifnot(k + 1 < max.bin)
+ xx <- seq(k + 1, max.bin)
+ yy <- dpois(xx, mu)
+ yy[length(yy)] <- ppois(max.bin - 1, mu, lower.tail = FALSE)
+ pp <- yy / sum(yy)
+ ecc <- length(x) * pp
+ if (any(ecc < 5.0))
+ warning("violates rule of thumb about > 5 expected in each cell")
+ cc <- tabulate(x, max.bin)
+ cc <- cc[xx]
+ cc[length(cc)] <- nsim - sum(cc[- length(cc)])
+ chisqstat <- sum((cc - ecc)^2 / ecc)
+ cat("chi squared statistic =", chisqstat, "\n")
+ cat("degrees of freedom =", length(ecc) - 1, "\n")
+ pval <- pchisq(chisqstat, length(ecc) - 1, lower.tail = FALSE)
+ cat("p-value =", pval, "\n")
+ if (exists("save.min.pval")) {
+ save.min.pval <<- min(pval, save.min.pval)
+ save.ntests <<- save.ntests + 1
+ } else {
+ save.min.pval <<- pval
+ save.ntests <<- 1
+ }
+ foo <- rbind(cc, ecc)
+ dimnames(foo) <- list(c("observed", "expected"), as.character(xx))
+ print(foo)
+ }
>
> set.seed(42)
> nsim <- 1e4
>
> mu <- 10
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 22)
chi squared statistic = 13.94277
degrees of freedom = 19
p-value = 0.7870174
3 4 5 6 7 8 9
observed 70.00000 195.0000 384.0000 659.0000 881.0000 1146.000 1269.000
expected 75.87668 189.6917 379.3834 632.3057 903.2938 1129.117 1254.575
10 11 12 13 14 15 16
observed 1258.000 1131.000 938.0000 678.0000 513.0000 363.0000 236.0000
expected 1254.575 1140.523 950.4354 731.1042 522.2173 348.1448 217.5905
17 18 19 20 21 22
observed 139.0000 63.00000 41.00000 17.00000 8.000000 11.000000
expected 127.9944 71.10802 37.42527 18.71264 8.910779 7.015935
>
> mu <- 3.5
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 11)
chi squared statistic = 11.41085
degrees of freedom = 8
p-value = 0.1794883
3 4 5 6 7 8 9 10
observed 3208.000 2715.000 1933.000 1127.000 578.0000 296.000 96.000 32.0000
expected 3177.274 2780.115 1946.080 1135.214 567.6068 248.328 96.572 33.8002
11
observed 15.00000
expected 15.00980
>
> mu <- 2.5
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 10)
chi squared statistic = 5.396952
degrees of freedom = 7
p-value = 0.6116408
3 4 5 6 7 8 9
observed 4755.000 2864.000 1464.000 598.0000 214.0000 81.00000 19.00000
expected 4685.865 2928.666 1464.333 610.1387 217.9067 68.09583 18.91551
10
observed 5.000000
expected 6.079791
>
> mu <- 1.5
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 8)
chi squared statistic = 8.190085
degrees of freedom = 5
p-value = 0.1460661
3 4 5 6 7 8
observed 6571.000 2393.000 789.0000 188.0000 50.00000 9.000000
expected 6565.976 2462.241 738.6723 184.6681 39.57173 8.870673
>
> mu <- 0.5
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 6)
chi squared statistic = 1.440299
degrees of freedom = 3
p-value = 0.6961162
3 4 5 6
observed 8803.000 1070.000 118.0000 9.000000
expected 8782.554 1097.819 109.7819 9.845187
>
> nsim <- 1e5
> mu <- 0.1
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 5)
chi squared statistic = 0.04216532
degrees of freedom = 2
p-value = 0.979138
3 4 5
observed 97504.0 2447.000 49.00000
expected 97512.6 2437.815 49.58066
>
> mu <- 0.01
> k <- 2
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 4)
chi squared statistic = 0.3912271
degrees of freedom = 1
p-value = 0.5316551
3 4
observed 99760.00 240.0000
expected 99750.13 249.8749
>
> mu <- 1.5
> xpred <- 0:10
> save.seed <- .Random.seed
> x <- rktp(xpred, k, mu, xpred)
> .Random.seed <- save.seed
> my.x <- rep(0, length(xpred))
> for (i in seq(along = xpred))
+ if (xpred[i] > 0)
+ for (j in 1:xpred[i])
+ my.x[i] <- my.x[i] + rktp(1, k, mu)
> all.equal(x, my.x)
[1] TRUE
>
> k <- 5
> mu <- pi
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 14)
chi squared statistic = 8.99717
degrees of freedom = 8
p-value = 0.3425347
6 7 8 9 10 11 12
observed 58140.00 26293.00 10369.00 3582.000 1187.000 329.0000 81.00000
expected 58366.81 26194.97 10286.74 3590.749 1128.067 322.1752 84.34528
13 14
observed 12.00000 7.000000
expected 20.38296 5.761997
>
> k <- 10
> mu <- exp(2)
> x <- rktp(nsim, k, mu)
> do.chisq.test(x, k, mu, 22)
chi squared statistic = 9.136909
degrees of freedom = 11
p-value = 0.6092566
11 12 13 14 15 16 17
observed 43437.00 26393.00 15018.00 8120.000 3890.000 1823.000 799.0000
expected 43218.94 26612.26 15126.12 7983.409 3932.657 1816.164 789.3963
18 19 20 21 22
observed 326.0000 115.0000 55.00000 16.00000 8.000000
expected 324.0497 126.0222 46.55924 16.38233 8.037665
>
> cat("number of tests:", save.ntests, "\n")
number of tests: 9
> cat("minimum p-value:", save.min.pval, "\n")
minimum p-value: 0.1460661
> cat("Bonferroni corrected minimum p-value:",
+ save.ntests * save.min.pval, "\n")
Bonferroni corrected minimum p-value: 1.314595
>
> #####
>
> set.seed(42)
> nind <- 25
> nnode <- 1
> ncoef <- 1
> k <- 2
>
> pred <- 0
> fam <- 4
>
> theta <- 4 / 3
> mu <- exp(theta)
> x <- rpois(100, mu)
> x <- x[x > k]
> x <- x[1:nind]
>
> modmat <- matrix(1, nrow = nind, ncol = 1)
> out <- mlogl(theta, pred, fam, x, modmat, modmat, deriv = 2,
+ type = "conditional")
> print(out)
$value
[1] -86.3567
$gradient
[1] -14.11330
$hessian
[,1]
[1,] 60.42326
>
> xxx <- seq(0, 100)
> ppp <- dpois(xxx, mu)
> ppp[xxx <= k] <- 0
> ppp <- ppp / sum(ppp)
> tau <- sum(xxx * ppp)
>
> sum(x - tau)
[1] 14.11330
> max(abs( sum(x - tau) + out$gradient ))
[1] 1.776357e-14
>
> length(x) * sum((xxx - tau)^2 * ppp)
[1] 60.42326
>
> epsilon <- 1e-6
> oute <- mlogl(theta + epsilon, pred, fam, x, modmat, modmat, deriv = 2,
+ type = "conditional")
> (oute$value - out$value) / epsilon
[1] -14.11327
> all.equal((oute$value - out$value) / epsilon, out$gradient,
+ tol = 10 * epsilon)
[1] TRUE
>
> (oute$gradient - out$gradient) / epsilon
[1] 60.42331
> all.equal((oute$gradient - out$gradient) / epsilon, as.numeric(out$hessian),
+ tol = 10 * epsilon)
[1] TRUE
>
>
