@@ -5,38 +5,44 @@ using namespace Rcpp;
5
5
6
6
#include < math.h>
7
7
8
- List markovchainFit (SEXP data, String method=" mle" , bool byrow=true , int nboot=10 , double laplacian=0
9
- , String name=" " , bool parallel=false , double confidencelevel=0.95 , bool confint = true
10
- , NumericMatrix hyperparam = NumericMatrix(), bool sanitize = false, CharacterVector possibleStates = CharacterVector());
8
+ List markovchainFit (SEXP data, String method = " mle" , bool byrow = true ,
9
+ int nboot = 10 , double laplacian = 0 , String name = " " ,
10
+ bool parallel = false , double confidencelevel = 0.95 , bool confint = true ,
11
+ NumericMatrix hyperparam = NumericMatrix(), bool sanitize = false,
12
+ CharacterVector possibleStates = CharacterVector());
11
13
12
14
// [[Rcpp::export]]
13
- List ctmcFit (List data, bool byrow=true , String name=" " , double confidencelevel = 0.95 )
14
- {
15
+ List ctmcFit (List data, bool byrow=true , String name=" " , double confidencelevel = 0.95 ) {
16
+
15
17
CharacterVector stateData (as<CharacterVector>(data[0 ]).size ());
16
- for (int i = 0 ; i < as<CharacterVector>(data[0 ]).size (); i++)
18
+
19
+ for (int i = 0 ; i < as<CharacterVector>(data[0 ]).size (); i++)
17
20
stateData[i] = as<CharacterVector>(data[0 ])[i];
21
+
18
22
NumericVector transData = data[1 ];
19
23
CharacterVector sortedStates = unique (as<CharacterVector>(data[0 ])).sort ();
20
24
NumericVector stateCount (sortedStates.size ());
21
25
NumericVector stateSojournTime (sortedStates.size ());
22
26
23
27
List dtmcData = markovchainFit (stateData, " mle" , byrow, 10 , 0 , name, false , confidencelevel);
24
28
25
- for (int i = 0 ; i < stateData.size () - 1 ; i++){
26
- int idx = std::find (sortedStates.begin (), sortedStates.end (), stateData[i]) - sortedStates.begin ();
29
+ for (int i = 0 ; i < stateData.size () - 1 ; i++){
30
+ int idx = std::find (sortedStates.begin (),
31
+ sortedStates.end (),
32
+ stateData[i]) - sortedStates.begin ();
27
33
stateCount[idx]++;
28
34
stateSojournTime[idx] += transData[i+1 ] - transData[i];
29
35
}
30
36
31
37
S4 dtmcEst = dtmcData[" estimate" ];
32
38
NumericMatrix gen = dtmcEst.slot (" transitionMatrix" );
33
39
34
- for (int i = 0 ; i < gen.nrow (); i++){
35
- for (int j = 0 ; j < gen.ncol (); j++){
36
- if (stateCount[i] > 0 )
40
+ for (int i = 0 ; i < gen.nrow (); i++){
41
+ for (int j = 0 ; j < gen.ncol (); j++){
42
+ if (stateCount[i] > 0 )
37
43
gen (i, j) *= stateCount[i] / stateSojournTime[i];
38
44
}
39
- if (stateCount[i] > 0 )
45
+ if (stateCount[i] > 0 )
40
46
gen (i, i) = - stateCount[i] / stateSojournTime[i];
41
47
else
42
48
gen (i, i) = -1 ;
@@ -45,12 +51,13 @@ List ctmcFit(List data, bool byrow=true, String name="", double confidencelevel
45
51
double zscore = stats::qnorm_0 (confidencelevel, 1.0 , 0.0 );
46
52
NumericVector lowerConfVecLambda (sortedStates.size ()), upperConfVecLambda (sortedStates.size ());
47
53
48
- for (int i = 0 ; i < sortedStates.size (); i++){
49
- if (stateCount[i] > 0 ){
50
- lowerConfVecLambda (i) = std::max (0 ., stateCount[i] / stateSojournTime[i] * (1 - zscore / sqrt (stateCount[i])));
51
- upperConfVecLambda (i) = std::min (1 ., stateCount[i] / stateSojournTime[i] * (1 + zscore / sqrt (stateCount[i])));
52
- }
53
- else {
54
+ for (int i = 0 ; i < sortedStates.size (); i++){
55
+
56
+ if (stateCount[i] > 0 ){
57
+ auto factor = stateCount[i] / stateSojournTime[i] * (1 - zscore / sqrt (stateCount[i]));
58
+ lowerConfVecLambda (i) = std::max (0 ., factor);
59
+ upperConfVecLambda (i) = std::min (1 ., factor);
60
+ } else {
54
61
lowerConfVecLambda (i) = 1 ;
55
62
upperConfVecLambda (i) = 1 ;
56
63
}
@@ -62,10 +69,13 @@ List ctmcFit(List data, bool byrow=true, String name="", double confidencelevel
62
69
outCtmc.slot (" name" ) = name;
63
70
64
71
return List::create (_[" estimate" ] = outCtmc,
65
- _[" errors" ] = List::create (_[" dtmcConfidenceInterval" ] = List::create (
66
- _[" confidenceLevel" ] = dtmcData[" confidenceLevel" ],
67
- _[" lowerEndpointMatrix" ] = dtmcData[" lowerEndpointMatrix" ],
68
- _[" upperEndpointMatrix" ] = dtmcData[" upperEndpointMatrix" ]),
69
- _[" lambdaConfidenceInterval" ] = List::create (_[" lowerEndpointVector" ] = lowerConfVecLambda,
70
- _[" upperEndpointVector" ] = upperConfVecLambda)));
71
- }
72
+ _[" errors" ] = List::create (
73
+ _[" dtmcConfidenceInterval" ] = List::create (
74
+ _[" confidenceLevel" ] = dtmcData[" confidenceLevel" ],
75
+ _[" lowerEndpointMatrix" ] = dtmcData[" lowerEndpointMatrix" ],
76
+ _[" upperEndpointMatrix" ] = dtmcData[" upperEndpointMatrix" ]),
77
+ _[" lambdaConfidenceInterval" ] = List::create (
78
+ _[" lowerEndpointVector" ] = lowerConfVecLambda,
79
+ _[" upperEndpointVector" ] = upperConfVecLambda))
80
+ );
81
+ }
0 commit comments