← Back to branch summary

~ubuntu-branches/ubuntu/vivid/r-cran-statmod/vivid

« back to all changes in this revision

Viewing changes to R/invgauss.R

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-10-02 15:40:39 UTC
  • Revision ID: package-import@ubuntu.com-20131002154039-lee9ml5cj2tkgjga
Tags: upstream-1.4.18
Import upstream version 1.4.18

Show diffs side-by-side

added added

removed removed

Lines of Context:
R/invgauss.R
 
1
dinvgauss <- function(x, mu, lambda = 1, log=FALSE)
 
2
#  Density of inverse Gaussian distribution
 
3
#  Gordon Smyth
 
4
#       15 Jan 1998.  Last revised 19 June 2009.
 
5
{
 
6
        if(any(mu<=0)) stop("mu must be positive")
 
7
        if(any(lambda<=0)) stop("lambda must be positive")
 
8
        d <- rep(-Inf,length(x))
 
9
        i <- x>0
 
10
        d[i] <- (log(lambda)-log(2*pi)-3*log(x[i]))/2-lambda*(x[i]-mu)^2/(2*mu^2*x[i])
 
11
        if(!log) d <- exp(d)
 
12
        if(!is.null(Names <- names(x))) names(d) <- rep(Names, length = length(d))
 
13
        d
 
14
}
 
15
 
 
16
pinvgauss <- function(q, mu, lambda = 1)
 
17
{
 
18
#  Inverse Gaussian distribution function
 
19
#  GKS  15 Jan 98
 
20
#
 
21
        if(any(mu<=0)) stop("mu must be positive")
 
22
        if(any(lambda<=0)) stop("lambda must be positive")
 
23
        n <- length(q)
 
24
        if(length(mu)>1 && length(mu)!=n) mu <- rep(mu,length=n)
 
25
        if(length(lambda)>1 && length(lambda)!=n) lambda <- rep(lambda,length=n)
 
26
        lq <- sqrt(lambda/q)
 
27
        qm <- q/mu
 
28
        p <- ifelse(q>0,pnorm(lq*(qm-1))+exp(2*lambda/mu)*pnorm(-lq*(qm+1)),0)
 
29
        if(!is.null(Names <- names(q)))
 
30
                names(p) <- rep(Names, length = length(p))
 
31
        p
 
32
}
 
33
 
 
34
rinvgauss <- function(n, mu, lambda = 1)
 
35
#       Random variates from inverse Gaussian distribution
 
36
#       Reference:  Chhikara and Folks, The Inverse Gaussian
 
37
#       Distribution, Marcel Dekker, 1989, page 53.
 
38
#       Gordon Smyth 15 Jan 98.
 
39
#       Revised by Trevor Park 14 June 2005.
 
40
{
 
41
        if(any(mu<=0)) stop("mu must be positive")
 
42
        if(any(lambda<=0)) stop("lambda must be positive")
 
43
        if(length(n)>1) n <- length(n)
 
44
        if(length(mu)>1 && length(mu)!=n) mu <- rep(mu,length=n)
 
45
        if(length(lambda)>1 && length(lambda)!=n) lambda <- rep(lambda,length=n)
 
46
        y2 <- rchisq(n,1)
 
47
        u <- runif(n)
 
48
        r2 <- mu/(2*lambda)*(2*lambda+mu*y2+sqrt(4*lambda*mu*y2+mu^2*y2^2))
 
49
        r1 <- mu^2/r2
 
50
        ifelse(u < mu/(mu+r1), r1, r2)
 
51
}
 
52
 
 
53
qinvgauss  <- function(p, mu, lambda = 1)
 
54
{
 
55
#  Quantiles of the inverse Gaussian distribution
 
56
#  Dr Paul Bagshaw
 
57
#  Centre National d'Etudes des Telecommunications (DIH/DIPS)
 
58
#  Technopole Anticipa, France
 
59
#  paul.bagshaw@cnet.francetelecom.fr
 
60
#  23 Dec 98
 
61
#
 
62
  if(any(mu <= 0.))
 
63
    stop("mu must be positive")
 
64
  if(any(lambda <= 0.))
 
65
    stop("lambda must be positive")
 
66
  n <- length(p)
 
67
  if(length(mu) > 1 && length(mu) != n)
 
68
    mu <- rep(mu, length = n)
 
69
  if(length(lambda) > 1 && length(lambda) != n)
 
70
    lambda <- rep(lambda, length = n)
 
71
  thi <- lambda / mu
 
72
  U <- qnorm (p)
 
73
  r1 <- 1 + U / sqrt (thi) + U^2 / (2 * thi) + U^3 / (8 * thi * sqrt(thi))
 
74
  x <- r1
 
75
  for (i in 1:10) {
 
76
    cum <- pinvgauss (x, 1., thi)
 
77
    dx <- (cum - p) / dinvgauss (x, 1., thi)
 
78
    dx <- ifelse (is.finite(dx), dx, ifelse (p > cum, -1, 1))
 
79
    dx[dx < -1] <- -1
 
80
    if (all(dx == 0.)) break
 
81
    x <- x - dx
 
82
  }
 
83
  x * mu
 
84
}