小编Inf*_*tor的帖子

我开始非常喜欢小提琴情节,因为当你有趣的发行时,它们给我一个更好的感觉.我喜欢自动化很多东西,因此遇到了一个问题:当一个变量的方差为0时,boxplot只会给你一条线.然而,Geom_violin以错误终止.我喜欢什么样的行为？好吧,要么排成一行,要么没有,但请给我其他变量的分布.

好的,快速的例子:

dff=data.frame(x=factor(rep(1:2,each=100)),y=c(rnorm(100),rep(0,100)))
ggplot(dff,aes(x=x,y=y)) + geom_violin()

产量

Error in `$<-.data.frame`(`*tmp*`, "n", value = 100L) : 
  replacement has 1 row, data has 0

但是,有效的是:

ggplot(dff,aes(x=x,y=y)) + geom_boxplot()

更新:

该问题从昨天开始解决:https://github.com/hadley/ggplot2/issues/972

更新2 :(来自问题作者)哇,哈德利自己回应了!geom_violin现在表现与geom_densityR基本一致density.

但是,我不认为这种行为是最优的.

(1)'零'问题

只需使用我的原始示例运行它:

dff=data.frame(x=factor(rep(1:2, each=100)), y=c(rnorm(100), rep(0,100)))
ggplot(dff,aes(x=x,y=y)) + geom_violin(trim=FALSE)

产生这个:

右边的情节是否是"全零"的适当表示？我不这么认为.最好是修剪产生一条线以显示数据没有变化.解决方法解决方案:添加+geom_boxplot()

(2)我可能真的想要TRIM=TRUE.

例:

dff=data.frame(x=factor(rep(1:2, each=100)), y=c(rgamma(100,1,1), rep(0,100)  ))
ggplot(dff,aes(x=x,y=y)) + geom_violin(trim=FALSE)

现在我有非零数据,标准内核密度估计不能正确处理.随着trim=T我可以很快看到数据是严格积极的.

我并不认为当前的行为是"错误的",因为它与其他功能一致.但是,geom_violin可以在不同的上下文中使用,用于探索具有异构数据类型的不同data.frames(例如,正面+倾斜或不正面).

更新:

我试图实施德克的建议.评论？我现在正忙于JSM,但我想在为画廊编织Rmd之前得到一些反馈.我从犰狳换回正常的Rcpp,因为它没有增加任何价值.带有R ::的标量版本非常好.如果将mean/sd作为标量输入,而不是作为所需输出长度的向量,我应该在参数n中输入绘制数量.

有许多MCMC应用程序需要从截断的Normal分布中绘制样本.我建立在TN的现有实现上并添加了并行计算.

问题:

1. 有没有人看到进一步的速度提升？在基准测试的最后一种情况下,rtruncnorm有时会更快.Rcpp实现总是比现有的包更快,但是可以进一步改进吗？
2. 我在一个我无法分享的复杂模型中运行它,我的R会话崩溃了.但是,我不能系统地重现它,所以它可能是代码的另一部分.如果有人在TN工作,请测试并告诉我.更新:我没有更新代码的问题,但请告诉我.

我如何把事情放在一起:据我所知,最快的实现不在CRAN上,但源代码可以下载OSU stat.在我的基准测试中,msm和trunco​​rm中的竞争实现较慢.诀窍是有效地调整提案分布,其中指数很好地适用于截断的Normal的尾部.所以我拿了Chris的代码,"Rcpp'ed"它并添加了一些openMP香料.动态调度在这里是最佳的,因为取样可以根据边界花费更多或更少的时间.我发现一件令人讨厌的事情:当我想使用双打时,许多统计分布基于NumericVector类型.我只是编写了我的方式.

继承人Rcpp代码:

#include <Rcpp.h>
#include <omp.h>


// norm_rs(a, b)
// generates a sample from a N(0,1) RV restricted to be in the interval
// (a,b) via rejection sampling.
// ======================================================================

// [[Rcpp::export]]

double norm_rs(double a, double b)
{
   double  x;
   x = Rf_rnorm(0.0, 1.0);
   while( (x < a) || (x > b) ) x = norm_rand();
   return x;
}

// half_norm_rs(a, b)
// generates …

我有很多代码使用base :: order()命令,我真的懒得在rcpp中编写代码.由于Rcpp只支持排序,但不支持订单,我花了2分钟创建这个功能:

// [[Rcpp::export]]
Rcpp::NumericVector order_cpp(Rcpp::NumericVector invec){
  int leng = invec.size();
  NumericVector y = clone(invec);
  for(int i=0; i<leng; ++i){
    y[sum(invec<invec[i])] = i+1;
  }
  return(y);
}

它有点工作.如果向量包含唯一数字,我得到与order()相同的结果.如果它们不是唯一的,结果是不同的,但没有错(实际上没有唯一的解决方案).

使用它:

c=sample(1:1000,500)
all.equal(order(c),order_cpp(c))
microbenchmark(order(c),order_cpp(c))

Unit: microseconds
         expr      min       lq   median       uq      max neval
     order(c)   33.507   36.223   38.035   41.356   78.785   100
 order_cpp(c) 2372.889 2427.071 2466.312 2501.932 2746.586   100

哎哟! 我需要一个有效的算法.好的,所以我挖出了一个bubbleort实现并对其进行了调整:

 // [[Rcpp::export]]
Rcpp::NumericVector bubble_order_cpp2(Rcpp::NumericVector vec){                                  
       double tmp = 0;
       int n = vec.size();
              Rcpp::NumericVector outvec = clone(vec);
       for (int i …

我在png图像文件上设置了一个带有透明叠加散点图的漂亮绘图.我希望我的绘图窗口和我的pdf输出与我的png-962x745完全相同.

但是,即使在关闭轴,注释和帧之后,R仍然会在图像周围留下边框.

这可以用一个简单的例子来显示:该图显示了两个点,它们应位于图的最外端.但它们不是:

plot(rbind(c(1,745),c(962,1)),bty ="n",axes=F,frame.plot=F, xaxt='n', ann=FALSE, yaxt='n', asp=745/962)

和PDF设备一起:

pdf(width=10.02,height=7.76)
par(mar=rep(0, 4),mai=rep(0, 4), xpd = NA) 
plot(rbind(c(1,745),c(962,1)),bty ="n",axes=F,frame.plot=F, xaxt='n', ann=FALSE, yaxt='n', asp=745/962)
dev.off()

该解决方案现已在Rcpp Gallery中联机

我从RcppArmadillo的mvtnorm包中重新实现了dmvnorm.我不知何故喜欢犰狳,但我想它也适用于普通的Rcpp.来自dmvnorm的方法基于马哈拉诺比斯距离,所以我有一个函数,然后是多元正态密度函数.

让我告诉你我的代码:

#include <RcppArmadillo.h>
#include <Rcpp.h>

// [[Rcpp::depends("RcppArmadillo")]]

// [[Rcpp::export]]
arma::vec mahalanobis_arma( arma::mat x ,  arma::mat mu, arma::mat sigma ){

  int n = x.n_rows;
  arma::vec md(n);
    for (int i=0; i<n; i++){
        arma::mat x_i = x.row(i) - mu;
        arma::mat Y = arma::solve( sigma, arma::trans(x_i) );
        md(i) = arma::as_scalar(x_i * Y);
    }
    return md;

    }



// [[Rcpp::export]]
arma::vec dmvnorm ( arma::mat x,  arma::mat mean,  arma::mat sigma, bool log){ 

arma::vec distval = mahalanobis_arma(x,  mean, sigma);

    double logdet …

在Rcpp中以数字方式计算以下内容的最佳方法是什么？

exp(-1500)/(exp(-1500)+exp(-1501))

在许多情况下,计算可能需要多精度(对于exp),但最终结果可以舍入到通常的double.

通过quadmath？通过提升？

如果你留在R(在Rcpp之外),那里有非常舒适的包装工作:

library(Rmpfr)

a = mpfr(-1500,100)
b = mpfr(-1501,100)

exp(a)/(exp(a)+exp(b))

但是如何使用rcpp访问？