r - 如何在满足特定均匀度的情况下生成一组随机生成的丰度？

Question

我想生成 100 个随机数(表示物种丰度)，条件是这些数字等于给定的物种均匀度。

例如，我想创建一个 100x1 矩阵 (M)，其中 -sum(p_i*log(p_i))/log(100) = 0.78

其中，p_i = M[i,1]/sum(M[,1]) 并且 0.78 是所需的物种均匀度。

感谢任何帮助!

Answer 1

我做了一些事情。我不知道它是否适合你。这是一个很小的蒙特卡洛。
但首先，我必须从香农指数开始。你的方程式似乎有错误。我用了这个等式

Jet 是我的H功能。

H = function(x) -sum(x/sum(x)/log(x/sum(x)))

我承认，我继续创建了一个略微延伸的特殊示例函数。

ssample = function(n, n10){
  if(sum(n10)>1) stop("sum n10>1")
  if(sum(n10<1)) n10=c(n10, 1-sum(n10))
  x = as.numeric()
  for(i in 1:length(n10)) {
    x = c(x, sample(1:(10^i), n*n10[i], replace = T))}
  x
}

如果我们想找到H大约等于 0.78，此函数必须接收一个特殊的权重参数 n10 .

n=100
n10 = c(0.5, 0.25, 0.1, 0.08, 0.05, 0.02)

就在前面。我们运行蒙特卡洛 10,000 次迭代。

library(tidyverse)
set.seed(1111)
df = tibble(n1 = 1:10000) %>%
  mutate(x = map(n1, ~ssample(n, n10))) %>%
  mutate(H = map(x, ~H(.x))) %>%
  unnest(H)

df %>% filter(H>0.7799 & H<0.7801) %>%
  arrange(H)

输出

# A tibble: 1 x 3
     n1 x               H
  <int> <list>      <dbl>
1  7550 <dbl [100]> 0.780

如您所见，一个向量带有 H在 0.7799 之间和 0.7801被画出来了。
让我们在情节上看到它。

fgat = function(x) tibble(n = x, gat = 1:length(x))
df %>% filter(H>0.7799 & H<0.7801) %>% 
  mutate(x = map(x, ~fgat(.x))) %>%
  unnest(x) %>%
  ggplot(aes(gat, n))+
  geom_point()+
  ggtitle(paste("Random species abundance with Shannon Index = 0.780"))+
  scale_y_continuous(trans='log10')

现在让我们检查一下我们绘制的向量在 H> 0.779 & H <0.781 范围内的样子.

df %>% filter(H>0.779 & H<0.781) %>% arrange(H) %>% 
  mutate(H = H %>% round(4) %>%  paste() %>% fct_inorder()) %>%
  mutate(x = map(x, ~fgat(.x))) %>% 
  unnest(x) %>% 
  ggplot(aes(gat, n, color=H))+
  geom_point()+
  facet_grid(vars(H))+
  scale_y_continuous(trans='log10')

它对你来说够好吗？

@Skaqqs 的更新

library(microbenchmark)

f = function(n, nmc, n10){
  df = tibble(n1 = 1:nmc) %>%
    mutate(x = map(n1, ~ssample(n, n10))) %>%
    mutate(H = map(x, ~H(.x))) %>%
    unnest(H)
}

ggplot2::autoplot(microbenchmark(f(100, 10000, n10), times=100))

@Dr_brachiopod 的更新

我不得不承认，寻找解决方案让我有些疲惫。然而，最终结果非常好。要回答你的问题，什么是 n10 .它是一个“权重”向量，它决定了多少百分比的抽签将在 1:10 范围内，有多少 1:100，有多少 1:1000 等。请参见下表。这应该清除一切。

该图显示了随机选择的分布密度 H预期值为 H = 0.78 的值标记。如您所见，并非 n10 的所有组合适合绘制此特定值。

当涉及到用某个 H 取样时值，使用简单的 filter ，例如像这样df %>% filter(H>0.779 and H<0.781) .