1
#include <Rcpp.h>
2

3
using namespace Rcpp;
4

5
// [[Rcpp::export]]
6 1
List mcmc_rank(IntegerMatrix P, 
7
               IntegerVector init_rank,
8
               int rp) {
9 1
  int t = 1;
10 1
  int n = P.rows();
11 1
  NumericMatrix rrp(n,n);
12 1
  NumericVector expected(n);
13 1
  IntegerVector elements=seq(0,n-1);
14

15
  //rrp
16 1
  for(int i=0; i<(n-1); ++i){
17 1
    int j=init_rank[i];
18 1
    IntegerVector tmp=init_rank[elements>i];
19 1
    for(int k=0;k<tmp.length();++k){
20 1
      rrp(j,tmp[k])=1;
21
    }
22
  }
23
  //expected
24 1
  for(int i=0;i<n; ++i){
25 1
      expected[init_rank[i]] =i;
26
  }
27
  
28
  //MC
29 1
  while(t<=rp){
30 1
    int p = floor(R::runif(0,1)*(n-1));//rand() % (n-1);
31 1
    int c = round(R::runif(0,1));//rand() % 2;
32 1
    int a = init_rank[p];
33 1
    int b = init_rank[p+1];
34 1
    t+=1;
35 1
    if((c==1) & (P(a,b)!=1)){
36 0
      init_rank[p]=b;
37 0
      init_rank[p+1]=a;
38
      //expected update
39 0
      for(int i=0;i<n; ++i){
40 0
        expected[init_rank[i]]=double(expected[init_rank[i]]*(t-1)+i)/double(t);
41
      }
42
      //rrp update
43 0
      for(int i=0; i<(n-1); ++i){
44 0
        int j=init_rank[i];
45 0
        IntegerVector tmp=init_rank[elements>i];
46 0
        for(int k=0;k<tmp.length();++k){
47 0
          rrp(j,tmp[k])=double(rrp(j,tmp[k])*(t-1)+1)/double(t);
48 0
          rrp(tmp[k],j)=1-rrp(j,tmp[k]);
49
        }
50
      }
51
    }
52
  }
53
  
54 1
  return Rcpp::List::create(Rcpp::Named("expected") = expected,
55 1
                            Rcpp::Named("rrp")=rrp);
56
}
57

Read our documentation on viewing source code .

Loading