<jupyter_start><jupyter_text>![image.png](attachment:image.png)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.integrate import quad
import seaborn as sns
from numba import jit
import time
plt.style.use('ggplot')
@jit
def Fokker_Planck_Equation(x, N, u, v, s):
    """ 
    returns a function for stationary distribution of Moran process.
    """
    return np.power(x, N*v - 1)*np.power((1 - x), (N*u - 1))*np.exp(N*s*x)

@jit
def Fokker_Planck_Mean(x, N, u, v, s):
    """
    returns the mean of Moran process at stationary using Fokker-Planck.
    """
    return x*np.power(x, N*v - 1)*np.power((1 - x), (N*u - 1))*np.exp(N*s*x)

@jit
def Fixation_Prob_Moran(N, w_a, w_A):
    """
    returns an approximation for fixation probability for Moran process.
    """
    return ((w_A/w_a) - 1)/(np.power(w_A/w_a, N) - 1)

@jit
def Fixation_Prob_Fisher(N, s):
    """
    returns an approximation for fixation probability for Wright-Fisher process.
    """
    return (1 - np.exp(-s))/(1 - np.exp(-N*s))

@jit
def Moran_Process(N, w_a, w_A, u, v,
                  num_a, num_A, max_steps,
                  dt, num_sim):
    """
    N -- population size,
    w_a, w_A -- weights,
    u, v -- mutation rates,
    dt -- time step.
    """
    step = 0
    a_counts = []
    while step < max_steps:
        prob_birth = ((w_a * num_a)*(1 - u) / (w_a * num_a + w_A * num_A) 
                      + (num_A * v) / (w_a * num_a + w_A * num_A)) * num_A / N
        prob_death = ((num_A)*(1 - v) / (w_a * num_a + w_A * num_A) 
                      + (w_a * num_a)*u / (w_a * num_a + w_A * num_A)) * num_a / N
        r = np.random.rand(1)
        if r <= prob_birth*dt:
            num_a += 1
            num_A -= 1
        elif r <= (prob_birth + prob_death)*dt:
            num_a -= 1
            num_A += 1
        step += 1    
        a_counts.append(num_a)
    return a_counts[-1]
start = time.time()

population_sizes = [10, 100, 1000]
mutation_rates = [0.2, 0.02, 0.002]
s_list = [0.2, 0.02, 0.002]
w_A = 1
max_steps = 500000
dt = 1
num_sims = 1000
frequency = []
Fokker_Planck = []

for i in range(len(population_sizes)):
    x = np.linspace(0, 1, population_sizes[i])
    P = Fokker_Planck_Equation(x, population_sizes[i], mutation_rates[i],
                               mutation_rates[i], s_list[i])
    print("Population Size: ", population_sizes[i])
    simulations = [Moran_Process(population_sizes[i], 1+s_list[i], w_A, 
                                 mutation_rates[i], mutation_rates[i], 
                                 population_sizes[i]//2, 
                                 population_sizes[i]-(population_sizes[i]//2),
                                 max_steps, dt, j) for j in range(num_sims)]
    Fokker_Planck.append(P)
    a_counts = np.asarray(simulations)
    frequency.append(a_counts)
    
end = time.time()
print('total running time = %1.2f' %(end-start), 'sec') 
fig, ax = plt.subplots(1,len(population_sizes), figsize=(20,5))
fig.suptitle('Stationary Distribution For Different Population Sizes',
             fontsize=18, y=1.12)
fig.text(0.5, -0.02, 'Frequency of a', ha='center', fontsize=16)
fig.text(-0.02, 0.5, 'Density', va='center', rotation='vertical', fontsize=16)
for i in range(len(population_sizes)):
    x = np.linspace(0, 1, population_sizes[i])
    data = frequency[i]/population_sizes[i]
    FP = Fokker_Planck[i]
    ax[i].plot(x, 2*FP, label='1', color='r')
    sns.kdeplot(data, ax=ax[i], label='2', color='b')
    ax[i].axvline(np.mean(data), color='teal', linestyle='dashed',
                  linewidth=2, alpha=0.75)
    ax[i].legend(['Fokker Planck solution', 'Simulations'], fontsize=14)
    ax[i].set_title('N = {0}'.format(population_sizes[i]), fontsize=14)
    ax[i].set_xlim([-0.01,1.01])
    plt.tight_layout()
start = time.time()

population_sizes = [10, 100, 1000, 10000]
mutation_rates = [0.01, 0.001, 0.0001, 0.00001]
s_list = [0.2, 0.02, 0.002, 0.0002]
w_A = 1
max_steps = 500000
dt = 1
num_sims = 1000
frequency = []
Fokker_Planck = []

for i in range(len(population_sizes)):
    x = np.linspace(0, 1, population_sizes[i])
    P = Fokker_Planck_Equation(x, population_sizes[i], mutation_rates[i],
                               mutation_rates[i], s_list[i])
    print("Population Size: ", population_sizes[i])
    simulations = [Moran_Process(population_sizes[i], 1+s_list[i], w_A,
                                 mutation_rates[i], mutation_rates[i], 
                                 population_sizes[i]//2,
                                 population_sizes[i]-(population_sizes[i]//2),
                                 max_steps, dt, j) for j in range(num_sims)]
    Fokker_Planck.append(P)
    a_counts = np.asarray(simulations)
    frequency.append(a_counts)
    
end = time.time()
print('total running time = %1.2f' %(end-start), 'sec') 
fig, ax = plt.subplots(1, 3, figsize=(20,5))
fig.suptitle('Stationary Distribution For Different Population Sizes',
             fontsize=18, y=1.12)
fig.text(0.5, -0.02, 'Frequency of a', ha='center', fontsize=16)
fig.text(-0.02, 0.5, 'Density', va='center', rotation='vertical', fontsize=16)

x = np.linspace(0, 1, population_sizes[-1])
integral = quad(Fokker_Planck_Equation, 0, 1, args=(1, 0.1, 0.1, 0.2))
first_moment = quad(Fokker_Planck_Mean, 0, 1, args=(1, 0.1, 0.1, 0.2))
analytic_mean = first_moment[0]/integral[0]
FP_density = Fokker_Planck[-1]/integral[0]

for i in range(len(population_sizes)-1):
    data = frequency[i]/population_sizes[i]
    print('Difference in mean =%1.2f' %(np.mean(data)-analytic_mean))
    ax[i].plot(x, FP_density, label='1', color='r')
    ax[i].axvline(x=analytic_mean, label='2', color='r', linestyle='dashed',
                  linewidth=2, alpha=0.75  )
    sns.distplot(data, ax=ax[i], bins=20, label='3', color='b')
    ax[i].axvline(np.mean(data), label='4', color='b', linestyle='dashed',
                  linewidth=2, alpha=0.75)
    ax[i].legend(['Fokker Planck solution','Fokker Planck mean','Simulations mean',
                  'Simulations'],fontsize=15,
                 loc='upper left')
    ax[i].set_title('N = {0}'.format(population_sizes[i]), fontsize=20)
    ax[i].set_xlim([-0.02,1])
    ax[i].set_ylim([0,10])
    plt.tight_layout()<jupyter_output>Difference in mean =0.28
Difference in mean =0.27
Difference in mean =0.11<jupyter_text>The plots illustrate how the Fokker-Planck approximation becomes more accurate as $N$ becomes larger. For a better comparison look at the difference in the mean of simulations ( the exact mean ) and the mean obtained by Fokker-Plank distribution. The difference becomes smaller as $N\rightarrow \infty$.<jupyter_code>start = time.time()

population_sizes = [10, 100, 1000, 10000]
mutation_rate = 0
s_list = [0.2, 0.02, 0.002, 0.0002]
w_A = 1
max_steps = 500000
dt = 1
num_sims = 1000
fix_prob_sim = []
fix_prob_Moran = []
fix_prob_wf = []

for i in range(len(population_sizes)):
    print("Population Size: ", population_sizes[i])
    f_M = Fixation_Prob_Moran(population_sizes[i], 1+s_list[i], w_A)
    f_WF = Fixation_Prob_Fisher(population_sizes[i], s_list[i])
    simulations = [Moran_Process(population_sizes[i], 1+s_list[i], w_A,
                                 mutation_rate, mutation_rate, 
                                 1, population_sizes[i] - 1, max_steps,
                                 dt, j) for j in range(num_sims)]
    fix_prob = simulations.count(population_sizes[i])/num_sims
    fix_prob_sim.append(fix_prob)
    fix_prob_Moran.append(f_M)
    fix_prob_wf.append(f_WF)
      
end = time.time()
print('total running time = %1.2f' %(end-start), 'sec') 
plt.scatter(population_sizes, fix_prob_sim, color='b', label='simulation')
plt.scatter(population_sizes, fix_prob_Moran, color='r', label='Moran-Analytical')
plt.scatter(population_sizes, fix_prob_wf, color='orange', label='WF-Analytical')
plt.xscale('Log')
plt.xlabel('population size')
plt.ylabel('fixation probability at N')
plt.ylim([-0.01,0.3])
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1539035305818001, 'fraction_numerical': 0.03891098955743411, 'mean_word_length': 2.970878578479763, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '481196', 'license': 'MIT', 'n_tokens_mistral': 3308, 'n_tokens_neox': 2990, 'n_words': 633, 'path': 'Wright-Fisher-Moran.ipynb', 'repo_name': 'mdaneshv/Stochastic-Simulations', 'size': '8044'}