Coin toss and Absorbing Markov Chain

http://www.mitbbs.com/mitbbs_bbsdoc_div_article.php?board=Quant&gid=31212871


A fair coin is flipped until the first time one of the following two
patterns appeared: TTH or HTH
ask: which one you should choose for a better chance to win


+ TT = .5 TTH + .5 TT
| HT = .5 HTH + .5 TT
| T = .5 TT + .5 H
| H = .5 HT + .5 H
+ S = .5 T + .5 H

If we let TTH = 1 and HTH = 0, we finally get S = 5/8 (the probability of reaching TTH first)
or otherwise
we let HTH = 1 and TTH = 0, we finally get S = 3/8 (the probability of reaching HTH first)



We start from the definition of Absorbing Markov Chain.
REF: Chap 11, Markov Chains, from www.dartmouth.edu/~chance/teaching


A state s_i of a Markov chain is called absorbing if it is impossible to leave it (i.e., p_ii = 1). A Markov chain is absorbing if it has at lease one absorbing state, and if from every state it is possible to go to an absorbing state (not necessarily in one step). The state that is not absorbing is called transient.

The interesting questions about an absorbing Markov chain are:

- What is the probability that the process will eventually reach a given absorbing state?
- On the average, how long will it take for the process to be absorbed?
- On the average, how many times will the process be in each transient state?

In general, the answers to all questions depends on
- the state from which the process starts
- the transition probabilities.
The answers to the above questions are the theorems listed below.

We first write the transition matrix in canonical Form by explicitly write Q and R.

      TT           T         S         H         HT
Q =

TT    0.5000         0         0         0         0
 T    0.5000         0         0    0.5000         0
 S         0    0.5000         0    0.5000         0
 H         0         0         0    0.5000    0.5000
HT    0.5000         0         0         0         0


       TTH       HTH
R =   

TT    0.5000         0
 T         0         0
 S         0         0
 H         0         0
HT         0    0.5000

Theorem: For an absorbing Markov chain the matrix I-Q has an inverse N = I + Q + Q^2 + .... The ij-entry of the matrix N is the expected number of times the chain is in state s_j, given that it starts in state s_i.

>> N = inv(eye(5)-Q)

N =

    2.0000         0         0         0         0
    1.5000    1.0000         0    1.0000    0.5000
    1.2500    0.5000    1.0000    1.5000    0.7500
    1.0000         0         0    2.0000    1.0000
    1.0000         0         0         0    1.0000

By reading the first line of N, we know that starting from TT, the expected number of times we still in TT is 2.

Theorem: The sum of row i of N is the expected number of steps before the chain is absorbed given that the chain starts in state s_i.

>> t = N*ones(5,1)

t =

     2
     4
     5
     4
     2

So starting from S, the expected number of steps before absorbing in either TTH or HTH is 5.

Theorem: Let b_ij be the probability that an absorbing chain will be absorbed in the absorbing state s_j if it starts in the transient state s_i. Let B be the matrix with entries b_ij. Then B is an txr matrix, and

B = NR,

where N is the fundamental matrix and R is as in the canonical form.

>> B = N*R

B =

    1.0000         0
    0.7500    0.2500
    0.6250    0.3750
    0.5000    0.5000
    0.5000    0.5000

Thus starting from S, we have probability of 0.625 = 5/8 to be absorbed into state TTH and probability of 0.375 = 3/8 to be absorbed into state HTH.



Now let us consider adding another absorbing state THT. The transition graph is as below:

Q =      S         T         H        TH        HT        TT

         0    0.5000    0.5000         0         0         0     S
         0         0         0    0.5000         0    0.5000     T
         0         0    0.5000         0    0.5000         0     H
         0         0    0.5000         0         0         0     TH
         0         0         0         0         0    0.5000     HT
         0         0         0         0         0    0.5000     TT

R =    THT       HTH       TTH

         0         0         0  S
         0         0         0  T
         0         0         0  H
    0.5000         0         0  TH
         0    0.5000         0  HT
         0         0    0.5000  TT

>> N = inv(eye(6)-Q)

N =

    1.0000    0.5000    1.2500    0.2500    0.6250    1.1250
         0    1.0000    0.5000    0.5000    0.2500    1.2500
         0         0    2.0000         0    1.0000    1.0000
         0         0    1.0000    1.0000    0.5000    0.5000
         0         0         0         0    1.0000    1.0000
         0         0         0         0         0    2.0000

>> N*ones(6,1)

ans =

    4.7500 S
    3.5000 T
    4.0000 H
    3.0000 TH
    2.0000 HT
    2.0000 TT

>> B = N*R

B =    THT       HTH       TTH

    0.1250    0.3125    0.5625  S
    0.2500    0.1250    0.6250  T
         0    0.5000    0.5000  H
    0.5000    0.2500    0.2500  TH
         0    0.5000    0.5000  HT
         0         0    1.0000  TT
 

gcd, binomial number

Problem: UVa 530 Binomial Showdown

这道题我用了cache + online computation, 一小部分数据先算好，不在表里的现算，在表内的直接查表

typedef unsigned long long uLL;
uLL gcd(uLL a, uLL b) {
    if(a % b == 0) 
        return b;
    else
        return gcd(b, a%b);
}

void divbygcd(uLL &a, uLL &b) {
    uLL g = gcd(a,b);
    a /= g;
    b /= g;
}

uLL combination(int n, int k) {
    uLL num = 1, den = 1, tomul, todiv, i;
    if(k > n/2) k = n-k; //use smaller k
    for(i = k; i; --i) {
        tomul = n-k+i;
        todiv = i;
        divbygcd(tomul, todiv);
        divbygcd(num, todiv);
        divbygcd(tomul, den);
        num *= tomul;
        den *= todiv;
    }
    return num / den;
}

Prime Numbers

Sieve of Eratosthenes: 

#include <cmath>
const int MAX = 1000000;
const int SQRT_MAX = (int)sqrt(MAX);

short primes[MAX];
void gen_primes() 
{ 
  for(int i=0;i<MAX;i++) 
    primes[i] = 1; 
  for(int i=2; i<= SQRT_MAX;i++) 
    if (primes[i]) 
      for(int j=i;j*i<MAX;j++) 
        primes[i*j] = 0; 
} 

or use our bitvector to store primes[]:

Problem: UVa 543: Goldbach's Conjecture

#include <iostream>
const int MAX = 1000001;
const int SQRT_MAX = 1000;

const int SHIFT = 5;
const int MASK = (1<<SHIFT)-1;

inline void SET(unsigned int a[], int i) {
	a[i>>SHIFT] |= (1 << (i & MASK));
}
inline void CLEAR(unsigned int a[], int i) {
	a[i>>SHIFT] &= ~(1 << (i & MASK));
}
inline bool TEST(unsigned int a[], int i) {
	return a[i>> SHIFT] & (1 << (i & MASK));
}

const unsigned int N = MAX/(sizeof(int)<<3);
unsigned int primes[N+1];

void gen_primes() 
{ 
    memset(primes, 0xAAAAAAAA, sizeof(primes));
    SET(primes, 2);
    for(int i = 3; i <= SQRT_MAX; ++ i) {
        if (TEST(primes, i))  {
            for(int j = i; j * i < MAX; ++ j)  {
                CLEAR(primes, i * j);
            }
        }
    }
} 

int main() {
    gen_primes();
    int n;
    while(scanf("%d ", &n) && n) {
        unsigned int i = 2, pos = 0;
        while(true) {
            unsigned int p = n-i;
            if(TEST(primes,i) && TEST(primes, p)) {             
                printf("%d = %d + %d\n", n, i, p);
                break;
            }
            while(!TEST(primes,++i)) {};
        }
    }
}

O(logn) Fibonacci using Matrix property

long long divide_conquer_fibonacci(long n) {
 long long i, h, j, k ,t;
 i = h = 1;
 j = k = 0;
 while( n > 0) {
     if(n % 2 == 1) {
         t = j * h;
         j = i * h + j * k + t;
         i = i * k + t;
     }
     t = h * h;
     h = 2 * k * h + t;
     k = k * k + t;
     n = (long) n / 2;
 }
 return j;
}

long long data type (-2^63-1 to 2^63-1] can hold until Fib(92).

1, 1
2, 1
3, 2
4, 3
5, 5
6, 8
7, 13
8, 21
9, 34
10, 55
11, 89
12, 144
13, 233
14, 377
15, 610
16, 987
17, 1597
18, 2584
19, 4181
20, 6765
.....
86, 420196140727489673
87, 679891637638612258
88, 1100087778366101931
89, 1779979416004714189
90, 2880067194370816120
91, 4660046610375530309
92, 7540113804746346429   ----- (19 digits)
93, -6246583658587674878  ----- (overflow)

Using the Fast exponentiation to solve linear recurrences, some similar problems follow here.

timeWatch:

#include "time.h"
unsigned int start = clock();
// .... your code
std::cout << "Time taken in millisecs: " << clock()-start;