【问题标题】:Calculating stats for random subsample using R loop使用 R 循环计算随机子样本的统计数据
【发布时间】:2015-08-31 08:23:34
【问题描述】:

我正在尝试在 R 中找到一种方法来随机子集一些数据(生态研究区域中合适栖息地的比例),计算值 > 0 的样本的平均值和比例,然后将这些值保存或附加到一个数据框。然后我想重复多次(例如 1000 次)。标准引导或重采样包不起作用,因为我需要计算出现频率以及子样本的平均值。我知道“应用”函数,但是这些函数会在整个数据帧上循环,而我试图在重复的子样本上执行它。我知道我需要一些代码来获取循环中的计算值以保存和输出但有问题。 “habprop”是数据框(“数据”)中的一列,我想计算并保存正值的平均值和比例。

for(i in 1000 {  
randsample=data[sample(1:nrow(data), 50, replace=FALSE),]
m=mean(randsample$habprop)
randsamplepos=subset(randsample, habprop > 0)
habfreq=(nrow(randsamplepos)/nrow(randsample))
})

【问题讨论】:

  • 为什么标准引导程序不起作用?您可以将任意函数传递给boot 并返回值列表
  • 在循环外创建一个空列表lst <- list(),然后在底部lst[i] <- habfreq。如果habfreq 是每次迭代一个值,您可以简化为v1 <- c() 和v1[i]
  • 另外,如果你不打算使用m,为什么要在for循环中计算1000次呢?
  • 最好创建一个reproducible example,其中包含样本输入和所需的输出。描述一下您所说的“有问题”是什么意思——究竟什么是行不通的?

标签: r for-loop subsampling


【解决方案1】:

使用boot 这应该是可能的

dat <- data.frame(habprop=rnorm(100))

## Function to return statistics from subsamples
stat <- function(dat, inds)
    with(dat, c(mu=mean(habprop[inds]), freq=sum(habprop[inds] > 0)/length(inds)))

library(boot)
boot(data=dat, statistic=stat, R=1000)

# Bootstrap Statistics :
#        original      bias    std. error
# t1* -0.06154533 -0.00324393  0.08377116
# t2*  0.52000000 -0.00073000  0.04853991

【讨论】:

  • 在这种情况下,我理解引导问题的方式是您从分布中提取并基于该分布进行计算(在本例中为正常分布。) - 设置第一个统计信息因此是负数,所有值都是正数。
【解决方案2】:

replicate 函数怎么样? This post 看起来很相似。

生成一些数据来处理

data &lt;- data.frame(x1=rpois(5000, 5), x2=runif(5000), x3=rnorm(5000))

定义一个函数来采样和获取均值和计数

sample_stats <- function(df, n=100){
  df <- df[sample(1:nrow(df), n, replace=F),]
  mx1 <- mean(df$x1[df$x1>0])
  x1pos <- sum(df$x1>0)
  return(c(mx1, x1pos))
}

运行一次只是为了查看输出

sample_stats(data)

运行 1000 次

results &lt;- replicate(1000, sample_stats(data, n=100))

【讨论】:

  • 感谢@ajb 这段代码似乎有效——你知道在函数中设置 df 时是否“sample”每次都会抽取一个随机样本吗?
  • @kevinT 是的sample 每次都会抽取一个随机样本。如果你想让结果可重复(确保你可以生成相同的随机样本),你可以在运行sample之前直接使用set.seed
猜你喜欢
  • 2018-10-26
  • 1970-01-01
  • 1970-01-01
  • 2021-12-19
  • 1970-01-01
  • 2019-01-24
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多