Início / coding / Perguntas / 79133654
Accepted
farrow90
Asked: 2024-10-28 21:42:21 +0800 CST

Gerando todas as sequências e probabilidades para um processo aleatório em R

  • 772

Tenho uma Cadeia de Markov em 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))

Parece algo assim:

          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

Para um número fixo de voltas, quero descobrir todas as sequências de estados possíveis que podem ocorrer e suas probabilidades correspondentes.

Tentei fazer isso manualmente usando loops para enumerar todas essas sequências:

# 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)
}

Então tentei chamar a função:

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

Existe algo que eu possa fazer para garantir que esse código seja executado mais rápido para um grande número de voltas

PS: Usei esse código para verificar se todas as probabilidades somam 1 em cada turno:

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 respostas

  1. Best Answer

    Acho que deve funcionar bem para o seu propósito (dado tm <- transition_matrixo nome da variável mais curto para "matriz de transição")

    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
    )
})

    e você pode ver

    > 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

    Outra opção ( pode ser mais eficiente, pois as linhas que resultam em 0probabilidades são filtradas antecipadamente ) é usar mergedentro de um repeatloop, para atualizar a tabela, iterativamente

    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
}

    tal que

    > 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

relate perguntas