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
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
/* actuar: Actuarial Functions and Heavy Tailed Distributions
*
* Functions to generate variates of phase-type distributions. This
* file is based on random.c with the following modifications:
*
* 1. support for a matrix argument;
* 2. no iteration over the parameters;
* 3. support for two parameter distributions only;
* 4. no support for integer random variates.
*
* For details, see random.c.
*
* AUTHOR: Vincent Goulet <vincent.goulet@act.ulaval.ca>
*/
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Memory.h>
#include "actuar.h"
#include "locale.h"
/* Prototypes of auxiliary function */
static Rboolean randomphtype2(double (*f)(double *, double **, double *, int),
double *, double *, int, double *, int);
static Rboolean randomphtype2(double (*f)(double *, double **, double *, int),
double *a, double *b, int na, double *x, int n)
{
int i, j;
double *rates, **Q;
Rboolean naflag = FALSE;
/* The sub-intensity matrix and initial probability vector never
* change, so compute the transition matrix of the underlying
* Markov chain and the vector of rate parameters before
* looping. */
rates = (double *) R_alloc(na, sizeof(double));
Q = (double **) R_alloc(na, sizeof(double));
for (i = 0; i < na; i++)
{
Q[i] = (double *) S_alloc(na, sizeof(double));
rates[i] = -b[i * (na + 1)];
for (j = 0; j < na; j++)
if (i != j)
Q[i][j] = b[i + j * na] / rates[i];
}
for (i = 0; i < n; i++)
{
x[i] = f(a, Q, rates, na);
if (!R_FINITE(x[i])) naflag = TRUE;
}
return(naflag);
}
#define RANDPHTYPE2(num, fun) \
case num: \
randomphtype2(fun, REAL(a), REAL(b), na, REAL(x), n); \
break
/* The function below retains a 'type' argument that is not actually
* used. This is to fit within the scheme of the other random
* generation functions of random.c and names.c. */
SEXP actuar_do_randomphtype2(int code, SEXP args, SEXPTYPE type /* unused */)
{
SEXP x, a, b, bdims;
int i, n, na, nrow, ncol;
Rboolean naflag = FALSE;
/* Check validity of arguments */
if (!isVector(CAR(args)) ||
!isNumeric(CADR(args)) ||
!isMatrix(CADDR(args)))
error(_("invalid arguments"));
/* Number of variates to generate */
if (LENGTH(CAR(args)) == 1)
{
n = asInteger(CAR(args));
if (n == NA_INTEGER || n < 0)
error(_("invalid arguments"));
}
else
n = LENGTH(CAR(args));
/* If n == 0, return numeric(0) */
PROTECT(x = allocVector(REALSXP, n));
if (n == 0)
{
UNPROTECT(1);
return(x);
}
/* Sanity checks of arguments. */
PROTECT(a = coerceVector(CADR(args), REALSXP));
PROTECT(b = coerceVector(CADDR(args), REALSXP));
bdims = getAttrib(b, R_DimSymbol);
nrow = INTEGER(bdims)[0];
ncol = INTEGER(bdims)[1];
if (nrow != ncol)
error(_("non-square sub-intensity matrix"));
na = LENGTH(a);
if (na != nrow)
error(_("non-conformable arguments"));
/* If length of parameters < 1, or either of the two parameters
* is NA return NA. */
if (na < 1 ||
(na == 1 && !(R_FINITE(REAL(a)[0]) && R_FINITE(REAL(b)[0]))))
{
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
}
/* Otherwise, dispatch to appropriate r* function */
else
{
naflag = FALSE;
GetRNGstate();
switch (code)
{
RANDPHTYPE2(1, rphtype);
default:
error(_("internal error in actuar_do_randomphtype2"));
}
if (naflag)
warning(R_MSG_NA);
PutRNGstate();
}
UNPROTECT(3);
return x;
}
/* Main function, the only one used by .External(). */
SEXP actuar_do_randomphtype(SEXP args)
{
int i;
const char *name;
/* Extract distribution name */
args = CDR(args);
name = CHAR(STRING_ELT(CAR(args), 0));
/* Dispatch to actuar_do_random{1,2,3,4} */
for (i = 0; random_tab[i].name; i++)
if (!strcmp(random_tab[i].name, name))
return random_tab[i].cfun(random_tab[i].code, CDR(args), random_tab[i].type);
/* No dispatch is an error */
error(_("internal error in actuar_do_randomphtype"));
return args; /* never used; to keep -Wall happy */
}
|