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, alpha, k, mu, max.bin) {
+      stopifnot(all(x > k))
+      stopifnot(k + 1 < max.bin)
+      xx <- seq(k + 1, max.bin)
+      yy <- dnbinom(xx, size = alpha, mu = mu)
+      yy[length(yy)] <- pnbinom(max.bin - 1, size = alpha, mu = 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
> 
>  alpha <- 2.222
>  mu <- 10
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 40)
chi squared statistic = 53.69477 
degrees of freedom = 37 
p-value = 0.03731892 
                3        4       5        6        7        8        9       10
observed 707.0000 724.0000 799.000 798.0000 742.0000 637.0000 597.0000 547.0000
expected 708.2678 756.5403 770.282 758.6016 729.0392 687.6125 638.9911 586.7091
               11       12       13       14      15       16       17       18
observed 522.0000 443.0000 410.0000 382.0000 336.000 326.0000 301.0000 235.0000
expected 533.3719 480.8435 430.4065 382.8955 338.806 298.3818 261.6840 228.6446
               19       20       21       22      23       24      25       26
observed 215.0000 183.0000 135.0000 140.0000 115.000 100.0000 95.0000 73.00000
expected 199.1082 172.8635 149.6666 129.2586 111.378  95.7691 82.1881 70.40659
               27       28       29       30       31       32       33
observed 74.00000 51.00000 39.00000 36.00000 40.00000 36.00000 18.00000
expected 60.21367 51.41676 43.84172 37.33230 31.74924 26.96912 22.88316
               34       35       36       37       38       39       40
observed 28.00000 12.00000 18.00000 17.00000 13.00000 11.00000 45.00000
expected 19.39585 16.42370 13.89395 11.74344  9.91743  8.36866 44.13479
> 
>  alpha <- 2.222
>  mu <- 3.5
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 20)
chi squared statistic = 15.95459 
degrees of freedom = 17 
p-value = 0.5270567 
                3        4        5        6        7        8        9
observed 2575.000 2045.000 1558.000 1152.000 860.0000 555.0000 411.0000
expected 2572.008 2053.853 1563.326 1151.002 826.9431 583.0838 405.0835
               10       11       12       13       14       15       16
observed 257.0000 205.0000 140.0000 82.00000 55.00000 41.00000 19.00000
expected 278.0577 188.9752 127.3623 85.22723 56.68162 37.49519 24.68649
               17       18       19       20
observed 16.00000 17.00000 3.000000  9.00000
expected 16.18552 10.57238 6.882776 12.57426
> 
>  alpha <- 2.222
>  mu <- 2.5
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 16)
chi squared statistic = 6.979199 
degrees of freedom = 13 
p-value = 0.9032203 
                3        4        5        6        7        8        9
observed 3465.000 2392.000 1565.000  983.000 634.0000 395.0000 247.0000
expected 3462.798 2393.415 1576.856 1004.876 624.8932 381.3774 229.3308
               10       11       12       13       14       15       16
observed 137.0000 82.00000 40.00000 31.00000 12.00000 10.00000  7.00000
expected 136.2532 80.15125 46.75631 27.08143 15.58940  8.92599 11.69563
> 
>  alpha <- 2.222
>  mu <- 1.5
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 12)
chi squared statistic = 6.382718 
degrees of freedom = 9 
p-value = 0.7010869 
                3        4        5        6        7        8        9
observed 4877.000 2621.000 1286.000 630.0000 329.0000 142.0000 69.00000
expected 4925.874 2591.642 1299.720 630.4796 298.4457 138.6486 63.46345
               10       11       12
observed 23.00000 11.00000 12.00000
expected 28.70178 12.85208 10.17242
> 
>  alpha <- 2.222
>  mu <- 0.5
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 8)
chi squared statistic = 4.75385 
degrees of freedom = 5 
p-value = 0.4466531 
                3        4       5        6        7        8
observed 7567.000 1863.000 445.000 93.00000 26.00000 6.000000
expected 7633.273 1830.499 418.419 92.51232 19.96002 5.336489
> 
>  alpha <- 2.222
>  mu <- 0.1
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 5)
chi squared statistic = 6.224807 
degrees of freedom = 2 
p-value = 0.04449389 
              3        4        5
observed 9424.0 558.0000 18.00000
expected 9439.3 530.7065 29.99429
> 
>  nsim <- 2e5
>  alpha <- 2.222
>  mu <- 0.01
>  k <- 2
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 5)
chi squared statistic = 3.500094 
degrees of freedom = 2 
p-value = 0.1737658 
                3        4       5
observed 198791.0 1206.000 3.00000
expected 198830.5 1162.963 6.51898
> 
>  alpha <- 2.222
>  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
> 
>  nsim <- 1e4
>  alpha <- 5.55
>  k <- 5
>  mu <- pi
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 16)
chi squared statistic = 9.004305 
degrees of freedom = 10 
p-value = 0.531695 
                6        7        8        9       10       11       12
observed 4321.000 2536.000 1495.000 782.0000 452.0000 207.0000 117.0000
expected 4324.703 2579.234 1462.497 795.8703 418.5581 213.8671 106.6133
               13       14       15       16
observed 54.00000 16.00000 11.00000 9.000000
expected 52.02302 24.91506 11.73729 9.981887
> 
>  alpha <- 5.55
>  k <- 10
>  mu <- exp(2)
>  x <- rktnb(nsim, alpha, k, mu)
>  do.chisq.test(x, alpha, k, mu, 29)
chi squared statistic = 30.63213 
degrees of freedom = 18 
p-value = 0.03173578 
               11       12       13       14       15       16      17       18
observed 2418.000 1866.000 1514.000 1168.000 867.0000 642.0000 438.000 310.0000
expected 2466.218 1942.382 1497.458 1133.071 843.3325 618.5532 447.777 320.3481
               19       20       21       22      23       24       25       26
observed 244.0000 151.0000 141.0000 82.00000 49.0000 38.00000 23.00000 12.00000
expected 226.7492 158.9475 110.4363 76.10968 52.0619 35.36725 23.87290 16.01878
               27       28       29
observed  5.00000 8.000000 24.00000
expected 10.68935 7.096279 13.51152
> 
>  cat("number of tests:", save.ntests, "\n")
number of tests: 9 
>  cat("minimum p-value:", save.min.pval, "\n")
minimum p-value: 0.03173578 
>  cat("Bonferroni corrected minimum p-value:",
+      save.ntests * save.min.pval, "\n")
Bonferroni corrected minimum p-value: 0.285622 
> 
>  #####
> 
>  set.seed(42)
>  nind <- 25
>  nnode <- 1
>  ncoef <- 1
>  alpha <- 3.333
>  k <- 2
> 
>  pred <- 0
>  fam <- 1
>  ifam <- fam.truncated.negative.binomial(size = alpha, trunc = k)
>  aster:::setfam(list(ifam))
>  theta.origin <- aster:::getfam()[[1]]$origin
> 
>  theta <- (- 4 / 3)
>  p <- 1 - exp(theta)
>  x <- rnbinom(1000, size = alpha, prob = p)
>  x <- x[x > k]
>  x <- x[1:nind]
>  modmat <- matrix(1, nrow = nind, ncol = 1)
> 
>  out <- mlogl(theta - theta.origin, pred, fam, x, modmat, modmat,
+      deriv = 2, type = "conditional", famlist = list(ifam))
>  print(out)
$value
[1] 84.61333

$gradient
[1] 9.847512

$hessian
         [,1]
[1,] 23.66414

> 
>  xxx <- seq(0, 100)
>  ppp <- dnbinom(xxx, size = alpha, prob = p)
>  ppp[xxx <= k] <- 0
>  ppp <- ppp / sum(ppp)
>  tau <- sum(xxx * ppp)
> 
>  my.grad.logl <- sum(x - tau)
>  all.equal(- out$gradient, my.grad.logl)
[1] TRUE
> 
>  my.fish.info <- length(x) * sum((xxx - tau)^2 * ppp)
>  all.equal(as.numeric(out$hessian), my.fish.info)
[1] TRUE
> 
>