-
Notifications
You must be signed in to change notification settings - Fork 0
/
sgrd.R
138 lines (138 loc) · 5.38 KB
/
sgrd.R
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#' Slashed Generalized Rayleigh Distribution
#' @export
#' @name sgrd
#' @param x,q vector of quantiles.
#' @param theta a scale parameter.
#' @param alpha a shape parameter.
#' @param beta a kurtosis parameter.
#' @param p vector of probabilities.
#' @param n number of observations. If \code{length(n) > 1}, the length is taken
#' to be the number required.
#' @param log,log.p logical; if TRUE, probabilities p are given as log(p).
#' @param lower.tail logical; if TRUE (default), probabilities are
#' \eqn{P\left[ X\leq x\right]}, otherwise, \eqn{P\left[ X>x\right] }.
#' @description
#' Density, distribution function, quantile function and random generation for
#' the Slashed generalized Rayleigh distribution with parameters \code{shape},
#' \code{scale} and \code{kurtosis}.
#' @return \code{dsgrd} gives the density, \code{psgrd} gives the distribution
#' function, \code{qsgrd} gives the quantile function and \code{rsgrd} generates
#' random deviates.
#' @details
#' The Slashed Generalized Rayleigh distribution with \code{shape} parameter
#' \eqn{\alpha}, \code{scale} parameter \eqn{\theta} and \code{kurtosis}
#' parameter \eqn{\beta}, has density
#' \deqn{f\left( x\right) =\frac{\beta x^{-\left( \beta+1\right)}}{\Gamma \left(
#' \alpha+1\right) \theta ^{\beta/2}}\Gamma \left( \frac{2\alpha +\beta +2}{2}
#' \right)F\left( \theta x^{2};\frac{2\alpha +\beta +2}{2},1\right), }
#' where F(.;a,b) is the cdf of the Gamma (a,b) distribution, and
#' \deqn{x>0,~\theta >0,~\alpha >-1~and~\beta >0}
#' @references Iriarte, Y. A., Vilca, F., Varela, H. and Gómez, H. W., 2017,
#' *Slashed generalized Rayleigh distribution*, Communications in Statistics-
#' Theory and Methods, 46 (10), 4686-4699.
#' @examples
#' library(new.dist)
#' dsgrd(2,theta=3,alpha=1,beta=4)
dsgrd <- function(x,theta,alpha,beta,log=FALSE)
{
if(any(alpha<=-1)) {stop("alpha must be > -1")}
if(any(theta<=0)) {stop("theta must be > 0")}
if(any(beta<=2)) {(stop("the function may not work when beta<2 since the
first and second moment does not exist!"))}
enuzun<-max(length(x),length(theta),length(alpha),length(beta))
x<-rep(x,enuzun/length(x)+1)[1:enuzun]
alpha<-rep(alpha,enuzun/length(alpha)+1)[1:enuzun]
beta<-rep(beta,enuzun/length(beta)+1)[1:enuzun]
theta<-rep(theta,enuzun/length(theta)+1)[1:enuzun]
pdf<-NULL
for (i in 1:enuzun)suppressWarnings(
{
if(x[i]<=0) {pdf[i]<-0} else
{pdf[i]<-beta[i]*x[i]^(-(beta[i]+1))/(gamma(alpha[i]+1)*
theta[i]^(beta[i]/2))*gamma((2*alpha[i]+beta[i]+2)/2)*
stats::pgamma(theta[i]*x[i]^2,(2*alpha[i]+beta[i]+2)/2,1)}
})
if(log==TRUE) pdf<-log(pdf)
return(pdf)
}
#' Slashed Generalized Rayleigh Distribution
#' @export
#' @rdname sgrd
#' @examples
#' psgrd(5,theta=3,alpha=1,beta=4)
psgrd<-function(q,theta,alpha,beta,lower.tail=TRUE,log.p=FALSE)
{
if(any(alpha<=-1)) {stop("alpha must be > -1")}
if(any(theta<=0)) {stop("theta must be > 0")}
if(any(beta<=2)) {(stop("the function may not work when beta<2 since the
first and second moment does not exist!"))}
enuzun <- max(length(q),length(theta),length(alpha),length(beta))
q<-rep(q,enuzun/length(q)+1)[1:enuzun]
theta<-rep(theta,enuzun/length(theta)+1)[1:enuzun]
alpha<-rep(alpha, enuzun/length(alpha)+1)[1:enuzun]
beta<-rep(beta,enuzun/length(beta)+1)[1:enuzun]
integral<-NULL
for (i in 1:enuzun)suppressWarnings(
{
if(q[i]>0){integral[i] <- stats::integrate(function(t){beta[i]*
t^(-(beta[i]+1))/(gamma(alpha[i]+1)*theta[i]^(beta[i]/2))*
gamma((2*alpha[i]+beta[i]+2)/2)*
stats::pgamma(theta[i]*t^2,(2*alpha[i]+beta[i]+2)/2,1)
},0,q[i])$value} else {integral[i]<-0}
})
if(lower.tail==FALSE) integral <- 1-integral
if(log.p==TRUE) integral <- log(integral)
return(integral)
}
#' Slashed Generalized Rayleigh Distribution
#' @export
#' @rdname sgrd
#' @examples
#' qsgrd(.4,theta=3,alpha=1,beta=4)
qsgrd<-function(p,theta,alpha,beta,lower.tail=TRUE)
{
if(any(p<0)|any(p>1)) {stop("p must be between >= 0 and <= 1")}
if(any(alpha<=-1)) {stop("alpha must be > -1")}
if(any(theta<=0)) {stop("theta must be > 0")}
if(any(beta<=2)) {(stop("the function may not work when beta<2 since
the first and second moment does not exist!"))}
enuzun <- max(length(p),length(theta),length(alpha),length(beta))
theta<-rep(theta,enuzun/length(theta)+1)[1:enuzun]
alpha<-rep(alpha, enuzun/length(alpha)+1)[1:enuzun]
beta<-rep(beta,enuzun/length(beta)+1)[1:enuzun]
kok<-NULL
for (i in 1:enuzun)suppressWarnings(
{
Ex<-(gamma(1/2+alpha[i]+1)*beta[i])/(gamma(alpha[i]+1)*theta[i]^(1/2)*
(beta[i]-1))
Y<-function(t)
{
abs(psgrd(t,theta[i],alpha[i],beta[i])-p[i])
}
if(lower.tail==FALSE)
{
Y<-function(t)
{
abs(psgrd(t,theta[i],alpha[i],beta[i])-(1-p[i]))
}
}
kok[i]<-(stats::optim(Ex,Y,method="Brent",lower=0,upper=20))$par
})
return(kok)
}
#' Slashed Generalized Rayleigh Distribution
#' @export
#' @rdname sgrd
#' @examples
#' rsgrd(10,theta=3,alpha=1,beta=4)
rsgrd<-function(n,theta,alpha,beta)
{
n<-floor(n)
if(any(n<1)) {stop("n must be >= 1")}
if(any(alpha<=-1)) {stop("alpha must be > -1")}
if(any(theta<=0)) {stop("theta must be > 0")}
if(any(beta<=2)) {(stop("the function may not work when beta<2 since
the first and second moment does not exist!"))}
rn<-qsgrd(stats::runif(n),theta,alpha,beta)
return(rn)
}