66766.py

import sys

from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
from scipy.stats import *

p17 = .05
p18 = .03

# 17, 18, any other state
pval = array([p17, p18, 1. - p17 - p18])
vals = array([17, 18, 0])

def main(argv):
	lengths = range(5, 251, 10)
	probs = []

	for T in lengths:
		counter = 0
		repetitions = 2000

		for _ in range(repetitions):
			past = 0

			for t in range(T):
				# draw current state
				current = vals[multinomial(1, pval) > .5]

				if current == 18:
					if past == 17:
						# 18, 17, 18 detected
						counter += 1
						break
					else:
						# 18 detected
						past = 18

				elif current == 17 and past == 18:
					# 18, 17 detected
					past = 17

		probs.append(counter / float(repetitions))

	ion()
	figure(figsize=(6, 6))
	plot(lengths, probs, linewidth=2)
	xlabel('Length of the sequence (T)')
	ylabel('Probability of subsequence')

	raw_input()

	return 0



if __name__ == '__main__':
	sys.exit(main(sys.argv))