rnegbinmix: The negative binomial mixture distribution with k components

In johannabertl/SelectionMix: Negative binomial mixture model

Description

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

Usage

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

# 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.