The gamma parameter if known (NOT IMPLEMENTED YET).
prob
The prob parameter if known (NOT IMPLEMENTED YET).
w
Vector of time length (durations).
sum.w
NOT IMPLEMENTED YET. The effective duration. If sum.w is
strictly inferior to sum(w), it is to be understood that
missing periods occur within the counts period. This can be taken
into account with a suitable algorithm (Expectation Maximisation,
etc.)
interval
Interval giving min and max values for gamma.
optim
If TRUE a one-dimensional optimisation is used. Else the zero
of the derivative of the (concentrated) log-likelihood is searched
for.
plot
Should a plot be drawn? May be removed in the future.
NegBinomial for the negative binomial distribution,
glm.nb from the MASS package for fitting Generalised
Linear Model of the negative binomial family.
Examples
## known parameters
nint <- 100
gam <- 6; prob <- 0.20
## draw random w, then the counts N
w <- rgamma(nint, shape = 3, scale = 1/5)
N <- rnbinom(nint, size = w * gam, prob = prob)
mu <- w * gam * (1 - prob) / prob
Res <- NBlevy(N = N, w = w)
## Use example data 'Brest'
## compute the number of event and duration of the non-skipped periods
gof1 <- gof.date(date = Brest$OTdata$date,
skip = Brest$OTmissing,
start = Brest$OTinfo$start,
end = Brest$OTinfo$end,
plot.type = "omit")
ns1 <- gof1$noskip
## fit the NBlevy
fit1 <- NBlevy(N = ns1$nevt, w = ns1$duration)
## use a higher threshold
OT2 <- subset(Brest$OTdata, Surge > 50)
gof2 <- gof.date(date = OT2$date,
skip = Brest$OTmissing,
start = Brest$OTinfo$start,
end = Brest$OTinfo$end,
plot.type = "omit")
ns2 <- gof2$noskip
## the NBlevy prob is now closer to 1
fit2 <- NBlevy(N = ns2$nevt, w = ns2$duration)
c(fit1$prob, fit2$prob)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Renext)
Loading required package: evd
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Renext/NBlevy.Rd_%03d_medium.png", width=480, height=480)
> ### Name: NBlevy
> ### Title: Negative Binomial Levy process
> ### Aliases: NBlevy
>
> ### ** Examples
>
> ## known parameters
> nint <- 100
> gam <- 6; prob <- 0.20
>
> ## draw random w, then the counts N
> w <- rgamma(nint, shape = 3, scale = 1/5)
> N <- rnbinom(nint, size = w * gam, prob = prob)
> mu <- w * gam * (1 - prob) / prob
> Res <- NBlevy(N = N, w = w)
>
> ## Use example data 'Brest'
> ## compute the number of event and duration of the non-skipped periods
> gof1 <- gof.date(date = Brest$OTdata$date,
+ skip = Brest$OTmissing,
+ start = Brest$OTinfo$start,
+ end = Brest$OTinfo$end,
+ plot.type = "omit")
> ns1 <- gof1$noskip
> ## fit the NBlevy
> fit1 <- NBlevy(N = ns1$nevt, w = ns1$duration)
>
> ## use a higher threshold
> OT2 <- subset(Brest$OTdata, Surge > 50)
> gof2 <- gof.date(date = OT2$date,
+ skip = Brest$OTmissing,
+ start = Brest$OTinfo$start,
+ end = Brest$OTinfo$end,
+ plot.type = "omit")
> ns2 <- gof2$noskip
> ## the NBlevy prob is now closer to 1
> fit2 <- NBlevy(N = ns2$nevt, w = ns2$duration)
>
> c(fit1$prob, fit2$prob)
[1] 0.3203404 0.5162115
>
>
>
>
>
> dev.off()
null device
1
>