r - 用 R 在对数尺度上转换可变密度

Question

我想绘制范围如下的变量的密度：

 Min.   :-1214813.0  
 1st Qu.:       1.0  
 Median :      40.0  
 Mean   :     303.2  
 3rd Qu.:     166.0  
 Max.   : 1623990.0

密度的线性图导致 [0,1000] 范围内的高柱，有两个非常长的尾向正无穷大和负无穷大。因此，我想将变量转换为对数刻度，以便我可以看到平均值附近发生了什么。例如，我正在考虑类似的事情：

log_values = c( -log10(-values[values<0]), log10(values[values>0]))

这导致：

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-6.085   0.699   1.708   1.286   2.272   6.211 

这样做的主要问题是它不包含这些0值。当然，我可以将所有值从0with移开values[values>=0]+1，但这会在数据中引入一些失真。

将这个变量转换为对数刻度的可接受且科学可靠的方法是什么？

Answer 1

你所拥有的基本上就是@James 所建议的。这对于 (-1,1) 中的值是有问题的，尤其是那些接近原点的值：

x <- seq(-2, 2, by=.01)
plot(x, sign(x)*log10(abs(x)), pch='.')

这样的事情可能会有所帮助：

y <- c(-log10(-x[x<(-1)])-1, x[x >= -1 & x <= 1], log10(x[x>1])+1)

plot(x, y, pch='.')

这是连续的。可以通过使用区间 (-1/log(10), 1/log(10)) 强制 C^1，该区间通过求解 d/dx log10(x) = 1 得到：

z <- c( -log10(-x[x<(-1/log(10))]) - 1/log(10)+log10(1/log(10)),
         x[x >= -1/log(10) & x <= 1/log(10)],
         log10(x[x>1/log(10)]) + 1/log(10)-log10(1/log(10))
       )
plot(x, z, pch='.')

Answer 2

除了转换之外，您还可以操纵直方图本身来了解您的数据。这为您提供了绘图本身保持可读性的优势，并且您可以立即了解中心的分布。假设我们模拟以下数据：

Data <- c(rnorm(1000,5,10),sample(-10000:10000,10))
> summary(Data)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-9669.000    -2.119     5.332    85.430    12.460  9870.000 

然后你有几种不同的方法。最容易看到数据中心正在发生的事情，就是绘制数据中心。在这种情况下，假设我对第一个和第三个四分位数之间发生的事情感兴趣，我可以绘制：

hist(Data,
     xlim=c(-30,30),
     breaks=c(min(Data),seq(-30,30,by=5),max(Data))
     main="Center of Data"
     )

如果您还想计算尾部，您可以转换数据以折叠尾部并更改轴以反映这一点，如下所示：

1. 您将感兴趣范围之外的所有值分配给该范围之外的值
2. 您绘制直方图，将所有极值合并到一个 bin 中
3. 您使用正确的标签构建 X 轴
4. 您axis.break()从包plotrix中使用在 X 轴上添加一些中断，表示不连续轴

为此，您可以使用以下代码：

 require(plotrix)
 # rearrange data
 plotdata <- Data
 id <- plotdata < -30 | plotdata > 30
 plotdata[id] <- sign(plotdata[id])*35
 # plot histogram
 hist(plotdata,
      xlim=c(-40,40),
      breaks=c(-40,seq(-30,30,by=5),40),
      main="Untailed Data",
      xaxt='n'   # leave the X axis away
      )
 # Construct the X axis
 axis(1,
      at=c(-40,seq(-30,30,by=10),40),
      labels=c(min(Data),seq(-30,30,by=10),max(Data))
 )
 # add axis breaks
 axis.break(axis=1,breakpos=-35)
 axis.break(axis=1,breakpos=35)

这给了你：

请注意，您可以通过添加freq=TRUE到hist()函数来获得原始频率。

Answer 3

我将此添加为另一个答案，因为尽管想法相似，但映射从根本上不同。

当对数标度图中包含小值 (<1) 时，它是典型图，log(1 + .)而不是log(.).

反思原点，我们得到一些有用的东西：

x <- seq(-2, 2, by=.01)   
w <- c( -log10(1-x[x<0]), x[x==0], log10(1+x[x>0]))

plot(x, w, pch='.')

应该清楚该函数是平滑的，因为在 0 附近的方向导数也会被反映。

在 x 中有更大的值：

x <- seq(-10000, 10000, by=.01)
w <- c( -log10(1-x[x<0]), x[x==0], log10(1+x[x>0]))
plot(x, w, pch='.')