<jupyter_start><jupyter_text>Problem 1Given: $N \sim Poisson(\lambda)$ and $X_1, \dots, X_n \sim \vec{\pi}$$X_k(t)$ is continous time MC with $X_k(0) = X_k$$N_t(a) = $\{k:X_k(t) = a\}$i.e. $N_t$ is the number of visits to state $a$ in time $t$.$\sum_a\pi(a)Q_{ab}=0$ for each $b$ with the constraint $\sum_a\pi(a)=1$$\sum_a\pi(a)Q_{ab}=0$ $\implies$ $\vec{\pi}^TQ=0$ $\implies$ $$\begin{align*}\vec{\pi}^TQ&=0\\\Longleftrightarrow \vec{\pi}^TQ^n&=0\ \  \forall n \geq 1\\ \Longleftrightarrow \sum_{n\geq 1}\vec{\pi}\frac{t^n}{n!}Q^n &=0 \ \  \forall t \geq 0\\\Longleftrightarrow \vec{\pi}\sum_{n\geq 0}\frac{t^n}{n!}Q^n &=\vec{\pi}\\\Longleftrightarrow \vec{\pi}P &=\vec{\pi}\\\Longleftrightarrow \vec{\pi}\  \text{is a stationary distribution}\end{align*}$$Now, $P(X_k(t)=a)=\pi(a)$ and $N_t(a) = \{k:X_k(t) = a\}$ $\implies$ $N_t(a)|N \sim Binom(N, \pi(a))$ and $N \sim Poisson(\lambda)$ then $\boxed{N_t \sim Poisson(\lambda \pi)}$  Problem 2<jupyter_code>%matplotlib inline
from __future__ import division
import pandas as pd
import matplotlib
import itertools
matplotlib.rcParams['figure.figsize'] = (16,12)
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(1)

def propose(S):
    r = np.random.choice(len(S), 2)
    rs = np.sort(r)
    j,k=rs[0],rs[1]
    y=np.copy(S)
    y[j:k+1] = y[j:k+1][::-1]
    return y

def count_cycles(S):
    sample_length = len(S)
    n_cycles = 0
    index = 0
    length_travelled = 0
    visited = []
    while length_travelled < sample_length:
        if S[index] == index and index < sample_length :
            index+=1
            n_cycles+=1
            length_travelled+=1
        else:
            visited.append(index)
            index = S[index]
            length_travelled+=1
            if index not in visited:
                n_cycles+=1
    return n_cycles
N = [2,3,4, 100]
alpha = 3
assert count_cycles([0,1]) == 2
assert count_cycles([0,2,1]) == 2
assert count_cycles([1,0]) == 1
N_iterations = 10000

def theoretical(S, alpha, denom):
    n_cycles = count_cycles(S)
    return n_cycles**alpha/denom


def run(n, show=True):
    oldS = np.arange(n)
    old_n_cycles = count_cycles(oldS)
    count_dict = {}
    if show:
        denom = sum([count_cycles(x)**alpha  for x in itertools.permutations(range(n))])
    for i in range(N_iterations):
        proposedS = propose(oldS)
        new_n_cycles = count_cycles(proposedS)
        pi_ab = new_n_cycles**alpha/(old_n_cycles**alpha)
        q = min(1,pi_ab)
        if q>= np.random.uniform():
            oldS = proposedS
            old_n_cycles = new_n_cycles
        tkey = ','.join([str(x+1) for x in oldS.tolist()])
        key="["+tkey+"]"
        if key not in count_dict:
            if show:
                count_dict[key] = [0,0,0]
                count_dict[key][1] = theoretical(oldS,alpha,denom)
                count_dict[key][2] = old_n_cycles
            else:
                count_dict[key] = [0,0]
                count_dict[key][1] = old_n_cycles


        count_dict[key][0]+=1
    df = pd.DataFrame(count_dict)
    df=df.transpose()
    if show:
        df.columns=[r'Simulated $\pi(s)$', 'Theoretical', 'c(s)']
        df[r'Simulated $\pi(s)$'] = df[r'Simulated $\pi(s)$']/N_iterations
        df['Percentage Error'] = 100*(df[r'Simulated $\pi(s)$']/df['Theoretical']-1)
    else:
        df.columns=[r'Simulated $\pi(s)$', 'c(s)']
        df[r'Simulated $\pi(s)$'] = df[r'Simulated $\pi(s)$']/N_iterations
        
    df.index.name='State'
    return df<jupyter_output><empty_output><jupyter_text>n=2<jupyter_code>df  = run(N[0])
df<jupyter_output><empty_output><jupyter_text>n=3<jupyter_code>df = run(N[1])
df<jupyter_output><empty_output><jupyter_text>n=4<jupyter_code>count_dict = run(N[2])
count_dict<jupyter_output><empty_output><jupyter_text>N=100<jupyter_code>df  = run(N[3], show=False)<jupyter_output><empty_output><jupyter_text>$$\sum_{s \in S_a}\pi(s)c(s)=E[c(s)]$$and similarly, $$\sum_{s \in S_a}\pi(s)c^2(s)=E[c^2(s)]=Var(c(s))+E^2[c(s)]$$<jupyter_code>expectation = sum(df[r'Simulated $\pi(s)$']*df['c(s)'])
expectation2 = sum(df[r'Simulated $\pi(s)$']*df['c(s)']*df['c(s)'])
t_expectation = np.mean(df['c(s)'])
t_expectation2 = np.var(df['c(s)'])+np.mean(df['c(s)'])**2
print 'Simulated E[c(s)] = {}\t\t Theoretical(M0M) E[c(s)] = {}'.format(expectation, t_expectation)
print 'Simulated E[c^2(s)] = {}\t\t Theoretical(M0M) E[c^2(s)] = {}'.format(expectation2, t_expectation2)<jupyter_output>Simulated E[c(s)] = 79.2362		 Theoretical(M0M) E[c(s)] = 74.9541157482
Simulated E[c^2(s)] = 6536.9272		 Theoretical(M0M) E[c^2(s)] = 5929.13534701<jupyter_text>I use method of moments to calculate the Theoretical values. They seem to be in sync with the simulated values. The estimates seem to be in sync even though MOM is just a first approximation because the sample size is large enough to capture the dynamics of the population distribution.<jupyter_code>cycles = df['c(s)']
plt.hist(cycles, normed=True)<jupyter_output><empty_output><jupyter_text>Problem 3  Part (A)<jupyter_code>def run():
    N = 1000
    N_iterations = 200
    chrom_length = 3*(10**9)
    transposon_length = 3*1000
    mu = 0.05
    t_positions = []

    n_initial = np.random.random_integers(N-1)
    x_initial = np.random.random_integers(chrom_length-1)

    offspring_positions = []
    all_positions = [[] for t in range(N)]
    all_positions[n_initial].append(x_initial)
    all_t_count =[]

    for nn in range(N_iterations):
        for i in range(N):
            indicator = np.random.binomial(1,mu,len(all_positions[i]))
            temp_indices = []
            for ind, ind_value in enumerate(indicator):
                if ind_value == 1:
                    temp_indices.append(ind)

            for j in temp_indices:
                x_temp = np.random.random_integers(chrom_length-1)
                all_positions[i][j] = x_temp
                all_positions[i].append(np.random.random_integers(chrom_length-1))
        offspring_positions = [[] for t in range(N)]
        for j in range(N):
            y,z = np.random.random_integers(0,N-1,2)
            y_parent = np.random.binomial(1,0.5,len(all_positions[y]))
            z_parent = np.random.binomial(1,0.5,len(all_positions[z]))
            temp_y = []
            temp_z = []
            for index,value in enumerate(y_parent):
                if value>=1:
                      temp_y.append(all_positions[y][index])
            for index,value in enumerate(z_parent):
                if value>=1:
                      temp_z.append(all_positions[z][index])
            for t_y in temp_y:
                offspring_positions[j].append(t_y)
            for t_z in temp_z:
                offspring_positions[j].append(t_z)
        all_positions = offspring_positions
        count_t = 0
        count_x = []
        for p in range(N):
            count_t += len(all_positions[p])
            count_x.append(all_positions[p])
        survived_t = np.unique(count_x, return_counts=True)[1]
        all_t_count.append((count_t, len(survived_t[survived_t>=N*mu])))
    return all_t_count
all_t_count = run()
die_out_transposons = all_t_count
fig, axs = plt.subplots(2,2)
axs[0][0].plot([x[0] for x in die_out_transposons])
axs[0][0].set_title('No. of Transpososn v/s Generations')
axs[0][0].set_xlabel('Generations')
axs[0][0].set_ylabel('No. of Transpososn')

axs[0][1].plot([x[1] for x in die_out_transposons])
axs[0][1].set_title('No. of Common Transpososn v/s Generations')
axs[0][1].set_xlabel('Generations')
axs[0][1].set_ylabel('No. of Common Transposons')


increasing_rate = []
for i in range(1,len(die_out_transposons)):
    increasing_rate.append(die_out_transposons[i][0]/(die_out_transposons[i-1][0]+0.000001))
axs[1][0].plot(increasing_rate)
axs[1][0].set_title('Increasing rate v/s Generations')
axs[1][0].set_xlabel('Generations')
axs[1][0].set_ylabel('Increasing Rate')<jupyter_output><empty_output><jupyter_text>The above example shows one case when the "the transposon does not spread"<jupyter_code>all_t_count = run()
nondie_out_transposons = all_t_count
fig, axs = plt.subplots(2,2)
axs[0][0].plot([x[0] for x in nondie_out_transposons])
axs[0][0].set_title('No. of Transpososn v/s Generations')
axs[0][0].set_xlabel('Generations')
axs[0][0].set_ylabel('No. of Transpososn')

axs[0][1].plot([x[1] for x in die_out_transposons])
axs[0][1].set_title('No. of Surviving Transpososn v/s Generations')
axs[0][1].set_xlabel('Generations')
axs[0][1].set_ylabel('No. of Common Transposons')


increasing_rate = []
for i in range(1,len(die_out_transposons)):
    increasing_rate.append(die_out_transposons[i][0]/(die_out_transposons[i-1][0]+0.000001))
axs[1][0].plot(increasing_rate)
axs[1][0].set_title('Increasing rate v/s Generations')
axs[1][0].set_xlabel('Generations')
axs[1][0].set_ylabel('Increasing Rate')<jupyter_output><empty_output><jupyter_text>The above example shows one case when the "the transposon does spread", with rate being exponential and the common transposons still being limited  Part (B)Treating the total number of transposons $N(t)$ at any time $t$ to be a branching process, then $N(t+1) = \sum_{i=1}^{N(t)} W_{t,i}$ Where $W_{t,i}$ is the number of locations of the $i^{th}$ transposon in the offspring.Now consider $E[N_t]$Claim: $E[N_t] = (1+\mu)^t$Proof: With probability of $\mu$ the transposon undergoes becomes 2 from 1. and hence $W_{t,i}$, the number  of locations of the $i^{th}$ transposon in the offspring is a poisson random varaible with mean $1+\mu$$W_k \sim Poisson(1+\mu)$$$\begin{align*}E[N_t] = E[\sum_{i=1}^{N(t-1)} W_{t,i}] &= E[E[\sum_{i=1}^{n} W_{t,i}|N(t-1)=n]]\\ &= E[N(t-1)] \times (1+\mu)\\&= N(1+\mu)^t\end{align*}$$Thus, the expected number of total transposons is an exponential.  Part (C)$P(X>0) \leq EX$Consider $X_t$ as the total number of trasposon copies at location $x$ at generation $t$ For each new generation, the new arrival at $x$ is a poisson process with mean = $2 \times \mu \times N(t) \times \frac{1}{L}$. Let $R(t)$ represent the new arrivals at x$N(t)$ represents the number of transposon copies of the transposon suriving for $t$ generations!Then, $E[N_1]=1+\mu$ By indution, $E[N_t] = (1+\mu)^t$Thus, $R(t) \sim \text{Poisson}(\frac{2\mu N(t)}{L})$Now Using a branching process model for number of transposon copies located at location $x$, the offspring mean number of offspring transposons is = $(1-\mu)*1 + \mu*2 = 1+\mu$Let $Z_{t,k}(u)$= Number of offspring copies of $k^{th}$ transposon at $x$ that occured at time $t$ inserted at time u ($u \leq t)Using branching process property, $E[Z_{t,k}(t+u)] = (1+\mu)^u$Then $$\begin{align*}EX_t &= \sum_{u \leq 0} \sum_{k=1}^{R(u)} Z_{u,k}(0)\\ &= \sum_{u \leq 0}E[R(u)]E[Z_{u,1}(0)]\\ &= \sum_{u \leq 0} (1+\mu)^t \frac{2 N \mu }{L} \times (1+\mu)^u \\&=  \sum_{u\geq 0}(1+\mu)^t \frac{2 N \mu }{L(1+\mu)^u}\\&\approx \frac{2 \mu }{L}\times(1+\frac{1}{\mu})\\&= \frac{2}{L}(1+\mu)^{t+1}\end{align*}$$Thus, $P(X>0) \leq \frac{2}{L}(1+\mu)^{t+1}$  Part (D) $\mu = 10^{-2}$$N = 10^7$For an individual $ EX = \frac{2}{L}(1+\mu)^{t+1} \times \frac{1}{N} = \frac{2}{NL}(1+\mu)^{t+1} $Now, $\frac{2}{NL}(1+\mu)^{t+1}=0.1L$ $\implies$ $(1+\mu)^{t+1}=0.1NL^2/2$<jupyter_code>from math import log
N=10**7
mu=0.01
L=3*(10**9)
t = log(0.1*N*L*L/2)/log(1+mu)
print(t)<jupyter_output>5704.94795262
{'fraction_non_alphanumeric': 0.18136086841408525, 'fraction_numerical': 0.03115347277380637, 'mean_word_length': 3.7673538073201516, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '500000', 'license': 'MIT', 'n_tokens_mistral': 4688, 'n_tokens_neox': 4423, 'n_words': 896, 'path': '2015_Fall/MATH-578B/Homework5/Homework5.ipynb', 'repo_name': 'NeveIsa/hatex', 'size': '11331'}