CompBio_2_Hw3/main.py at master · Umbriona/CompBio_2_Hw3 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
import numpy as np
import matplotlib.pyplot as plt
import Exercises as ex
import scipy.special
import math

def main():
    populationSize = 2*187
    ex.Question1(populationSize)

    N = 1000
    r = 0.01
    R =1
    listTimeA = []
    listTimeB = []
    numberOfRepeats = 100000
    for i in range(numberOfRepeats):
        timeA,timeB = ex.Question2(2,R)
        listTimeA.append(timeA)
        listTimeB.append(timeB)
    listTimeA = np.asarray(listTimeA)
    listTimeB = np.asarray(listTimeB)
    meanTimeA = np.mean(listTimeA)
    meanTimeB = np.mean(listTimeB)
    mean2TimeA = np.mean(listTimeA**2)
    mean2TimeB = np.mean(listTimeB**2)

    fig, axes = plt.subplots(nrows=2, ncols=1)
    ax0, ax1 = axes.flatten()
    #n_binsA = len(set(listTimeA))
    #n_binsB = len(set(listTimeB))
    def analytT(x):
        recombinationRate = 0
        coalescentRate = scipy.special.comb(2, 2)
        Lambda = recombinationRate + coalescentRate
        return Lambda * np.exp(-Lambda*x)
    x = np.linspace(0,17,100)
    meanTime = np.mean(analytT(x))
    mean2Time = np.mean(analytT(x) ** 2)
    print('<Ta>', meanTimeA, '\n<Tb>', meanTimeB, '\n<Ta^2', mean2TimeA, '\n<Tb^2', mean2TimeB,'\n<T>',meanTime,'\n<T^2>',mean2Time)
    ax0.hist(listTimeA, bins=100,density = True, histtype = 'bar', facecolor='red',label='Simulation')
    ax0.plot(x,analytT(x),label='Analytical',color='blue')
    ax0.set_title('T_a R = 1')
    ax0.legend(loc='upper center', shadow=True, fontsize='x-large')

    ax1.hist(listTimeB, bins=100, density = True , facecolor='red',label='Simulation')
    ax1.plot(x, analytT(x),label='Analytical',color='blue')
    ax1.legend(loc='upper center', shadow=True, fontsize='x-large')

    ax1.set_title('T_b R = 1')
    plt.show()

    R = np.array([0,10**-5,10**-4,10**-3,10**-2,0.1,1,10])
    n = np.array([2])
    R_a = np.linspace(0,10,100)
    analyt = lambda x: (x+18)/(x**2+13*x+18)
    numberOfRepeats = 100000
    fig, axes = plt.subplots(nrows=2, ncols=1)
    ax = axes.flatten()
    for k in range(len(n)):
        Pab = []
        Pab_c = []
        for j in range(len(R)):
            listTimeA = []
            listTimeB = []
            for i in range(numberOfRepeats):
                timeA, timeB = ex.Question2(n[k], R[j])
                listTimeA.append(timeA)
                listTimeB.append(timeB)
            listTimeA = np.asarray(listTimeA)
            listTimeB = np.asarray(listTimeB)
            meanTimeA = np.mean(listTimeA)
            meanTimeB = np.mean(listTimeB)
            meanTimeAB = np.mean(listTimeA*listTimeB)
            varTimeA = np.var(listTimeA)
            varTimeB = np.var(listTimeB)
            Pab.append((meanTimeAB-meanTimeA*meanTimeB)/(varTimeA*varTimeB))
            Pab_c.append(np.cov([listTimeA,listTimeB]))
        norm = np.amax(Pab)
        Pab = Pab/norm
        ax[k].set_xscale('log')
        ax[k].plot(R,Pab,'*',label='sim')
        ax[k].plot(R, analyt(R),label='analytic')
        ax[k].set_ylabel('Pab')
        ax[k].set_xlabel('R')
        ax[k].legend(loc='upper right', shadow=True)
        ax[k].set_title('N = '+str(n[k]))
    plt.show()


if __name__ == '__main__':
    main()