forked from r-lib/cpp11
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmatrix.cpp
More file actions
87 lines (75 loc) · 1.93 KB
/
matrix.cpp
File metadata and controls
87 lines (75 loc) · 1.93 KB
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
#include "cpp11/matrix.hpp"
#include "Rmath.h"
#include "cpp11/doubles.hpp"
using namespace cpp11;
[[cpp11::register]] SEXP gibbs_cpp(int N, int thin) {
cpp11::writable::doubles_matrix mat(N, 2);
double x = 0, y = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < thin; j++) {
x = Rf_rgamma(3., 1. / double(y * y + 4));
y = Rf_rnorm(1. / (x + 1.), 1. / sqrt(2. * (x + 1.)));
// REprintf("x: %f y: %f\n", x, y);
}
mat[i][0] = x;
mat[i][1] = y;
}
return mat;
}
[[cpp11::register]] cpp11::doubles_matrix gibbs_cpp2(int N, int thin) {
cpp11::writable::doubles_matrix mat(N, 2);
double x = 0, y = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < thin; j++) {
x = Rf_rgamma(3., 1. / double(y * y + 4));
y = Rf_rnorm(1. / (x + 1.), 1. / sqrt(2. * (x + 1.)));
}
mat(i, 0) = x;
mat(i, 1) = y;
}
return mat;
}
#include <Rcpp.h>
using namespace Rcpp;
[[cpp11::register]] NumericMatrix gibbs_rcpp(int N, int thin) {
NumericMatrix mat(N, 2);
double x = 0, y = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < thin; j++) {
x = rgamma(1, 3, 1 / (y * y + 4))[0];
y = rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))[0];
}
mat(i, 0) = x;
mat(i, 1) = y;
}
return (mat);
}
[[cpp11::register]] NumericMatrix gibbs_rcpp2(int N, int thin) {
NumericMatrix mat(N, 2);
double x = 0, y = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < thin; j++) {
x = Rf_rgamma(3., 1. / (y * y + 4));
y = Rf_rnorm(1. / (x + 1), 1 / sqrt(2 * (x + 1)));
}
mat(i, 0) = x;
mat(i, 1) = y;
}
return (mat);
}
[[cpp11::register]] cpp11::doubles row_sums(cpp11::doubles_matrix x) {
cpp11::writable::doubles sums(x.nrow());
int i = 0;
for (auto& row : x) {
sums[i] = 0.;
for (auto&& val : row) {
if (cpp11::is_na(val)) {
sums[i] = NA_REAL;
break;
}
sums[i] += val;
}
++i;
}
return sums;
}