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
|
/* actuar: Actuarial Functions and Heavy Tailed Distributions
*
* Functions to compute density, cumulative distribution and quantile
* functions, and to simulate random variates for the loggamma
* distribution. See ../R/Loggamma.R for details.
*
* AUTHORS: Mathieu Pigeon and Vincent Goulet <vincent.goulet@act.ulaval.ca>
*/
#include <R.h>
#include <Rmath.h>
#include "locale.h"
#include "dpq.h"
double dlgamma(double x, double shapelog, double ratelog, int give_log)
{
#ifdef IEEE_754
if (ISNAN(x) || ISNAN(shapelog) || ISNAN(ratelog))
return x + shapelog + ratelog;
#endif
if (!R_FINITE(shapelog) ||
!R_FINITE(ratelog) ||
shapelog <= 0.0 ||
ratelog < 0.0)
return R_NaN;;
if (!R_FINITE(x) || x < 1.0)
return ACT_D__0;
return ACT_D_exp(dgamma(log(x), shapelog, 1.0/ratelog, 1) - log(x));
}
double plgamma(double q, double shapelog, double ratelog, int lower_tail,
int log_p)
{
#ifdef IEEE_754
if (ISNAN(q) || ISNAN(shapelog) || ISNAN(ratelog))
return q + shapelog + ratelog;
#endif
if (!R_FINITE(shapelog) ||
!R_FINITE(ratelog) ||
shapelog <= 0.0 ||
ratelog < 0.0)
return R_NaN;;
if (q <= 1.0)
return ACT_DT_0;
return pgamma(log(q), shapelog, 1.0/ratelog, lower_tail, log_p);
}
double qlgamma(double p, double shapelog, double ratelog, int lower_tail,
int log_p)
{
#ifdef IEEE_754
if (ISNAN(p) || ISNAN(shapelog) || ISNAN(ratelog))
return p + shapelog + ratelog;
#endif
if (!R_FINITE(shapelog) ||
!R_FINITE(ratelog) ||
shapelog <= 0.0 ||
ratelog <= 0.0)
return R_NaN;;
ACT_Q_P01_boundaries(p, 1, R_PosInf);
p = ACT_D_qIv(p);
return exp(qgamma(p, shapelog, 1.0/ratelog, lower_tail, 0));
}
double rlgamma(double shapelog, double ratelog)
{
if (!R_FINITE(shapelog) ||
!R_FINITE(ratelog) ||
shapelog <= 0.0 ||
ratelog <= 0.0)
return R_NaN;;
return exp(rgamma(shapelog, 1.0/ratelog));
}
double mlgamma(double order, double shapelog, double ratelog, int give_log)
{
#ifdef IEEE_754
if (ISNAN(order) || ISNAN(shapelog) || ISNAN(ratelog))
return order + shapelog + ratelog;
#endif
if (!R_FINITE(shapelog) ||
!R_FINITE(ratelog) ||
!R_FINITE(order) ||
shapelog <= 0.0 ||
ratelog <= 0.0)
return R_NaN;
if (order >= ratelog)
return R_PosInf;
return R_pow(1.0 - order / ratelog, -shapelog);
}
double levlgamma(double limit, double shapelog, double ratelog, double order,
int give_log)
{
#ifdef IEEE_754
if (ISNAN(limit) || ISNAN(shapelog) || ISNAN(ratelog) || ISNAN(order))
return limit + shapelog + ratelog + order;
#endif
if (!R_FINITE(shapelog) ||
!R_FINITE(ratelog) ||
!R_FINITE(limit) ||
!R_FINITE(order) ||
shapelog <= 0.0 ||
ratelog <= 0.0 ||
limit <= 0.0)
return R_NaN;
if (order >= ratelog)
return R_PosInf;
if (limit <= 1.0)
return 0.0;
double u = log(limit);
return R_pow(1.0 - order / ratelog, -shapelog)
* pgamma(u * (ratelog - order), shapelog, 1.0, 1, 0)
+ ACT_DLIM__0(limit, order) * pgamma(u * ratelog, shapelog, 1.0, 0, 0);
}
|