-
Notifications
You must be signed in to change notification settings - Fork 0
/
loglikC.cpp
92 lines (87 loc) · 2.82 KB
/
loglikC.cpp
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
#include <Rcpp.h>
using namespace Rcpp;
// parameterization type C
// parm[0] is log(-log(S[0])) corresponding to first visit time
// parm[j] is log(-log(S[j])) - log(-log(S[j-1])) corresponding to change
// in log(-log(S)) in that time period
// This is inspired by representing survival function as exp(-exp(lambda + beta*Z))
// [[Rcpp::export]]
double loglikC0(NumericVector parm, NumericMatrix Dm) {
int nsub = Dm.nrow(), J = Dm.ncol() - 1, i, j;
double result = 0, temp, templik;
for (i = 0; i < nsub; i++) {
templik = Dm(i, 0);
temp = 0;
for (j = 0; j < J; j++) {
temp += parm[j];
templik += Dm(i, j+1)*exp(-exp(temp));
}
result += log(templik);
}
return -result;
}
// [[Rcpp::export]]
NumericVector gradlikC0(NumericVector parm, NumericMatrix Dm) {
int nsub = Dm.nrow(), J = Dm.ncol() - 1, i, j;
NumericVector result(J), Deriv(J);
double temp, templik, likj;
for (i = 0; i < nsub; i++) {
templik = Dm(i, 0);
temp = 0;
Deriv.fill(0);
for (j = 0; j < J; j++) {
temp += parm[j];
likj = Dm(i, j+1)*exp(-exp(temp));
templik += likj;
for (int k = 0; k <= j; k++) Deriv[k] += likj*exp(temp);
}
for (j = 0; j < J; j++) result[j] += Deriv[j]/templik;
}
return result;
}
// [[Rcpp::export]]
double loglikC(NumericVector parm, NumericMatrix Dm, NumericMatrix Xmat) {
int nsub = Dm.nrow(), J = Dm.ncol() - 1, nbeta = Xmat.ncol(), i, j, k;
double result = 0, temp, templik, b;
NumericVector lamb(J), beta(nbeta);
for (i = 0; i < J; i++) lamb[i] = parm[i];
for (i = 0; i < nbeta; i++) beta[i] = parm[J + i];
for (i = 0; i < nsub; i++) {
templik = Dm(i, 0);
b = 0;
for (k = 0; k < nbeta; k++) b += Xmat(i, k)*beta[k];
temp = b;
for (j = 0; j < J; j++) {
temp += parm[j];
templik += Dm(i, j+1)*exp(-exp(temp));
}
result += log(templik);
}
return -result;
}
// [[Rcpp::export]]
NumericVector gradlikC(NumericVector parm, NumericMatrix Dm, NumericMatrix Xmat) {
int nsub = Dm.nrow(), J = Dm.ncol() - 1, nbeta = Xmat.ncol(), i, j, k;
double temp, templik, b, likj;
NumericVector lamb(J), beta(nbeta), Dlamb(J), Dbeta(nbeta), result(J + nbeta);
for (i = 0; i < J; i++) lamb[i] = parm[i];
for (i = 0; i < nbeta; i++) beta[i] = parm[J + i];
for (i = 0; i < nsub; i++) {
templik = Dm(i, 0);
b = 0;
for (k = 0; k < nbeta; k++) b += Xmat(i, k)*beta[k];
temp = b;
Dlamb.fill(0);
Dbeta.fill(0);
for (j = 0; j < J; j++) {
temp += parm[j];
likj = Dm(i, j+1)*exp(-exp(temp));
templik += likj;
for (k = 0; k <= j; k++) Dlamb[k] += likj*exp(temp);
for (k = 0; k < nbeta; k++) Dbeta[k] += likj*exp(temp)*Xmat(i, k);
}
for (j = 0; j < J; j++) result[j] += Dlamb[j]/templik;
for (j = 0; j < nbeta; j++) result[J + j] += Dbeta[j]/templik;
}
return result;
}