基于 R 中的不同概率模拟 "random walk" 类模型
Simulating a "random walk"-like model based on varying probabilities in R
我对 R 编程还比较陌生。我想模拟一个人在 5x5 网格中的移动,因为网格的环境条件各不相同,而且从一个单元格到另一个单元格的移动是基于环境的他们的近邻的情况。 我想要的这个模拟的最终结果是个体在 x 个时间步之后的位置。
首先,我制作了一个数据框,其中包含网格的 x,y 坐标及其环境条件。然后我根据我的随机环境条件 (v1,v2) 计算了运动阻力及其倒数。
env_cond<-data.frame(x=rep(1:5,5),y=rep(1:5,each=5),v1=rnorm(25),v2=rnorm(25))
env_cond$resistance<- res_surf<- (env_cond [1,3] - env_cond [,3])^2 + (env_cond [1,4]- env_cond [,4])^2
env_cond$inv_res <- 1/env_cond$resistance #where movement is based on inverse resistance
env_cond$cell_num <- 1:25
head (env_cond)
x y v1 v2 resistance inv_res cell_num
1 1 1 1.233266019 0.3554372 0.0000000 Inf 1
2 2 1 0.499331993 0.3780565 0.5391708 1.8546999 2
3 3 1 1.633103368 0.7464020 0.3127234 3.1977142 3
4 4 1 -0.583125893 0.6591043 3.3914933 0.2948554 4
5 5 1 0.929743728 -0.7338991 1.2787793 0.7819958 5
6 1 2 0.009317203 0.2060074 1.5203800 0.6577303 6
>
接下来,我创建了一个邻居矩阵。我假设一个人只能移动到它的 4 个直接邻居,而不能移动到网格上的其他任何地方。它显示了对应于一个单元格的直接 4 个邻居的网格的单元格编号。例如,单元格 1(对应于 x=1,y=1)为北方提供 NA,因为它不能移动到网格 space 上方。
north <- ifelse (env_cond$y==1, NA, env_cond$cell_num-5) #y+1
south <- ifelse (env_cond$y==5, NA, env_cond$cell_num+5) #y-1
west <- ifelse (env_cond$x==1, NA, env_cond$cell_num-1) #x-1
east <- ifelse (env_cond$x==5, NA, env_cond$cell_num+1) #x+1
neighbours <- data.frame(north, south, west, east)
head (neighbours)
north south west east
1 NA 6 NA 2
2 NA 7 1 3
3 NA 8 2 4
4 NA 9 3 5
5 NA 10 4 NA
6 1 11 NA 7
>
我通过首先将邻居的逆电阻值分配给单元格编号来创建概率矩阵。我将 NA 替换为 0 以说明运动的不可能,并将无限大的任意值替换为 10。然后我将这些值转换为概率:
prob_mat <- cbind (env_cond$inv_res [neighbours$north], env_cond$inv_res [neighbours$south],env_cond$inv_res [neighbours$west], env_cond$inv_res [neighbours$east])
colnames(prob_mat) <- c("y+1", "y-1", "x-1", "x+1") #renamed the columns of prob matrix
#changing NA to O
prob_mat[is.na(prob_mat)]<-0
#changing inf to 10
prob_mat [6, 1] <- 10
prob_mat [2, 3] <- 10
prob_mat1 <- matrix (nrow = nrow(prob_mat), ncol=4)
for (i in 1:nrow (prob_mat)) {
prob_mat1 [i,]<- prob_mat[i,]/sum(prob_mat[i,])
head (prob_mat1)
[,1] [,2] [,3] [,4]
[1,] 0.0000000 0.26179048 0.0000000 0.73820952
[2,] 0.0000000 0.01556767 0.7459112 0.23852109
[3,] 0.0000000 0.06208574 0.8092602 0.12865408
[4,] 0.0000000 0.10119069 0.7221972 0.17661214
[5,] 0.0000000 0.39156264 0.6084374 0.00000000
[6,] 0.9246218 0.05608074 0.0000000 0.01929748
此概率矩阵显示每个单元格编号移动到其相邻单元格的概率(不显示那些单独邻居的单元格编号)。从这里开始,我有点卡住了。 我不知道如何从单元格 1 实际模拟一个人的移动(假设每个选择都是独立于前一步做出的,有点像马尔可夫链,其中有不同的移动概率在你当前的步骤)。 我怀疑它与索引有关,但我还没有弄清楚如何管理每个单元格的不同概率。这是我第一次在这里发帖,所以希望这能让 sense/is 可重现。非常感谢任何帮助!
最好的方法可能是编写代码将您的矩阵转换为 25x25 转换矩阵并使用马尔可夫链库,但按原样使用您的设置相当简单:
rand_walk <- function(start,steps){
walk = numeric(steps)
walker = start
for(i in 1:steps){
walk[i] <- walker
walker <- walker + sample(c(-5,5,-1,1),1,prob = prob_mat1[walker,])
}
walk
}
基本思想是,向上或向下移动是对当前单元格数加或减5,向右或向左移动是加或减1,因此从向量c(-5,5,-1,1)
中采样就足够了这 4 个步骤的概率由概率矩阵的相应行给出。
典型输出:
> rand_walk(1,100)
[1] 1 2 1 6 1 2 1 2 1 2 1 2 1 6 1 2 1 2 1 6 1 6
[23] 1 6 1 2 1 2 3 8 9 8 13 12 7 8 7 8 3 8 7 8 7 8
[45] 7 8 7 8 7 8 7 8 7 8 7 12 17 22 21 22 17 12 7 12 7 8
[67] 3 8 13 8 7 12 7 8 9 8 9 8 7 6 7 8 7 2 1 6 1 2
[89] 1 6 1 2 1 2 1 2 1 2 1 6
在这段代码中,我给出了完整的遍历(这对于调试目的很有用)但是你当然可以完全放弃累积矩阵 walk
而只是 return 最后的 walker
.另外,请注意,在此代码中,我记录了 steps
个位置,因此只有 steps - 1
个转换。
我对 R 编程还比较陌生。我想模拟一个人在 5x5 网格中的移动,因为网格的环境条件各不相同,而且从一个单元格到另一个单元格的移动是基于环境的他们的近邻的情况。 我想要的这个模拟的最终结果是个体在 x 个时间步之后的位置。
首先,我制作了一个数据框,其中包含网格的 x,y 坐标及其环境条件。然后我根据我的随机环境条件 (v1,v2) 计算了运动阻力及其倒数。
env_cond<-data.frame(x=rep(1:5,5),y=rep(1:5,each=5),v1=rnorm(25),v2=rnorm(25))
env_cond$resistance<- res_surf<- (env_cond [1,3] - env_cond [,3])^2 + (env_cond [1,4]- env_cond [,4])^2
env_cond$inv_res <- 1/env_cond$resistance #where movement is based on inverse resistance
env_cond$cell_num <- 1:25
head (env_cond)
x y v1 v2 resistance inv_res cell_num
1 1 1 1.233266019 0.3554372 0.0000000 Inf 1
2 2 1 0.499331993 0.3780565 0.5391708 1.8546999 2
3 3 1 1.633103368 0.7464020 0.3127234 3.1977142 3
4 4 1 -0.583125893 0.6591043 3.3914933 0.2948554 4
5 5 1 0.929743728 -0.7338991 1.2787793 0.7819958 5
6 1 2 0.009317203 0.2060074 1.5203800 0.6577303 6
>
接下来,我创建了一个邻居矩阵。我假设一个人只能移动到它的 4 个直接邻居,而不能移动到网格上的其他任何地方。它显示了对应于一个单元格的直接 4 个邻居的网格的单元格编号。例如,单元格 1(对应于 x=1,y=1)为北方提供 NA,因为它不能移动到网格 space 上方。
north <- ifelse (env_cond$y==1, NA, env_cond$cell_num-5) #y+1
south <- ifelse (env_cond$y==5, NA, env_cond$cell_num+5) #y-1
west <- ifelse (env_cond$x==1, NA, env_cond$cell_num-1) #x-1
east <- ifelse (env_cond$x==5, NA, env_cond$cell_num+1) #x+1
neighbours <- data.frame(north, south, west, east)
head (neighbours)
north south west east
1 NA 6 NA 2
2 NA 7 1 3
3 NA 8 2 4
4 NA 9 3 5
5 NA 10 4 NA
6 1 11 NA 7
>
我通过首先将邻居的逆电阻值分配给单元格编号来创建概率矩阵。我将 NA 替换为 0 以说明运动的不可能,并将无限大的任意值替换为 10。然后我将这些值转换为概率:
prob_mat <- cbind (env_cond$inv_res [neighbours$north], env_cond$inv_res [neighbours$south],env_cond$inv_res [neighbours$west], env_cond$inv_res [neighbours$east])
colnames(prob_mat) <- c("y+1", "y-1", "x-1", "x+1") #renamed the columns of prob matrix
#changing NA to O
prob_mat[is.na(prob_mat)]<-0
#changing inf to 10
prob_mat [6, 1] <- 10
prob_mat [2, 3] <- 10
prob_mat1 <- matrix (nrow = nrow(prob_mat), ncol=4)
for (i in 1:nrow (prob_mat)) {
prob_mat1 [i,]<- prob_mat[i,]/sum(prob_mat[i,])
head (prob_mat1)
[,1] [,2] [,3] [,4]
[1,] 0.0000000 0.26179048 0.0000000 0.73820952
[2,] 0.0000000 0.01556767 0.7459112 0.23852109
[3,] 0.0000000 0.06208574 0.8092602 0.12865408
[4,] 0.0000000 0.10119069 0.7221972 0.17661214
[5,] 0.0000000 0.39156264 0.6084374 0.00000000
[6,] 0.9246218 0.05608074 0.0000000 0.01929748
此概率矩阵显示每个单元格编号移动到其相邻单元格的概率(不显示那些单独邻居的单元格编号)。从这里开始,我有点卡住了。 我不知道如何从单元格 1 实际模拟一个人的移动(假设每个选择都是独立于前一步做出的,有点像马尔可夫链,其中有不同的移动概率在你当前的步骤)。 我怀疑它与索引有关,但我还没有弄清楚如何管理每个单元格的不同概率。这是我第一次在这里发帖,所以希望这能让 sense/is 可重现。非常感谢任何帮助!
最好的方法可能是编写代码将您的矩阵转换为 25x25 转换矩阵并使用马尔可夫链库,但按原样使用您的设置相当简单:
rand_walk <- function(start,steps){
walk = numeric(steps)
walker = start
for(i in 1:steps){
walk[i] <- walker
walker <- walker + sample(c(-5,5,-1,1),1,prob = prob_mat1[walker,])
}
walk
}
基本思想是,向上或向下移动是对当前单元格数加或减5,向右或向左移动是加或减1,因此从向量c(-5,5,-1,1)
中采样就足够了这 4 个步骤的概率由概率矩阵的相应行给出。
典型输出:
> rand_walk(1,100)
[1] 1 2 1 6 1 2 1 2 1 2 1 2 1 6 1 2 1 2 1 6 1 6
[23] 1 6 1 2 1 2 3 8 9 8 13 12 7 8 7 8 3 8 7 8 7 8
[45] 7 8 7 8 7 8 7 8 7 8 7 12 17 22 21 22 17 12 7 12 7 8
[67] 3 8 13 8 7 12 7 8 9 8 9 8 7 6 7 8 7 2 1 6 1 2
[89] 1 6 1 2 1 2 1 2 1 2 1 6
在这段代码中,我给出了完整的遍历(这对于调试目的很有用)但是你当然可以完全放弃累积矩阵 walk
而只是 return 最后的 walker
.另外,请注意,在此代码中,我记录了 steps
个位置,因此只有 steps - 1
个转换。