主页 / coding / 问题 / 79133654
Accepted
farrow90
Asked: 2024-10-28 21:42:21 +0800 CST

在 R 中生成随机过程的所有序列和概率

  • 772

我在 R 中有一个马尔可夫链：

set.seed(123)

n_states <- 5

matrix <- matrix(runif(n_states^2), nrow=n_states)

# set some transitions to 0 
matrix[1, 4:5] <- 0 
matrix[5, 1:3] <- 0  
matrix[2, 5] <- 0    

transition_matrix <- t(apply(matrix, 1, function(x) x/sum(x)))

rownames(transition_matrix) <- paste0("S", 1:n_states)
colnames(transition_matrix) <- paste0("S", 1:n_states)

print(round(transition_matrix, 3))

它看起来像这样：

          S1         S2        S3         S4        S5
S1 0.2229340 0.03531601 0.7417500 0.00000000 0.0000000
S2 0.3910569 0.26197885 0.2248868 0.12207747 0.0000000
S3 0.1536622 0.33530265 0.2545791 0.01580275 0.2406533
S4 0.2652280 0.16563210 0.1719994 0.09849610 0.2986444
S5 0.0000000 0.00000000 0.0000000 0.59278229 0.4072177

对于固定的回合数，我想找出所有可能出现的状态序列及其对应的概率。

我尝试使用循环手动枚举所有这样的序列：

# Function to generate sequences for multiple turn lengths
find_sequences_all_turns <- function(transition_matrix, start_state = 1, max_turns = 5) {
    n_states <- nrow(transition_matrix)
    
   
    all_sequences <- list()
    all_probabilities <- numeric()
    all_turns <- numeric()
    seq_counter <- 1
    
    generate_sequence <- function(current_seq, current_prob, steps_left, total_steps) {
        if(length(current_seq) > 1) {
            all_sequences[[seq_counter]] <<- current_seq
            all_probabilities[seq_counter] <<- current_prob
            all_turns[seq_counter] <<- total_steps - steps_left
            seq_counter <<- seq_counter + 1
        }
        
        if(steps_left == 0) {
            return()
        }
        
        current_state <- current_seq[length(current_seq)]
        possible_next_states <- which(transition_matrix[current_state,] > 0)
        
        for(next_state in possible_next_states) {
            prob <- transition_matrix[current_state, next_state]
            generate_sequence(
                c(current_seq, next_state),
                current_prob * prob,
                steps_left - 1,
                total_steps
            )
        }
    }
    
    generate_sequence(c(start_state), 1, max_turns - 1, max_turns)
    
    result_df <- data.frame(
        turn = all_turns,
        sequence_no = 1:length(all_sequences),
        sequence = sapply(all_sequences, paste, collapse=""),
        probability = all_probabilities
    )
    
    result_df <- result_df[order(result_df$turn, -result_df$probability),]
    rownames(result_df) <- NULL
    
    return(result_df)
}

然后我尝试调用该函数：

sequences_df <- find_sequences_all_turns(transition_matrix)

> sequences_df
    turn sequence_no sequence  probability
1      2         154       13 7.417500e-01
2      2           1       11 2.229340e-01
3      2          65       12 3.531601e-02
4      3         171      132 2.487108e-01
5      3         193      133 1.888341e-01
6      3         243      135 1.785046e-01
7      3          40      113 1.653613e-01
8      3         155      131 1.139789e-01
9      3           2      111 4.969955e-02
10     3          66      121 1.381057e-02
11     3         218      134 1.172169e-02
12     3          82      122 9.252048e-03
13     3         104      123 7.942105e-03
14     3          18      112 7.873137e-03
15     3         129      124 4.311289e-03

我可以做些什么来确保此代码在大量回合中运行得更快

PS：我使用此代码来验证每次所有概率的总和是否为 1：

library(dplyr)

probability_sums <- sequences_df %>%
    group_by(turn) %>%
    summarise(
        total_probability = sum(probability),
        num_sequences = n(),
        check_sum_to_one = abs(total_probability - 1) < 1e-10
    )

print(probability_sums)

1 个回答

  1. Best Answer

    我认为它应该能很好地满足你的目的（给定tm <- transition_matrix“转换矩阵”的较短变量名）

    nms <- row.names(tm)
N <- 5
res <-lapply(2:N, function(k) {
    u <- expand.grid(rep(list(nms), k))
    subset(
        transform(
            u,
            Freq = Reduce(`*`, Map(\(x, y) tm[cbind(x, y)], u[-length(u)], u[-1]))
        ),
        Freq > 0
    )
})

    你可以看到

    > lapply(res, head, n = 10)
[[1]]
   Var1 Var2       Freq
1    S1   S1 0.22293395
2    S2   S1 0.39105686
3    S3   S1 0.15366217
4    S4   S1 0.26522803
6    S1   S2 0.03531601
7    S2   S2 0.26197885
8    S3   S2 0.33530265
9    S4   S2 0.16563210
11   S1   S3 0.74175004
12   S2   S3 0.22488682

[[2]]
   Var1 Var2 Var3       Freq
1    S1   S1   S1 0.04969955
2    S2   S1   S1 0.08717985
3    S3   S1   S1 0.03425651
4    S4   S1   S1 0.05912833
6    S1   S2   S1 0.01381057
7    S2   S2   S1 0.10244863
8    S3   S2   S1 0.13112240
9    S4   S2   S1 0.06477157
11   S1   S3   S1 0.11397892
12   S2   S3   S1 0.03455660

[[3]]
   Var1 Var2 Var3 Var4        Freq
1    S1   S1   S1   S1 0.011079716
2    S2   S1   S1   S1 0.019435348
3    S3   S1   S1   S1 0.007636940
4    S4   S1   S1   S1 0.013181713
6    S1   S2   S1   S1 0.003078844
7    S2   S2   S1   S1 0.022839277
8    S3   S2   S1   S1 0.029231635
9    S4   S2   S1   S1 0.014439781
11   S1   S3   S1   S1 0.025409771
12   S2   S3   S1   S1 0.007703838

[[4]]
   Var1 Var2 Var3 Var4 Var5         Freq
1    S1   S1   S1   S1   S1 0.0024700449
2    S2   S1   S1   S1   S1 0.0043327990
3    S3   S1   S1   S1   S1 0.0017025332
4    S4   S1   S1   S1   S1 0.0029386513
6    S1   S2   S1   S1   S1 0.0006863789
7    S2   S2   S1   S1   S1 0.0050916503
8    S3   S2   S1   S1   S1 0.0065167238
9    S4   S2   S1   S1   S1 0.0032191175
11   S1   S3   S1   S1   S1 0.0056647005
12   S2   S3   S1   S1   S1 0.0017174471

    另一个选项（可能更有效，因为导致0概率的行被提前过滤掉）是merge在repeat循环中使用，迭代地更新表

    N <- 5
dp <- d <- subset(as.data.frame.table(tm), Freq > 0)
res <- list(dp)
repeat {
    if (length(res) == N - 1) break
    p <- names(dp[-length(dp)])
    q <- paste0("Var", as.integer(sub("\\D+", "", p)) + 1)
    dnew <- setNames(dp, c(q, "yFreq"))
    dp <- subset(
        transform(
            merge(d, dnew)[unique(c(p, names(dnew), "Freq"))],
            Freq = Freq * yFreq
        ),
        select = -yFreq
    )
    res[[length(res) + 1]] <- dp
}

    使得

    > lapply(res, head, n = 10)
[[1]]
   Var1 Var2       Freq
1    S1   S1 0.22293395
2    S2   S1 0.39105686
3    S3   S1 0.15366217
4    S4   S1 0.26522803
6    S1   S2 0.03531601
7    S2   S2 0.26197885
8    S3   S2 0.33530265
9    S4   S2 0.16563210
11   S1   S3 0.74175004
12   S2   S3 0.22488682

[[2]]
   Var1 Var2 Var3        Freq
1    S1   S1   S1 0.049699546
2    S1   S1   S2 0.007873137
3    S1   S1   S3 0.165361267
4    S2   S1   S1 0.087179851
5    S2   S1   S2 0.013810568
6    S2   S1   S3 0.290066443
7    S3   S1   S1 0.034256514
8    S3   S1   S2 0.005426735
9    S3   S1   S3 0.113978919
10   S4   S1   S1 0.059128333

[[3]]
   Var1 Var2 Var3 Var4         Freq
1    S1   S1   S1   S1 0.0110797161
2    S1   S1   S1   S2 0.0017551896
3    S1   S1   S1   S3 0.0368646403
4    S1   S1   S3   S2 0.0554460705
5    S1   S1   S3   S3 0.0420975207
6    S1   S1   S3   S4 0.0026131624
7    S1   S1   S3   S5 0.0397947423
8    S1   S1   S2   S4 0.0009611327
9    S1   S1   S2   S1 0.0030788444
10   S1   S1   S2   S2 0.0020625955

[[4]]
   Var1 Var2 Var3 Var4 Var5         Freq
1    S1   S1   S1   S1   S1 0.0024700449
2    S1   S1   S1   S1   S2 0.0003912914
3    S1   S1   S1   S1   S3 0.0082183799
4    S1   S1   S1   S3   S2 0.0123608115
5    S1   S1   S1   S3   S3 0.0093849666
6    S1   S1   S1   S3   S4 0.0005825626
7    S1   S1   S1   S3   S5 0.0088715991
8    S1   S1   S1   S2   S4 0.0002142691
9    S1   S1   S1   S2   S1 0.0006863789
10   S1   S1   S1   S2   S2 0.0004598226
    • 2

相关问题