-
Notifications
You must be signed in to change notification settings - Fork 0
/
tpmd.R
126 lines (126 loc) · 4.26 KB
/
tpmd.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
#' Power Muth Distribution
#' @export
#' @name tpmd
#' @param x,q vector of quantiles.
#' @param beta a scale 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 alpha a shape parameter.
#' @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 Power Muth distribution with parameters \code{shape} and \code{scale}.
#' @return \code{dtpmd} gives the density, \code{ptpmd} gives the distribution
#' function, \code{qtpmd} gives the quantile function and \code{rtpmd} generates
#' random deviates.
#' @details
#' The Power Muth distribution with \code{shape} parameter \eqn{\alpha} and
#' \code{scale} parameter \eqn{\beta} has density
#' \deqn{f\left( x\right) =\frac{\alpha }{\beta ^\alpha }x^{\alpha -1}
#' \left( e^{\left(x/\beta \right) ^{\alpha }}-1\right)
#' \left( e^{\left( x/\beta \right) ^{\alpha }-
#' \left( e^{\left( x/\beta \right) ^{\alpha }}-1\right) }\right), }
#' where
#' \deqn{x>0,~\alpha ,\beta>0.}
#' @references Jodra, P., Gomez, H. W., Jimenez-Gamero,
#' M. D., and Alba-Fernandez, M. V. (2017).
#' *The power Muth distribution* . Mathematical Modelling and Analysis, 22(2),
#' 186-201.
#' @note
#' Hazard function;
#' \deqn{h\left( \beta ,\alpha \right) =\frac{\alpha }{\beta ^{\alpha }}
#' \left(e^{\left( x/\beta \right) ^{\alpha }}-1\right) x^{\alpha -1}}
#' @examples
#' library(new.dist)
#' dtpmd(1, beta=2, alpha=3)
dtpmd<-function(x,beta=1,alpha,log=FALSE)
{
if(any(beta<=0)) {stop("beta must be > 0")}
if(any(alpha<=0)) {stop("alpha must be > 0")}
enuzun <- max(length(x), 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]
pdf<-NULL
for (i in 1:enuzun)
{
if(x[i]<=0) {pdf[i]<-0} else
{
pdf[i]<-(alpha[i]/(beta[i]^alpha[i]))*(x[i]^(alpha[i]-1))*
(exp((x[i]/beta[i])^alpha[i])-1)*exp(((x[i]/beta[i])^alpha[i])-
(exp((x[i]/beta[i])^alpha[i])-1))
}
}
if(log==TRUE) pdf<-log(pdf)
return(pdf)
}
#' Power Muth Distribution
#' @export
#' @rdname tpmd
#' @examples
#' ptpmd(1,beta=2,alpha=3)
ptpmd<-function(q,beta=1,alpha,lower.tail=TRUE,log.p=FALSE)
{
if(any(beta<=0)) {stop("beta must be > 0")}
if(any(alpha<=0)) {stop("alpha must be > 0")}
enuzun<-max(length(q), length(alpha), length(beta))
q<-rep(q,enuzun/length(q)+1)[1:enuzun]
alpha<-rep(alpha, enuzun/length(alpha)+1)[1:enuzun]
beta<-rep(beta,enuzun/length(beta)+1)[1:enuzun]
cdf<-NULL
for (i in 1:enuzun)
{
if(q[i]>0) cdf[i]<-1-exp(((q[i]/beta[i])^alpha[i])-(exp((q[i]/
beta[i])^alpha[i])-1)) else cdf[i]<-0
}
if(lower.tail==FALSE) cdf<-1-cdf
if(log.p==TRUE) cdf<-log(cdf)
return(cdf)
}
#' Power Muth Distribution
#' @export
#' @rdname tpmd
#' @examples
#' qtpmd(.5,beta=2,alpha=3)
qtpmd<-function(p,beta=1,alpha,lower.tail=TRUE)
{
if(any(p<0)|any(p>1)) {stop("p must be between >= 0 and <= 1")}
if(any(beta<=0)) {stop("beta must be > 0")}
if(any(alpha<=0)) {stop("alpha must be > 0")}
enuzun <- max(length(p), length(alpha), length(beta))
p<-rep(p,enuzun/length(p)+1)[1:enuzun]
alpha<-rep(alpha, enuzun/length(alpha)+1)[1:enuzun]
beta<-rep(beta,enuzun/length(beta)+1)[1:enuzun]
quant<-NULL
for (i in 1:enuzun)
{
quant[i]<-exp(log(-pracma::lambertWn((-1+p[i])/exp(1))-1+log(1-p[i]))/
alpha[i])*beta[i]
}
if(lower.tail==FALSE)
{
quant[i]<-(exp(log(-pracma::lambertWn((-1+(1-p[i]))/exp(1))-1+
log(1-(1-p[i])))/alpha[i])*beta[i])
}
{
return(quant)
}
}
#' Power Muth Distribution
#' @export
#' @rdname tpmd
#' @examples
#' rtpmd(10,beta=2,alpha=3)
rtpmd<-function(n,beta=1,alpha)
{
n<-floor(n)
if(any(n<1)) {stop("n must be >= 1")}
if(any(beta<=0)) {stop("beta must be > 0")}
if(any(alpha<=0)) {stop("alpha must be > 0")}
suppressWarnings({
rn<-qtpmd(stats::runif(n),beta,alpha)})
return(rn)
}