2017-03-23 48 views
0

对不起,这是一个混乱,但寻找一些建议。基于向量中的概率应用条件

我试图模拟环境随机性的人口增长。我创建了一个概率向量,我使用的10个补丁大小中的任何一个都会遇到一个突变事件,这会使人口规模减少75%。我想在一个函数中应用它。所以基本上我需要让函数运行并确定下一步的总体大小,但是应用该分块大小的概率来确定在新的总体大小存储在矩阵中之前是否会发生灾难性事件。

所以我基本上要像做了的if/then,但不是定义“如果”参数我想应用存储的概率。我四处搜寻没有太多的运气,但我认为这不可能很难做到。谢谢!

d0 <- c(0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05) # chance that a catastrophic 
# disturbance will reduce population for each patch size, assuming that rates of disturbance 
# are much higher in small patches. 

cat <- .25 # disturbance factor, assume that during catastrophic event 2/3 of animals are 
#removed, I want to multiply this by the population size within each time step at the frequency 
#defined in d0. I could probably make this into a function but I still need to know 
# how to use the probability to decide when to apply it. 

# Ricker model (N_t+1=N_t*exp(r*(1-N_t/K))) 

Ricker = function(nt, r, k0, d0) { #setup the Ricker function 
    nt1 = (nt*exp(r*(1-nt/k0))) # Run Ricker model 
    nt1 = ((nt1*cat)) ### Here I would apply the probability, and when necessary the 
# disturbance factor. I.E. Breeding season happens then there is a very harsh winter 
# and many individuals die. 
    return(nt1) #return the value of (Nt+1) 
} 

for(t in 1:(tf-1)) { #loop through time 
    n[t+1,] = Ricker(n[t,],r,k0) #step through Ricker 
} 

我最终做了类似于@Marius建议的事情,它似乎工作得很好,谢谢大家的意见!

Ricker = function(nt, r, k0, d0) { #setup the Ricker function 
    nt1 = (nt*exp(r*(1-nt/k0))) # Run Ricker model 
for(d in 1:(length(d0))) { # Create loop to test each patch for disturbance probability 
dice_rolls = runif(length(d0)) # Generate random uniforms for each element in d0. 
nt1 = ifelse(dice_rolls < d0, nt1 * cat, nt1) # If the'dice roll' is less than the corresponding element of d0 # the patch experiences the disturbance 
    } 
    return(nt1) #return the value of (Nt+1) 
} 

回答

0

如果我理解你的模型正确,你想:

Ricker = function(nt, r, k0, d0) { #setup the Ricker function 
    nt1 = (nt*exp(r*(1-nt/k0))) # Run Ricker model 
    # Generate random uniforms for each element in d0. If the 
    # 'dice roll' is less than the corresponding element of d0, 
    # the patch experiences the disturbance 
    dice_rolls = runif(length(d0)) 
    nt1 = ifelse(dice_rolls < d0, nt1 * cat, nt1) 
    return(nt1) #return the value of (Nt+1) 
} 

要知道如何工作,看看如何dice_rolls < d0运行一个简单的模拟:你可以看到,它把每个元素d0分别与 产生长期平均水平近似于所需的概率:

d0 <- c(0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05) 

n_catastrophes = numeric(length = length(d0)) 
for (sim_num in 1:1000) { 
    dice_rolls = runif(length(d0)) 
    n_catastrophes = n_catastrophes + (dice_rolls < d0) 
} 
# Number of times each patch had a catastrophe in the simulation 
print(n_catastrophes) 
# Simulated probabilities 
print(n_catastrophes/1000) 
+0

我用你的代码在for循环中,它吐出结果这看起来似乎有道理,但我不完全确定它实际上是通过各自的比例循环。 –

+0

@KatieSmith:尝试在你的'Ricker'函数中添加'print(dice_rolls Marius

0

如果我解释这对,你需要一个代表无灾难事件或灾难性事件的值向量,概率等于d0?

可以品尝矢量c(1.0,0.25)与基于D0加权的概率,即第一个条目是50/50,最后的入口是只有5%灾难性损失...等。

对于每一个循环,stochasticly与画你的猫向量:

cat <- unlist(lapply(d0, function(x) sample(c(1.0, 0.25), size=1, prob=c(1-x, x)))) 

,并坚持在你的循环运行里克模型后