API tools faq
paste
Guest User

Untitled

a guest
Sep 29th, 2016
57
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.13 KB | None | 0 0
raw download clone embed print report
  1. #include <RcppArmadillo.h>
  2. #include <algorithm>
  3. #include <numeric>
  4. using namespace Rcpp;
  5. using namespace arma;
  6. // [[Rcpp::depends(RcppArmadillo)]]
  7.  
  8.  
  9. // [[Rcpp::export]]
  10. arma::uvec intersection(const arma::uvec& x_,const arma::uvec& y_) {
  11.  
  12. std::vector<double> x(x_.begin(),x_.end()),y(y_.begin(),y_.end());
  13.  
  14. std::sort(x.begin(),x.end());std::sort(y.begin(),y.end());
  15.  
  16. std::vector<double> temp;
  17.  
  18. std::set_intersection(x.begin(),x.end(),y.begin(),y.end(),std::back_inserter(temp));
  19.  
  20. arma::uvec out=arma::conv_to<arma::uvec>::from(temp);
  21.  
  22. return out;
  23. }
  24.  
  25. // [[Rcpp::export]]
  26. arma::vec groupSum(const arma::vec& x_,const arma::vec& group_) {
  27.  
  28. uvec index=arma::sort_index(group_);
  29. vec x=x_.elem(index);
  30. vec group=group_.elem(index);
  31. std::vector<double> temp;
  32. int row=0;
  33. temp.push_back(0);
  34. for (int i=0;i<x.size();i++) {
  35. temp[row]+=x[i];
  36. if (group[i]!=group[i+1] && i!=(x.size()-1)) {
  37. row++;
  38. temp.push_back(0);
  39. }
  40. }
  41. return arma::conv_to<arma::vec>::from(temp);
  42. }
  43.  
  44.  
  45. // [[Rcpp::export]]
  46. DataFrame tsTotal1(const NumericVector& value_,const NumericVector& weight_,
  47. const NumericVector& strata_,const NumericVector& cluster_) {
  48. // calculates totals, standard errors for totals, percent, standard error for
  49. // percent for one way contigency table
  50.  
  51. // make arma vecs
  52. vec value(value_), weight(weight_), strata(strata_), cluster(cluster_);
  53. // unique levels and strata
  54. vec levels=unique(value);
  55. vec strats=unique(strata);
  56. // allocate memory for output
  57. vec totals(levels.size());
  58. vec variances(levels.size());
  59. vec percent(levels.size());
  60. vec percentSE(levels.size());
  61. vec outVal(levels.size());
  62. Rcpp::Rcout << "# of Strata: " << strats.size() << std::endl;
  63.  
  64. double total=sum(weight);
  65. Rcpp::Rcout << "Sum of weights: " << total << std::endl;
  66. int totalclusters=0;
  67. // allocate memorey for objects in loop
  68. int clustnum;
  69. double variance,pvar,cn;
  70. vec clustsum, subweights, pclustsum;
  71. // for each level
  72. for (int i=0;i<levels.size();i++) {
  73. // total summed weight
  74. uvec valdelta=arma::find(value==levels[i]);
  75. totals[i]=sum(weight.elem(valdelta));
  76. percent[i]=totals[i]/total;
  77. outVal[i]=levels[i];
  78. // init variance sum
  79. variance=0;pvar=0;
  80. // for each strata
  81. for (int j=0;j<strats.size();j++) {
  82. // clusters nested in strata
  83. uvec stratdelta=arma::find(strata==strats[j]);
  84. vec subclust=unique(cluster.elem(stratdelta)); //shouldn't define vec in loop
  85. clustnum=subclust.size();
  86. if (clustnum==1) {
  87. continue;
  88. }
  89. totalclusters+=clustnum;
  90. cn=clustnum;
  91. clustsum.set_size(clustnum);
  92. pclustsum.set_size(clustnum);
  93. subweights.set_size(stratdelta.size());
  94.  
  95. subweights=weight.elem(stratdelta);
  96. subweights.elem(arma::find(value.elem(stratdelta)!=levels[i])).zeros();
  97. clustsum=groupSum(subweights,cluster.elem(stratdelta));
  98. pclustsum=(clustsum-percent[i]*groupSum(weight.elem(stratdelta),cluster.elem(stratdelta)))/total;
  99. // add variance of strata j
  100. variance+=(cn/(cn-1))*sum(pow(clustsum-mean(clustsum),2));
  101. pvar+=(cn/(cn-1))*sum(pow(pclustsum-mean(pclustsum),2));
  102. }
  103. // take sqrt and append to output vec
  104. variances[i]=sqrt(variance);
  105. percentSE[i]=sqrt(pvar);
  106. }
  107. Rcpp::Rcout << "# of Clusters: " << totalclusters/levels.size() << std::endl;
  108. return DataFrame::create(_["values"]=outVal,
  109. _["totals"]=totals,
  110. _["SE"]=variances,
  111. _["percent"]=percent*100,
  112. _["percentSE"]=percentSE*100);
  113. }
  114.  
  115.  
  116. // [[Rcpp::export]]
  117. DataFrame tsTotal2(const NumericVector& value_,const NumericVector& value2_, const NumericVector& weight_,
  118. const NumericVector& strata_,const NumericVector& cluster_) {
  119. // calculates totals, standard errors for totals, percent, standard error for percent
  120. // column and row percent and their standard errors for a two way contigency table
  121.  
  122. // make arma vecs
  123. vec value(value_), value2(value2_),weight(weight_), strata(strata_), cluster(cluster_);
  124. // unique levels and strata
  125. vec levels=unique(value);
  126. vec levels2=unique(value2);
  127. vec strats=unique(strata);
  128. // allocate memory for output
  129. vec totals(levels.size()*levels2.size());
  130. vec variances(levels.size()*levels2.size());
  131. vec percent(levels.size()*levels2.size());
  132. vec percentSE(levels.size()*levels2.size());
  133. vec rpercent(levels.size()*levels2.size());
  134. vec rpercentSE(levels.size()*levels2.size());
  135. vec cpercent(levels.size()*levels2.size());
  136. vec cpercentSE(levels.size()*levels2.size());
  137. vec outVal(levels.size()*levels2.size());
  138. vec outVal2(levels.size()*levels2.size());
  139. Rcpp::Rcout << "# of Strata: " << strats.size() << std::endl;
  140.  
  141. double total=sum(weight);
  142. Rcpp::Rcout << "Sum of weights: " << total << std::endl;
  143. int totalclusters=0;
  144. // allocate memorey for objects in loop
  145. int clustnum, count;
  146. double variance,pvar,rpvar,cpvar,cn,rtotal,ctotal;
  147. vec clustsum, subweights, pclustsum, rpclustsum,cpclustsum, rdelta, cdelta;
  148. // for each level
  149. count=0;
  150. for (int i=0;i<levels.size();i++) {
  151. for (int k=0;k<levels2.size();k++) {
  152. // total summed weight
  153. uvec valdelta1=arma::find(value==levels[i]);
  154. uvec valdelta2=arma::find(value2==levels2[k]);
  155. totals[count]=sum(weight.elem(intersection(valdelta1,valdelta2)));
  156. percent[count]=totals[count]/total;
  157. rtotal=sum(weight.elem(valdelta1));
  158. ctotal=sum(weight.elem(valdelta2));
  159. rpercent[count]=totals[count]/rtotal;
  160. cpercent[count]=totals[count]/ctotal;
  161. outVal[count]=levels[i];
  162. outVal2[count]=levels2[k];
  163. // init variance sum
  164. variance=0;pvar=0;rpvar=0;cpvar=0;
  165. // for each strata
  166. for (int j=0;j<strats.size();j++) {
  167. // clusters nested in strata
  168. uvec stratdelta=arma::find(strata==strats[j]);
  169. vec subclust=unique(cluster.elem(stratdelta)); //shouldn't define vec in loop
  170. clustnum=subclust.size();
  171. if (clustnum==1) {
  172. continue;
  173. }
  174. totalclusters+=clustnum;
  175. cn=clustnum;
  176. clustsum.set_size(clustnum);
  177. pclustsum.set_size(clustnum);
  178. subweights.set_size(stratdelta.size());
  179. rdelta.set_size(stratdelta.size());
  180. cdelta.set_size(stratdelta.size());
  181.  
  182. subweights=weight.elem(stratdelta);
  183. subweights.elem(arma::find(value.elem(stratdelta)!=levels[i])).zeros();
  184. subweights.elem(arma::find(value2.elem(stratdelta)!=levels2[k])).zeros();
  185. clustsum=groupSum(subweights,cluster.elem(stratdelta));
  186. pclustsum=(clustsum-percent[count]*groupSum(weight.elem(stratdelta),cluster.elem(stratdelta)))/total;
  187. rdelta=weight.elem(stratdelta);
  188. rdelta.elem(arma::find(value.elem(stratdelta)!=levels[i])).zeros();
  189. cdelta=weight.elem(stratdelta);
  190. cdelta.elem(arma::find(value2.elem(stratdelta)!=levels2[k])).zeros();
  191.  
  192. rpclustsum=(clustsum-rpercent[count]*groupSum(rdelta,cluster.elem(stratdelta)))/rtotal;
  193. cpclustsum=(clustsum-cpercent[count]*groupSum(cdelta,cluster.elem(stratdelta)))/ctotal;
  194.  
  195. // add variance of strata j
  196. variance+=(cn/(cn-1))*sum(pow(clustsum-mean(clustsum),2));
  197. pvar+=(cn/(cn-1))*sum(pow(pclustsum-mean(pclustsum),2));
  198. rpvar+=(cn/(cn-1))*sum(pow(rpclustsum-mean(rpclustsum),2));
  199. cpvar+=(cn/(cn-1))*sum(pow(cpclustsum-mean(cpclustsum),2));
  200. }
  201. // take sqrt and append to output vec
  202. variances[count]=sqrt(variance);
  203. percentSE[count]=sqrt(pvar);
  204. rpercentSE[count]=sqrt(rpvar);
  205. cpercentSE[count]=sqrt(cpvar);
  206. count+=1;
  207. }
  208. }
  209. Rcpp::Rcout << "# of Clusters: " << totalclusters/(levels.size()*levels2.size()) << std::endl;
  210. return DataFrame::create(_["values"]=outVal,
  211. _["byValues"]=outVal2,
  212. _["totals"]=totals,
  213. _["SE"]=variances,
  214. _["percent"]=percent*100,
  215. _["percentSE"]=percentSE*100,
  216. _["colPercent"]=cpercent*100,
  217. _["colPercentSE"]=cpercentSE*100,
  218. _["rowPercent"]=rpercent*100,
  219. _["rowPercentSE"]=rpercentSE*100);
  220. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!