R/theta.R

   1 ## theta.R (2002-08-28)
   2
   3 ##   Population Parameter THETA
   4
   5 ## theta.h: using homozigosity
   6 ## theta.k: using expected number of alleles
   7 ## theta.s: using segregating sites in DNA sequences
   8
   9 ## Copyright 2002 Emmanuel Paradis
  10
  11 ## This file is part of the R-package `ape'.
  12 ## See the file ../COPYING for licensing issues.
  13
  14 theta.h <- function(x, standard.error = FALSE)
  15 {
  16     HE <- H(x, variance = TRUE)
  17     sdH <- HE[2]
  18     HE <- HE[1]
  19     f <- function(th) HE - th * (1 + (2 * (1 + th)) / ((2 + th) * (3 + th)))
  20     th <- uniroot(f, interval = c(0, 1))$root
  21     if (standard.error) {
  22         SE <- (2 + th)^2 * (2 + th)^3 * sdH /
  23           HE^2 * (1 + th) * ((2 + th) * (3 + th) * (4 + th) + 10 * (2 + th) + 4)
  24         return(c(th, SE))
  25     }
  26     else return(th)
  27 }
  28
  29 theta.k <- function(x, n = NULL, k = NULL)
  30 {
  31     if (is.null(n)) {
  32         if (!is.factor(x)) {
  33             if (is.numeric(x)) {
  34                 n <- sum(x)
  35                 k <- length(x)
  36             }
  37             else x <- factor(x)
  38         }
  39         if (is.factor(x)) { # ne pas remplacer par `else'...
  40             n <- length(x)
  41             k <- nlevels(x)
  42         }
  43     }
  44     f <- function(th) th * sum(1 / (th + (0:(n - 1)))) - k
  45     th <- uniroot(f, interval = c(1e-8, 100))$root
  46     return(th)
  47 }
  48
  49 theta.s <- function(s, n, variance = FALSE)
  50 {
  51     a1 <- sum(1 / (1:(n - 1)))
  52     th <- s / a1
  53     if (variance) {
  54         a2 <- sum(1 / (1:(n - 1))^2)
  55         var.th <- (a1^2 * s + a2 * s^2) / (a1^2 * (a1^2 + a2))
  56         return(c(th, var.th))
  57     }
  58     else return(th)
  59 }