Browse code

Updates C to RCpp, begins using Roxygen2

Max McGrath authored on 22/07/2021 15:30:05
Showing8 changed files

... ...
@@ -1,4 +1,5 @@
1 1
 Package: discordant
2
+Type: Package
2 3
 Version: 1.17.0
3 4
 Date: 2016-10-21
4 5
 Title: The Discordant Method: A Novel Approach for Differential
... ...
@@ -17,12 +18,22 @@ Maintainer: Max McGrath <max.mcgrath@ucdenver.edu>
17 18
 Description: Discordant is a method to determine differential
18 19
         correlation of molecular feature pairs from -omics data using
19 20
         mixture models. Algorithm is explained further in Siska et al.
20
-Encoding: latin1
21
-biocViews: ImmunoOncology, BiologicalQuestion, StatisticalMethod, mRNAMicroarray,
22
-        Microarray, Genetics, RNASeq
23
-Suggests: BiocStyle, knitr
24
-Imports: Biobase, stats, biwt, gtools, MASS, tools
21
+Encoding: UTF-8
22
+biocViews: ImmunoOncology, BiologicalQuestion, StatisticalMethod, 
23
+    mRNAMicroarray,
24
+    Microarray, Genetics, RNASeq
25
+Suggests: BiocStyle, knitr, testthat
26
+Imports: 
27
+    Rcpp,
28
+    Biobase,
29
+    stats,
30
+    biwt,
31
+    gtools,
32
+    MASS,
33
+    tools
25 34
 License: GPL (>= 2)
26 35
 URL: https://github.com/siskac/discordant
27 36
 NeedsCompilation: yes
37
+LinkingTo: Rcpp
28 38
 VignetteBuilder: knitr
39
+RoxygenNote: 7.1.1
... ...
@@ -1,10 +1,14 @@
1
-useDynLib(discordant, .registration = TRUE)
1
+# Generated by roxygen2: do not edit by hand
2 2
 
3
+export(createVectors)
4
+export(discordantRun)
5
+export(fishersTrans)
6
+export(splitMADOutlier)
7
+import(Biobase)
3 8
 import(MASS)
4
-import(gtools)
5 9
 import(biwt)
6
-import(tools)
10
+import(gtools)
7 11
 import(stats)
8
-import(Biobase)
9
-
10
-export(createVectors, discordantRun, fishersTrans, splitMADOutlier)
12
+import(tools)
13
+importFrom(Rcpp,sourceCpp)
14
+useDynLib(discordant, .registration = TRUE)
11 15
new file mode 100644
... ...
@@ -0,0 +1,11 @@
1
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
2
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
3
+
4
+em_normal_partial_concordant_cpp <- function(x, y, zxy, n, pi, mu, sigma, nu, tau, g, loglik, tol, restriction, constrain, iteration, convergence) {
5
+    .Call(`_discordant_em_normal_partial_concordant_cpp`, x, y, zxy, n, pi, mu, sigma, nu, tau, g, loglik, tol, restriction, constrain, iteration, convergence)
6
+}
7
+
8
+subsampling_cpp <- function(x, y, zxy, n, pi, mu, sigma, nu, tau, g) {
9
+    .Call(`_discordant_subsampling_cpp`, x, y, zxy, n, pi, mu, sigma, nu, tau, g)
10
+}
11
+
0 12
new file mode 100644
... ...
@@ -0,0 +1,5 @@
1
+## usethis namespace: start
2
+#' @useDynLib discordant, .registration = TRUE
3
+#' @importFrom Rcpp sourceCpp
4
+## usethis namespace: end
5
+NULL
0 6
\ No newline at end of file
... ...
@@ -158,9 +158,11 @@ SparCC.frac <- function(x, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
158 158
 
159 159
 
160 160
 #modified from package mclust
161
-unmap <- function(classification){
161
+unmap <- function(classification, components){
162 162
     n <- length(classification)
163
-    u <- sort(unique(classification))
163
+    # u <- sort(unique(classification)) # OG Code
164
+    # u <- 0:max(classification) # Max's potential fix
165
+    u <- 0:(components - 1)
164 166
     labs <- as.character(u)
165 167
     k <- length(u)
166 168
     z <- matrix(0, n, k)
... ...
@@ -171,7 +173,7 @@ unmap <- function(classification){
171 173
 
172 174
 
173 175
 em.normal.partial.concordant <- function(data, class, tol=0.001, 
174
-    restriction=0, constrain=0, iteration=1000){
176
+    restriction=0, constrain=0, iteration=1000, components){
175 177
     n <- as.integer(dim(data)[1])
176 178
     g <- as.integer(nlevels(as.factor(class)))
177 179
 
... ...
@@ -182,8 +184,8 @@ em.normal.partial.concordant <- function(data, class, tol=0.001,
182 184
         return( c(diag(z[k,])) )
183 185
     }
184 186
 
185
-    zx <- unmap(class[,1])
186
-    zy <- unmap(class[,2])
187
+    zx <- unmap(class[,1], components = components)
188
+    zy <- unmap(class[,2], components = components)
187 189
     zxy <- sapply(1:dim(zx)[1], yl.outer, zx, zy)
188 190
 
189 191
     pi <- double(g*g)
... ...
@@ -193,10 +195,21 @@ em.normal.partial.concordant <- function(data, class, tol=0.001,
193 195
     tau <- double(g)
194 196
     loglik <- double(1)
195 197
     convergence <- integer(1)
196
-    results <- .C("em_normal_partial_concordant", as.double(data[,1]), 
197
-        as.double(data[,2]), as.double(t(zxy)), n, pi, mu, sigma, nu, tau, g, 
198
-        loglik, as.double(tol), as.integer(restriction), 
199
-        as.integer(constrain), as.integer(iteration), convergence)
198
+    # results <- .C("em_normal_partial_concordant", as.double(data[,1]), 
199
+    #     as.double(data[,2]), as.double(t(zxy)), n, pi, mu, sigma, nu, tau, g, 
200
+    #     loglik, as.double(tol), as.integer(restriction), 
201
+    #     as.integer(constrain), as.integer(iteration), convergence)
202
+    
203
+    results <- em_normal_partial_concordant_cpp(as.double(data[,1]), 
204
+                                                as.double(data[,2]), 
205
+                                                as.double(t(zxy)), n, pi, mu, 
206
+                                                sigma, nu, tau, g, loglik, 
207
+                                                as.double(tol), 
208
+                                                as.integer(restriction), 
209
+                                                as.integer(constrain), 
210
+                                                as.integer(iteration), 
211
+                                                convergence)
212
+    
200 213
     return(list(model="PCD", convergence=results[[16]], 
201 214
         pi=t(array(results[[5]], dim=c(g,g))), mu_sigma=rbind(results[[6]], 
202 215
         results[[7]]), nu_tau=rbind(results[[8]], results[[9]]), 
... ...
@@ -204,7 +217,7 @@ em.normal.partial.concordant <- function(data, class, tol=0.001,
204 217
         order,decreasing=TRUE)[1,], z=array(results[[3]], dim=c(n,g*g))))
205 218
 }
206 219
 
207
-subSampleData <- function(pdata, class, mu, sigma, nu, tau, pi) {
220
+subSampleData <- function(pdata, class, mu, sigma, nu, tau, pi, components) {
208 221
     n <- as.integer(dim(pdata)[1])
209 222
     g <- as.integer(nlevels(as.factor(class)))
210 223
 
... ...
@@ -215,19 +228,26 @@ subSampleData <- function(pdata, class, mu, sigma, nu, tau, pi) {
215 228
         return( c(diag(z[k,])) )
216 229
     }
217 230
 
218
-    zx <- unmap(class[,1])
219
-    zy <- unmap(class[,2])
231
+    zx <- unmap(class[,1], components = components)
232
+    zy <- unmap(class[,2], components = components)
220 233
     zxy <- sapply(1:dim(zx)[1], yl.outer, zx, zy)
221 234
 
222
-    results <- .C("subsampling", as.double(pdata[,1]), as.double(pdata[,2]), 
223
-        as.double(t(zxy)), n, as.double(pi), as.double(mu), as.double(sigma), 
224
-        as.double(nu), as.double(tau), g)
235
+    # results <- .C("subsampling", as.double(pdata[,1]), as.double(pdata[,2]), 
236
+    #     as.double(t(zxy)), n, as.double(pi), as.double(mu), as.double(sigma), 
237
+    #     as.double(nu), as.double(tau), g)
238
+    
239
+    results <- subsampling_cpp(as.double(pdata[,1]), as.double(pdata[,2]), 
240
+                               as.double(t(zxy)), n, as.double(pi), 
241
+                               as.double(mu), as.double(sigma), as.double(nu), 
242
+                               as.double(tau), g)
243
+    
225 244
     return(list(pi=t(array(results[[5]],dim=c(g,g))), 
226 245
         mu_sigma=rbind(results[[6]], results[[7]]), nu_tau=rbind(results[[8]],
227 246
         results[[9]]), class=apply(array(results[[3]], dim=c(n,g*g)),1,order,
228 247
         decreasing=TRUE)[1,], z=array(results[[3]], dim=c(n,g*g))))
229 248
 } 
230 249
 
250
+#' @export
231 251
 fishersTrans <- function(rho) {
232 252
     r = (1+rho)/(1-rho)
233 253
     z = 0.5*log(r,base = exp(1))
... ...
@@ -280,6 +300,14 @@ checkInputs <- function(x,y,groups = NULL) {
280 300
     return(issue)
281 301
 }
282 302
 
303
+#' @import Biobase
304
+#' @import biwt
305
+#' @import gtools
306
+#' @import MASS
307
+#' @import stats
308
+#' @import tools
309
+#' 
310
+#' @export
283 311
 createVectors <- function(x, y = NULL, groups, cor.method = c("spearman")) {
284 312
     print(x)
285 313
     if(checkInputs(x,y,groups)) {
... ...
@@ -362,6 +390,7 @@ createVectors <- function(x, y = NULL, groups, cor.method = c("spearman")) {
362 390
     return(list(v1 = statVector1, v2 = statVector2))
363 391
 }
364 392
 
393
+#' @export
365 394
 discordantRun <- function(v1, v2, x, y = NULL, transform = TRUE, 
366 395
     subsampling = FALSE, subSize = dim(x)[1], iter = 100, components = 3) {
367 396
 
... ...
@@ -458,7 +487,7 @@ discordantRun <- function(v1, v2, x, y = NULL, transform = TRUE,
458 487
             }
459 488
             pd <- em.normal.partial.concordant(sub.pdata, sub.class, 
460 489
                 tol=0.001, restriction=0, constrain=c(0,-sd(pdata),sd(pdata)),
461
-                iteration=1000)
490
+                iteration=1000, components=components)
462 491
             total_mu <- total_mu + pd$mu_sigma[1,]
463 492
             total_sigma <- total_sigma + pd$mu_sigma[2,]
464 493
             total_nu <- total_nu + pd$nu_tau[1,]
... ...
@@ -472,13 +501,13 @@ discordantRun <- function(v1, v2, x, y = NULL, transform = TRUE,
472 501
         tau <- total_tau/iter
473 502
         pi <- total_pi/iter
474 503
 
475
-        finalResult <- subSampleData(pdata, class, mu, sigma, nu, tau, pi)
504
+        finalResult <- subSampleData(pdata, class, mu, sigma, nu, tau, pi, components)
476 505
         zTable <- finalResult$z
477 506
         classVector <- finalResult$class
478 507
     } else {
479 508
         pd <- em.normal.partial.concordant(pdata, class, tol=0.001, 
480 509
             restriction=0, constrain=c(0,-sd(pdata),sd(pdata)), 
481
-            iteration=1000)
510
+            iteration=1000, components = components)
482 511
         zTable <- pd$z
483 512
         classVector <- pd$class
484 513
     }
... ...
@@ -522,6 +551,7 @@ discordantRun <- function(v1, v2, x, y = NULL, transform = TRUE,
522 551
         probMatrix = zTable, loglik = pd$loglik))
523 552
 }
524 553
 
554
+#' @export
525 555
 splitMADOutlier <- function(mat, filter0 = TRUE, threshold = 2) {
526 556
     if(mode(mat) != "S4") {
527 557
         stop("data matrix mat must be type ExpressionSet")
528 558
new file mode 100644
... ...
@@ -0,0 +1,64 @@
1
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
2
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
3
+
4
+#include <Rcpp.h>
5
+
6
+using namespace Rcpp;
7
+
8
+// em_normal_partial_concordant_cpp
9
+Rcpp::List em_normal_partial_concordant_cpp(Rcpp::NumericVector x, Rcpp::NumericVector y, Rcpp::NumericVector zxy, int n, Rcpp::NumericVector pi, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector nu, Rcpp::NumericVector tau, int g, double loglik, double tol, int restriction, Rcpp::NumericVector constrain, int iteration, int convergence);
10
+RcppExport SEXP _discordant_em_normal_partial_concordant_cpp(SEXP xSEXP, SEXP ySEXP, SEXP zxySEXP, SEXP nSEXP, SEXP piSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP nuSEXP, SEXP tauSEXP, SEXP gSEXP, SEXP loglikSEXP, SEXP tolSEXP, SEXP restrictionSEXP, SEXP constrainSEXP, SEXP iterationSEXP, SEXP convergenceSEXP) {
11
+BEGIN_RCPP
12
+    Rcpp::RObject rcpp_result_gen;
13
+    Rcpp::RNGScope rcpp_rngScope_gen;
14
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP);
15
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type y(ySEXP);
16
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type zxy(zxySEXP);
17
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
18
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type pi(piSEXP);
19
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type mu(muSEXP);
20
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type sigma(sigmaSEXP);
21
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type nu(nuSEXP);
22
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type tau(tauSEXP);
23
+    Rcpp::traits::input_parameter< int >::type g(gSEXP);
24
+    Rcpp::traits::input_parameter< double >::type loglik(loglikSEXP);
25
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
26
+    Rcpp::traits::input_parameter< int >::type restriction(restrictionSEXP);
27
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type constrain(constrainSEXP);
28
+    Rcpp::traits::input_parameter< int >::type iteration(iterationSEXP);
29
+    Rcpp::traits::input_parameter< int >::type convergence(convergenceSEXP);
30
+    rcpp_result_gen = Rcpp::wrap(em_normal_partial_concordant_cpp(x, y, zxy, n, pi, mu, sigma, nu, tau, g, loglik, tol, restriction, constrain, iteration, convergence));
31
+    return rcpp_result_gen;
32
+END_RCPP
33
+}
34
+// subsampling_cpp
35
+Rcpp::List subsampling_cpp(Rcpp::NumericVector x, Rcpp::NumericVector y, Rcpp::NumericVector zxy, int n, Rcpp::NumericVector pi, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector nu, Rcpp::NumericVector tau, int g);
36
+RcppExport SEXP _discordant_subsampling_cpp(SEXP xSEXP, SEXP ySEXP, SEXP zxySEXP, SEXP nSEXP, SEXP piSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP nuSEXP, SEXP tauSEXP, SEXP gSEXP) {
37
+BEGIN_RCPP
38
+    Rcpp::RObject rcpp_result_gen;
39
+    Rcpp::RNGScope rcpp_rngScope_gen;
40
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP);
41
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type y(ySEXP);
42
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type zxy(zxySEXP);
43
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
44
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type pi(piSEXP);
45
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type mu(muSEXP);
46
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type sigma(sigmaSEXP);
47
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type nu(nuSEXP);
48
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type tau(tauSEXP);
49
+    Rcpp::traits::input_parameter< int >::type g(gSEXP);
50
+    rcpp_result_gen = Rcpp::wrap(subsampling_cpp(x, y, zxy, n, pi, mu, sigma, nu, tau, g));
51
+    return rcpp_result_gen;
52
+END_RCPP
53
+}
54
+
55
+static const R_CallMethodDef CallEntries[] = {
56
+    {"_discordant_em_normal_partial_concordant_cpp", (DL_FUNC) &_discordant_em_normal_partial_concordant_cpp, 16},
57
+    {"_discordant_subsampling_cpp", (DL_FUNC) &_discordant_subsampling_cpp, 10},
58
+    {NULL, NULL, 0}
59
+};
60
+
61
+RcppExport void R_init_discordant(DllInfo *dll) {
62
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
63
+    R_useDynamicSymbols(dll, FALSE);
64
+}
0 65
deleted file mode 100644
... ...
@@ -1,269 +0,0 @@
1
-#include <R.h>
2
-#include <R_ext/Rdynload.h>
3
-#include <Rinternals.h>
4
-
5
-void em_normal_partial_concordant(double *x, double *y, double *zxy, int *n, double *pi, double *mu, double *sigma, double *nu, double *tau, int *g, double *loglik, double *tol, int *restriction, int *constrain, int *iteration, int *convergence) {
6
-	int i, j, k, flag, iter;
7
-	double temp, loglik2;
8
-
9
-
10
-	/*theoretical restriction on the null*/
11
-	if((*restriction)>0){
12
-	/*initialization*/
13
-		flag=1;
14
-		(*loglik) = 0;
15
-		for(k=0;k<*n;k++){	
16
-			(*loglik) = (*loglik) - 0.5*x[k]*x[k] - 0.5*log(2*PI) - 0.5*y[k]*y[k] - 0.5*log(2*PI);
17
-		}
18
-		/*iteration*/
19
-		
20
-		iter = 0;
21
-		while(flag>0 & iter<(*iteration)){
22
-		/*M step for pi, mu and sigma*/
23
-			for(i=0;i<(*g);i++){
24
-				pi[i] = 0;
25
-				for(k=0;k<(*n);k++){
26
-					for(j=0;j<(*g);j++){
27
-						pi[i] = pi[i] + zxy[(j*(*g)+i)*(*n)+k];
28
-					}
29
-				}
30
-				mu[i] = 0;
31
-				for(k=0;k<(*n);k++){
32
-					for(j=0;j<(*g);j++){
33
-						mu[i] = mu[i] + zxy[(j*(*g)+i)*(*n)+k]*x[k];
34
-					}
35
-				}
36
-				mu[i] = mu[i] / pi[i];
37
-				sigma[i] = 0;
38
-				for(k=0;k<(*n);k++){
39
-					for(j=0;j<(*g);j++){
40
-						sigma[i] = sigma[i] + zxy[(j*(*g)+i)*(*n)+k]*(x[k]-mu[i])*(x[k]-mu[i]);
41
-					}
42
-				}
43
-				sigma[i] = sigma[i] / pi[i];
44
-			}
45
-			for(j=0;j<(*g);j++){
46
-				pi[j] = 0;
47
-				for(k=0;k<(*n);k++){
48
-					for(i=0;i<(*g);i++){
49
-						pi[j] = pi[j] + zxy[(j*(*g)+i)*(*n)+k];
50
-					}
51
-				}
52
-				nu[j] = 0;
53
-				for(k=0;k<(*n);k++){
54
-					for(i=0;i<(*g);i++){
55
-						nu[j] = nu[j] + zxy[(j*(*g)+i)*(*n)+k]*y[k];
56
-					}
57
-				}
58
-				nu[j] = nu[j] / pi[j];
59
-				tau[j] = 0;
60
-				for(k=0;k<(*n);k++){
61
-					for(i=0;i<(*g);i++){
62
-						tau[j] = tau[j] + zxy[(j*(*g)+i)*(*n)+k]*(y[k]-nu[j])*(y[k]-nu[j]);
63
-					}
64
-				}
65
-				tau[j] = tau[j] / pi[j];
66
-			}
67
-			for(i=0;i<(*g);i++){
68
-				for(j=0;j<(*g);j++){
69
-					pi[i*(*g)+j] = 0;
70
-					for(k=0;k<(*n);k++){
71
-						pi[i*(*g)+j] = pi[i*(*g)+j] + zxy[(j*(*g)+i)*(*n)+k];
72
-					}
73
-					pi[i*(*g)+j] = pi[i*(*g)+j] / (double)(*n);
74
-				}
75
-			}
76
-			/*constrain*/
77
-			for(i=0;i<(*g);i++){
78
-				if(constrain[i]==0-1){
79
-					if(mu[i]>0.0){
80
-						mu[i]=0.0;
81
-					}
82
-					if(nu[i]>0.0){
83
-						nu[i]=0.0;
84
-					}
85
-				}
86
-				if(constrain[i]==0){
87
-					mu[i] = 0.0;
88
-					sigma[i] = 1.0;
89
-					nu[i] = 0.0;
90
-					tau[i] = 1.0;
91
-				}
92
-				if(constrain[i]==0+1){
93
-					if(mu[i]<0.0){
94
-						mu[i]=0.0;
95
-					}
96
-					if(nu[i]<0.0){
97
-						nu[i]=0.0;
98
-					}
99
-				}
100
-			}
101
-			/*check convergence*/
102
-			loglik2 = (*loglik);
103
-			(*loglik) = 0;
104
-			for(k=0;k<(*n);k++){
105
-				temp = 0;
106
-				for(i=0;i<(*g);i++){
107
-					for(j=0;j<(*g);j++){
108
-						temp = temp + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) );
109
-					}
110
-				}
111
-				(*loglik) = (*loglik) + log(temp);
112
-			}
113
-			if(fabs((*loglik)-loglik2)<(*tol)){
114
-				flag = 0;
115
-			}
116
-			/*E step for z*/
117
-			for(k=0;k<(*n);k++){
118
-				temp = 0;
119
-				for(i=0;i<(*g);i++){
120
-					for(j=0;j<(*g);j++){
121
-						temp = temp + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) );
122
-					}
123
-				}
124
-				for(i=0;i<(*g);i++){
125
-					for(j=0;j<(*g);j++){
126
-						zxy[(j*(*g)+i)*(*n)+k] = pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) )/temp;
127
-					}
128
-				}
129
-			}
130
-			iter = iter + 1;
131
-		}
132
-		(*convergence) = 0;
133
-		if(iter<(*iteration)){
134
-			(*convergence) = 1;
135
-		}
136
-	}
137
-	/*no restriction on the null*/
138
-	else{
139
-	/*initialization*/
140
-		flag=1;
141
-		(*loglik) = 0;
142
-		for(k=0;k<(*n);k++){
143
-			(*loglik) = (*loglik) - 0.5*x[k]*x[k] - 0.5*log(2*PI) - 0.5*y[k]*y[k] - 0.5*log(2*PI);
144
-		}
145
-		/*iteration*/
146
-		iter = 0;
147
-		while(flag>0 & iter<(*iteration)){
148
-		/*M step for pi, mu and sigma*/
149
-			for(i=0;i<(*g);i++){
150
-				pi[i] = 0;
151
-				for(k=0;k<(*n);k++){
152
-					for(j=0;j<(*g);j++){
153
-						pi[i] = pi[i] + zxy[(j*(*g)+i)*(*n)+k];
154
-					}
155
-				}
156
-				mu[i] = 0;
157
-				for(k=0;k<(*n);k++){
158
-					for(j=0;j<(*g);j++){
159
-						mu[i] = mu[i] + zxy[(j*(*g)+i)*(*n)+k]*x[k];
160
-					}
161
-				}
162
-				mu[i] = mu[i] / pi[i];
163
-				sigma[i] = 0;
164
-				for(k=0;k<(*n);k++){
165
-					for(j=0;j<(*g);j++){
166
-						sigma[i] = sigma[i] + zxy[(j*(*g)+i)*(*n)+k]*(x[k]-mu[i])*(x[k]-mu[i]);
167
-					}
168
-				}
169
-				sigma[i] = sigma[i] / pi[i];
170
-			}
171
-			for(j=0;j<(*g);j++){
172
-				pi[j] = 0;
173
-				for(k=0;k<(*n);k++){
174
-					for(i=0;i<(*g);i++){
175
-						pi[j] = pi[j] + zxy[(j*(*g)+i)*(*n)+k];
176
-					}
177
-				}
178
-				nu[j] = 0;
179
-				for(k=0;k<(*n);k++){
180
-					for(i=0;i<(*g);i++){
181
-						nu[j] = nu[j] + zxy[(j*(*g)+i)*(*n)+k]*y[k];
182
-					}
183
-				}
184
-				nu[j] = nu[j] / pi[j];
185
-				tau[j] = 0;
186
-				for(k=0;k<(*n);k++){
187
-					for(i=0;i<(*g);i++){
188
-						tau[j] = tau[j] + zxy[(j*(*g)+i)*(*n)+k]*(y[k]-nu[j])*(y[k]-nu[j]);
189
-					}
190
-				}
191
-				tau[j] = tau[j] / pi[j];
192
-			}
193
-			for(i=0;i<(*g);i++){
194
-				for(j=0;j<(*g);j++){
195
-					pi[i*(*g)+j] = 0;
196
-					for(k=0;k<(*n);k++){
197
-						pi[i*(*g)+j] = pi[i*(*g)+j] + zxy[(j*(*g)+i)*(*n)+k];
198
-					}
199
-					pi[i*(*g)+j] = pi[i*(*g)+j] / (double)(*n);
200
-				}
201
-			}
202
-			/*check convergence*/
203
-			loglik2 = (*loglik);
204
-			(*loglik) = 0;
205
-			for(k=0;k<(*n);k++){
206
-				temp = 0;
207
-				for(i=0;i<(*g);i++){
208
-					for(j=0;j<(*g);j++){
209
-						temp = temp + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) );
210
-					}
211
-				}
212
-				(*loglik) = (*loglik) + log(temp);
213
-			}
214
-			if(fabs((*loglik)-loglik2)<(*tol)){
215
-				flag = 0;
216
-			}
217
-			/*E step for z*/
218
-			for(k=0;k<(*n);k++){
219
-				temp = 0;
220
-				for(i=0;i<(*g);i++){
221
-					for(j=0;j<(*g);j++){
222
-						temp = temp + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) );
223
-					}
224
-				}
225
-				for(i=0;i<(*g);i++) {
226
-					for(j=0;j<(*g);j++){
227
-						zxy[(j*(*g)+i)*(*n)+k] = zxy[(j*(*g)+i)*(*n)+k] + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) )/temp;
228
-					}
229
-				}
230
-			}
231
-			iter = iter + 1;
232
-		}
233
-		(*convergence) = 0;
234
-		if(iter<(*iteration)){
235
-			(*convergence) = 1;
236
-		}
237
-	}
238
-}
239
-
240
-void subsampling(double *x, double *y, double *zxy, int *n, double *pi, double *mu, double *sigma, double *nu, double *tau, int *g) {
241
-        int i, j, k;
242
-        double temp;
243
-        for(k=0;k<(*n);k++){
244
-                temp = 0;
245
-                for(i=0;i<(*g);i++){
246
-                        for(j=0;j<(*g);j++){
247
-                                temp = temp + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) );
248
-                        }
249
-                }
250
-                for(i=0;i<(*g);i++) {
251
-                        for(j=0;j<(*g);j++){
252
-                                zxy[(j*(*g)+i)*(*n)+k] = zxy[(j*(*g)+i)*(*n)+k] + pi[i*(*g)+j] * ( exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i]) )*( exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) )/temp;
253
-                        }
254
-                }
255
-        }
256
-
257
-}
258
-
259
-R_CMethodDef cMethods[]  = {
260
-  {"em_normal_partial_concordant", (DL_FUNC) &em_normal_partial_concordant, 16},
261
-  {"subsampling", (DL_FUNC) &subsampling, 11},
262
-  {NULL, NULL, 0}
263
-};
264
-
265
-void R_init_discordant(DllInfo *info) {
266
-  R_registerRoutines(info, cMethods, NULL, NULL, NULL);
267
-  R_useDynamicSymbols(info, FALSE);
268
-}
269
-
270 0
new file mode 100644
... ...
@@ -0,0 +1,190 @@
1
+#include <Rcpp.h>
2
+#include <math.h>
3
+#include <iostream>
4
+using namespace Rcpp;
5
+
6
+// [[Rcpp::export]]
7
+Rcpp::List em_normal_partial_concordant_cpp(Rcpp::NumericVector x,
8
+                                        Rcpp::NumericVector y,
9
+                                        Rcpp::NumericVector zxy,
10
+                                        int n,
11
+                                        Rcpp::NumericVector pi,
12
+                                        Rcpp::NumericVector mu,
13
+                                        Rcpp::NumericVector sigma,
14
+                                        Rcpp::NumericVector nu,
15
+                                        Rcpp::NumericVector tau,
16
+                                        int g,
17
+                                        double loglik,
18
+                                        double tol,
19
+                                        int restriction,
20
+                                        Rcpp::NumericVector constrain,
21
+                                        int iteration,
22
+                                        int convergence) {
23
+    int i, j, k, flag, iter;
24
+    double temp, loglik2;
25
+    
26
+    /*initialization*/
27
+    flag = 1;
28
+    loglik = 0;
29
+    for (k = 0; k < n; k++) {
30
+        loglik = loglik - 0.5*x[k]*x[k] - 0.5*log(2*PI) - 
31
+            0.5*y[k]*y[k] - 0.5*log(2*PI);
32
+    }
33
+    
34
+    /*iteration*/
35
+    iter = 0;
36
+    while (flag > 0 & iter < iteration) {
37
+
38
+        /*M step for pi, mu and sigma*/
39
+        for (i = 0; i < g; i++) {
40
+            pi[i] = 0;
41
+            for (k = 0; k < n; k++) {
42
+                for (j = 0; j < g; j++){
43
+                    pi[i] = pi[i] + zxy[(j*g+i)*n+k];
44
+                }
45
+            }
46
+            
47
+            mu[i] = 0;
48
+            for (k = 0; k < n; k++) {
49
+                for (j = 0; j < g; j++) {
50
+                    mu[i] = mu[i] + zxy[(j*g+i)*n+k] * x[k];
51
+                }
52
+            }
53
+            mu[i] = mu[i] / pi[i];
54
+
55
+            sigma[i] = 0;
56
+            for (k = 0; k < n; k++) {
57
+                for (j = 0; j < g; j++) {
58
+                    sigma[i] = sigma[i] + zxy[(j*g+i)*n+k] * 
59
+                        (x[k]-mu[i]) * (x[k]-mu[i]);
60
+                }
61
+            }
62
+            sigma[i] = sigma[i] / pi[i];
63
+        }
64
+
65
+        for (j = 0; j < g; j++) {
66
+            pi[j] = 0;
67
+            for (k = 0; k < n; k++) {
68
+                for (i = 0; i < g; i++) {
69
+                    pi[j] = pi[j] + zxy[(j*g+i)*n+k];
70
+                }
71
+            }
72
+            
73
+            nu[j] = 0;
74
+            for (k = 0; k < n; k++) {
75
+                for (i = 0; i < g; i++) {
76
+                    nu[j] = nu[j] + zxy[(j*g+i)*n+k] * y[k];
77
+                }
78
+            }
79
+            nu[j] = nu[j] / pi[j];
80
+
81
+            tau[j] = 0;
82
+            for (k = 0; k < n; k++) {
83
+                for (i = 0; i < g; i++){
84
+                    tau[j] = tau[j] + zxy[(j*g+i)*n+k] * 
85
+                        (y[k]-nu[j]) * (y[k]-nu[j]);
86
+                }
87
+            }
88
+            tau[j] = tau[j] / pi[j];
89
+        }
90
+
91
+        for (i = 0; i < g; i++) {
92
+            for (j = 0; j < g; j++) {
93
+                pi[i*g+j] = 0;
94
+                for (k = 0; k < n; k++) {
95
+                    pi[i*g+j] = pi[i*g+j] + zxy[(j*g+i)*n+k];
96
+                }
97
+                pi[i*g+j] = pi[i*g+j] / static_cast<double>(n);
98
+            }
99
+        }
100
+
101
+        /* check convergence */
102
+        loglik2 = loglik;
103
+        loglik = 0;
104
+        for (k = 0; k < n; k++) {
105
+            temp = 0;
106
+            for (i = 0; i < g; i++) {
107
+                for (j = 0; j < g; j++) {
108
+                    temp = temp + pi[i*g+j] * 
109
+                        (exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i])) * 
110
+                        (exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]));
111
+                }
112
+            }
113
+            loglik = loglik + log(temp);
114
+        }
115
+        
116
+        if (fabs(loglik-loglik2) < tol){
117
+            flag = 0;
118
+        }
119
+        
120
+        /*E step for z*/
121
+        for (k = 0; k < n; k++){
122
+            temp = 0;
123
+            for (i = 0; i < g; i++) {
124
+                for (j = 0; j < g; j++) {
125
+                    temp = temp + pi[i*g+j] * 
126
+                        (exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i])) *
127
+                        (exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]));
128
+                }
129
+            }
130
+            for (i = 0; i < g;i++) {
131
+                for (j = 0; j < g; j++) {
132
+                    zxy[(j*g+i)*n+k] = zxy[(j*g+i)*n+k] + pi[i*g+j] * 
133
+                        (exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i])) *
134
+                        (exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j])) / temp;
135
+                }
136
+            }
137
+        }
138
+        iter = iter + 1;
139
+    
140
+    }
141
+
142
+    convergence = 0;
143
+    if(iter < iteration) {
144
+        convergence = 1;
145
+    }
146
+    
147
+    Rcpp::List rtn = Rcpp::List::create(x, y, zxy, n, pi, mu, sigma, nu,
148
+                                        tau, g, loglik, tol, restriction,
149
+                                        constrain, iteration, convergence);
150
+    
151
+    return rtn;
152
+}
153
+
154
+// [[Rcpp::export]]
155
+Rcpp::List subsampling_cpp(Rcpp::NumericVector x, 
156
+                       Rcpp::NumericVector y,
157
+                       Rcpp::NumericVector zxy,
158
+                       int n, 
159
+                       Rcpp::NumericVector pi, 
160
+                       Rcpp::NumericVector mu, 
161
+                       Rcpp::NumericVector sigma, 
162
+                       Rcpp::NumericVector nu, 
163
+                       Rcpp::NumericVector tau, 
164
+                       int g) {
165
+    int i, j, k;
166
+    double temp;
167
+    
168
+    for (k = 0; k < n; k++) {
169
+        temp = 0;
170
+        for (i = 0; i < g; i++) {
171
+            for (j = 0; j < g; j++) {
172
+                temp = temp + pi[i*g+j] * 
173
+                    (exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i])) *
174
+                    (exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]));
175
+            }
176
+        }
177
+        for (i = 0; i < g; i++) {
178
+            for (j = 0; j < g; j++) {
179
+                zxy[(j*g+i)*n+k] = zxy[(j*g+i)*n+k] + pi[i*g+j] * 
180
+                    (exp(0-0.5*(x[k]-mu[i])*(x[k]-mu[i])/sigma[i])/sqrt(2*PI*sigma[i])) *
181
+                    (exp(0-0.5*(y[k]-nu[j])*(y[k]-nu[j])/tau[j])/sqrt(2*PI*tau[j]) ) / temp;
182
+            }
183
+        }
184
+    }
185
+    
186
+    Rcpp::List rtn = Rcpp::List::create(x, y, zxy, n, pi, mu, 
187
+                                        sigma, nu, tau, g);
188
+    
189
+    return rtn;
190
+}