CppBUGS

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Introduction

Date:	January 20, 2012
Authors:	Whit Armstrong
Contact:	armstrong.whit@gmail.com
Web site:	http://github.com/armstrtw/CppBugs
License:	GPL-3

Purpose

CppBugs is a c++ library designed for MCMC sampling.

Features

CppBugs attempts to make writing mcmc models as painless as possible. It incorporates features from both WinBugs and PyMC and requires users only to implment an update method which resembles the model section of a WinBUGS script.

• CppBugs is extremely fast. Typically between 5x to 10x faster than equivalent WinBugs and 3x to 5x faster than PyMC models.
• Common statistical distributions are supported drawing heavily on Boost libraries. Many more will be implemented to eventually be as feature complete as WinBugs/PyMC.
  

Usage

1. vec b(randn<vec>(4));
   
2. vec b_herd(randn<vec>(N_herd));
   
3. vec overdisp(randn<vec>(N));
   
4. vec phi;
   
5. double tau_overdisp(1), tau_b_herd(1), sigma_overdisp(1), sigma_b_herd(1);
   
6. 
   
7. 
   
8. std::function<void ()> model = [&]() {
   
9.   phi = fixed*b + indicator_matrix*b_herd + overdisp;
   
10.   phi = 1/(1+exp(-phi));
    
11.   sigma_overdisp = 1/sqrt(tau_overdisp);
    
12.   sigma_b_herd = 1/sqrt(tau_b_herd);
    
13. };
    
14. 
    
15. MCModel m(model);
    
16. m.normal(b).dnorm(0,0.001);
    
17. m.uniform(tau_overdisp).dunif(0,1000);
    
18. m.uniform(tau_b_herd).dunif(0,100);
    
19. m.normal(b_herd).dnorm(0, tau_b_herd);
    
20. m.normal(overdisp).dnorm(0,tau_overdisp);
    
21. m.binomial(incidence).dbinom(size,phi);
    
22. m.deterministic(sigma_overdisp);
    
23. m.deterministic(sigma_b_herd);
    
24. m.deterministic(phi);
    
25. That's it. The model can be compiled and run as follows:
    
26. 
    
27. #include <iostream>
    
28. #include <vector>
    
29. #include <armadillo>
    
30. #include <boost/random.hpp>
    
31. #include <cppbugs/cppbugs.hpp>
    
32. #include <cppbugs/mcmc.model.hpp>
    
33. 
    
34. using namespace arma;
    
35. using namespace cppbugs;
    
36. using std::cout;
    
37. using std::endl;
    
38. 
    
39. int main() {
    
40.   int incidence_raw[] = {2,3,4,0,3,1,1,8,2,0,2,2,0,2,0,5,0,0,1,3,0,0,1,8,1,3,0,12,2,0,0,0,1,1,0,2,0,5,3,1,2,1,0,0,1,2,0,0,11,0,0,0,1,1,1,0};
    
41.   int size_raw[] = {14,12,9,5,22,18,21,22,16,16,20,10,10,9,6,18,25,24,4,17,17,18,20,16,10,9,5,34,9,6,8,6,22,22,18,22,25,27,22,22,10,8,6,5,21,24,19,23,19,2,3,2,19,15,15,15};
    
42.   int herd_raw[] = {1,1,1,1,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,9,9,9,9,10,10,10,10,11,11,11,11,12,12,12,12,13,13,13,13,14,14,14,14,15,15,15,15};
    
43.   double period2_raw[] = {0,1,0,0,0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0};
    
44.   double period3_raw[] = {0,0,1,0,0,0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0};
    
45.   double period4_raw[] = {0,0,0,1,0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1};
    
46. 
    
47.   int N = 56;
    
48.   int N_herd = 15;
    
49. 
    
50.   const ivec incidence(incidence_raw,N);
    
51.   const ivec size(size_raw,N);
    
52.   ivec herd(herd_raw,N); herd -= 1;
    
53.   const vec period2(period2_raw,N);
    
54.   const vec period3(period3_raw,N);
    
55.   const vec period4(period4_raw,N);
    
56. 
    
57.   mat indicator_matrix(N,N_herd);
    
58.   indicator_matrix.fill(0.0);
    
59.   for(uint i = 0; i < herd.n_elem; i++) {
    
60.     indicator_matrix(i,herd[i]) = 1.0;
    
61.   }
    
62. 
    
63.   mat fixed(N,4);
    
64.   fixed.col(0).fill(1);
    
65.   fixed.col(1) = period2;
    
66.   fixed.col(2) = period3;
    
67.   fixed.col(3) = period4;
    
68. 
    
69.   vec b(randn<vec>(4));
    
70.   vec b_herd(randn<vec>(N_herd));
    
71.   vec overdisp(randn<vec>(N));
    
72.   vec phi;
    
73.   double tau_overdisp(1), tau_b_herd(1), sigma_overdisp(1), sigma_b_herd(1);
    
74. 
    
75.   std::function<void ()> model = [&]() {
    
76.     phi = fixed*b + indicator_matrix*b_herd + overdisp;
    
77.     phi = 1/(1+exp(-phi));
    
78.     sigma_overdisp = 1/sqrt(tau_overdisp);
    
79.     sigma_b_herd = 1/sqrt(tau_b_herd);
    
80.   };
    
81. 
    
82.   MCModel m(model);
    
83.   m.normal(b).dnorm(0,0.001);
    
84.   m.uniform(tau_overdisp).dunif(0,1000);
    
85.   m.uniform(tau_b_herd).dunif(0,100);
    
86.   m.normal(b_herd).dnorm(0, tau_b_herd);
    
87.   m.normal(overdisp).dnorm(0,tau_overdisp);
    
88.   m.binomial(incidence).dbinom(size,phi);
    
89.   m.deterministic(sigma_overdisp);
    
90.   m.deterministic(sigma_b_herd);
    
91.   m.deterministic(phi);
    
92.   m.sample(1e6,1e5,1e4,50);
    
93. 
    
94.   cout << "samples: " << m.getNode(b).history.size() << endl;
    
95.   cout << "b: " << endl << m.getNode(b).mean() << endl;
    
96.   cout << "tau_overdisp: " << m.getNode(tau_overdisp).mean() << endl;
    
97.   cout << "tau_b_herd: " << m.getNode(tau_b_herd).mean() << endl;
    
98.   cout << "sigma_overdisp: " << m.getNode(sigma_overdisp).mean() << endl;
    
99.   cout << "sigma_b_herd: " << m.getNode(sigma_b_herd).mean() << endl;
    
100.   cout << "b_herd: " << endl << m.getNode(b_herd).mean() << endl;
     
101.   cout << "acceptance_ratio: " << m.acceptance_ratio() << endl;
     
102. 
     
103.   return 0;
     
104. }

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：BUGS bug introduction distribution Incorporate designed features possible library between