rnegbinmix: The negative binomial mixture distribution with k components

Description Usage Arguments Details Author(s) Examples

Description

rnegbinmix simulates data from a negative binomial mixture distribution and dnegbinmix computes the density.

Usage

1
2
3
rnegbinmix(n, alpha, beta, cvec, p)

dnegbinmix(k, alpha, beta, cvec, p)

Arguments

n

scalar, positive. number of samples

alpha

scalar, positive

beta

positive

cvec

numeric, positive, length k

p

numeric, probability vector, length k

Details

dnegbinmix can handle non-integer counts (without warning), because the binomial coefficient is implemented with a gamma function.

Author(s)

Johanna Bertl

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
# Test the negative binomial density function and how well it fits to simulated values

k = 0:20
density = dnegbin.alphabeta(k=0:20, alpha=10, beta=5)

simulated = rnegbin.alphabeta(1000, alpha=10, beta=5)
plot(table(simulated)/1000)

lines(k, density, col="purple")


# Same for a mixture of two negative binomial densities

k = 0:50
c1 = 0.5; c2 = 10
p1 = 0.2; p2 = 0.8
alpha = 10
beta = 5

density1 = dnegbin.alphabeta(k = k, alpha = alpha, beta = beta/c1)
plot(k, density1, t="b")
density2 = dnegbin.alphabeta(k = k, alpha = alpha, beta = beta/c2)
plot(k, density2, t="b")
density.mixture = p1 * density1 + p2 * density2
plot(k, density.mixture, t="b")

density.mixture2 = dnegbinmix(k = k, alpha = alpha, beta = beta, cvec = c(c1, c2), p = c(p1, p2))
lines(k, density.mixture2, col="red")

simulated.mixture = rnegbinmix(1000, alpha = alpha, beta = beta, cvec = c(c1, c2), p = c(p1, p2))
plot(table(simulated.mixture)/1000)
lines(k, density.mixture, col="purple")


# Mixture with 3 components

k = 0:500
c1 = 0.01; c2 = 1; c3=100
p1 = 0.33; p2 = 0.33; p3 = 0.34
alpha = 10
beta = 5

density1 = dnegbin.alphabeta(k = k, alpha = alpha, beta = beta/c1)
plot(k, density1, t="l")
density2 = dnegbin.alphabeta(k = k, alpha = alpha, beta = beta/c2)
plot(k, density2, t="l")
density3 = dnegbin.alphabeta(k = k, alpha = alpha, beta = beta/c3)
plot(k, density3, t="l")
density.mixture = p1 * density1 + p2 * density2 + p3 * density3
plot(k, density.mixture, t="l")

simulated.mixture = rnegbinmix(1000, alpha = alpha, beta = beta, cvec = c(c1, c2, c3), p = c(p1, p2, p3))
plot(table(simulated.mixture)/1000)
lines(k, density.mixture, col="purple")

johannabertl/SelectionMix documentation built on May 3, 2019, 4:03 p.m.