cran
如题,我要做一个Zero inflated negative binomial model
argument:
p--- prob of zero
alpha
mu
这个是我做的
########R code##########
r.zinb <- function(p,alpha,mu) {
N <- 100
pr <- rep(0,(N+1))
sc <- 100000000
for (i in 0:N){
pr[i+1] <- sc*p.zinb(p,i,alpha,mu)
}
cp <- cumsum(pr)
ran <- runif(1,0,sc)
i <- 1
while(ran >= cp ) {
i <- i + 1
}
return (i-1)
}
###############
还有一个ZINB的pdf
######################
p.zinb <- function(p,y,alpha,mu) {
if (y == 0)
return(p+(1-p)*exp(-alpha*log(1+mu/alpha)))
return((1-p)*exp(y*log(mu/alpha)-(y+alpha)*log(1+mu/alpha)
+lgamma(y+alpha)-lgamma(alpha)-lgamma(y+1)))
}
###############################
由于技术有限,我只能做出这样的。请rtist等高手指导一下,谢谢。
cran
马的,弄了一晚上,今天我supervisor给我一个超级简单以及快的方法
present <- 1-rbinom(100,1,0.2)
library(MASS)
temp <- rnegbin(100,1,0.8)
y <- present * temp
这样就好了