How to randomly change the transition matrix using by Markov chain?

Question

I'm using a 1st step Transition matrix to generate the DNA sequences. Now I need to give a probability to the transition matrix to change every 1000 steps. Let's say, every 1000 steps, there is 40% probability the transition matrix will change. Every row should add to 1 after the change. Now I don't know how to access the value in the nested dictionary data in python, and how to implement the 40% probability changing. I attached my code here, any suggestion is appreciated ~

#!/usr/bin/env python

import sys, random


length = 10000

tran_matrix = {'a': {'a':0.495,'c':0.113,'g':0.129,'t':0.263},
               'c': {'a':0.129,'c':0.063,'g':0.413,'t':0.395},
               't': {'a':0.213,'c':0.495,'g':0.263,'t':0.029},
               'g': {'a':0.263,'c':0.129,'g':0.295,'t':0.313}}

initial_p = {'a':0.25,'c':0.25,'t':0.25,'g':0.25}             

def choose(dist):
    r = random.random()
    sum = 0.0
    keys = dist.keys()
    for k in keys:
        sum += dist[k]
        if sum > r:
        return k
    return keys[-1]
c = choose(initial_p)
for i in range(length):
    sys.stdout.write(c) 
    c = choose(tran_matrix[c])

Marius · Answer 1 · 2012-12-10T04:39:53.060

EDIT: Added a quick implementation of some code that generates new transition frequencies. You may have to play around to find which random number generator works best for your case, and to see if you can get more sensible generation using some constraints on the random probabilities.

import sys, random


LENGTH = 10000
CHANGE_EVERY = 1000
CHANGE_PROB = 0.4

tran_matrix = {'a': {'a':0.495,'c':0.113,'g':0.129,'t':0.263},
               'c': {'a':0.129,'c':0.063,'g':0.413,'t':0.395},
               't': {'a':0.213,'c':0.495,'g':0.263,'t':0.029},
               'g': {'a':0.263,'c':0.129,'g':0.295,'t':0.313}}

initial_p = {'a':0.25,'c':0.25,'t':0.25,'g':0.25}             


def choose(dist):
    r = random.random()
    sum = 0.0
    keys = dist.keys()
    for k in keys:
        sum += dist[k]
        if sum > r:
            return k
    return keys[-1]


def new_probs(precision=2):
    """
    Generate a dictionary of random transition frequencies, of the form
    {'a':0.495,'c':0.113,'g':0.129,'t':0.263}
    """
    probs = []
    total_prob = 0
    # Choose a random probability p1 from a uniform distribution in
    # the range (0, 1), then choose p2 in the range (0, 1 - p1), etc.
    for i in range(3):
        up_to = 1 - total_prob
        p = round(random.uniform(0, up_to), precision)
        probs.append(p)
        total_prob += p
    # Final probability is 1 - (sum of first 3 probabilities)
    probs.append(1 - total_prob)
    # Assign randomly to bases
    # If you don't shuffle the order of the bases each time, 't'
    # would end up with consistently lower probabilities
    bases = ['a', 'c', 'g', 't']
    random.shuffle(bases)
    new_prob_dict = {}
    for base, prob in zip(bases, probs):
        new_prob_dict[base] = prob
    return new_prob_dict

c = choose(initial_p)
for i in range(LENGTH):
    if i % CHANGE_EVERY == 0:
        dice_roll = random.random()
        if dice_roll < CHANGE_PROB:
            for base in tran_matrix:
                # Generate a new probability dictionary for each
                # base in the transition matrix
                tran_matrix[base] = new_probs()
    sys.stdout.write(c) 
    c = choose(tran_matrix[c])

Hi Marius, thanks for the help, yes this should works for me, the change_matrix() function should randomly generate 4 float numbers sum up to 1. For instance 'a': {'a':0.495,'c':0.113,'g':0.129,'t':0.263}, the four numbers add to 1, I just need 4 new random number to replace the former 4 numbers. — Frank, Dec 10 '12 at 04:10
@Frank: Have a look at the `new_probs()` function I've added- you might have to try different ways to generate the random numbers, as my implementation can give long strings of the same base. — Marius, Dec 10 '12 at 04:35
That works perfectly, thank you so much Marius, you are fantastic, I will click the useful button when I got enough reputation. Have a nice day!!!! — Frank, Dec 10 '12 at 04:47

How to randomly change the transition matrix using by Markov chain?

1 Answers1