In [1]:
%pylab
In [2]:
mpl.rcParams.update({'font.family': 'serif',
'font.size': 12,
'font.serif': [],})
In [3]:
from mpl_toolkits.axes_grid1 import ImageGrid
import json
import os.path
$A\Omega = 100$
In [107]:
ddj_2 = np.load('output/ddjd_csweep0_wa100.npz')
with open('params/ddjd_csweep0_wa100.json') as param_file:
params_2 = json.load(param_file)
In [108]:
pdist_params = params_2['pdist_params']
binrange = pdist_params['binrange']
nbins = pdist_params['nbins']
bin_lefts = np.linspace(*binrange, num=nbins, endpoint=False)
bin_ctrs = bin_lefts + (bin_lefts[1] - bin_lefts[0]) / 2.0
In [118]:
rxn_params = params_2['rxn_params']
A_val = rxn_params['k_plus']
B_val = rxn_params['k_minus']
C_vals = ddj_2['C_vals']
j_pdists = ddj_2['j_pdists']
tau = rxn_params['tau_delay']
In [110]:
C_vals
Out[110]:
In [119]:
fig = figure()
ax = fig.add_subplot(111, frameon=False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel('$n(t)$')
ax.set_ylabel('$n(t - \\tau)$')
ax.set_aspect(0.75)
ax.set_position((0.05, 0.05, 0.9, 0.9))
tiles = ImageGrid(fig, 111, nrows_ncols=(3, 4), axes_pad=0.1)
for gridx in range(12):
tiles[gridx].pcolormesh(bin_lefts, bin_lefts, j_pdists[gridx,...].T, cmap='hot')
fig.suptitle('Parameter sweep, $\\Omega A=100$, $B=3$, $C \in [1, 5]$')
Out[119]:
$A\Omega = 10^3$
In [113]:
ddj_3 = np.load('output/ddjd_csweep0_wa1k.npz')
with open('params/ddjd_csweep0_wa1k.json') as param_file:
params_3 = json.load(param_file)
In [114]:
pdist_params = params_3['pdist_params']
binrange = pdist_params['binrange']
nbins = pdist_params['nbins']
bin_lefts = np.linspace(*binrange, num=nbins, endpoint=False)
bin_ctrs = bin_lefts + (bin_lefts[1] - bin_lefts[0]) / 2.0
In [126]:
rxn_params = params_3['rxn_params']
A_val = rxn_params['k_plus']
B_val = rxn_params['k_minus']
C_vals = ddj_3['C_vals']
j_pdists = ddj_3['j_pdists']
tau = rxn_params['tau_delay']
In [121]:
fig = figure()
ax = fig.add_subplot(111, frameon=False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel('$n(t)$')
ax.set_ylabel('$n(t - \\tau)$')
ax.set_aspect(0.75)
ax.set_position((0.05, 0.05, 0.9, 0.9))
tiles = ImageGrid(fig, 111, nrows_ncols=(3, 4), axes_pad=0.1)
for gridx in range(12):
tiles[gridx].pcolormesh(bin_lefts, bin_lefts, j_pdists[gridx,...].T, cmap='hot')
fig.suptitle('Parameter sweep, $\\Omega A=1000$, $B=3$, $C \in [1, 5]$')
Out[121]:
$A\Omega = 10^4$
In [122]:
ddj_4 = np.load('output/ddjd_csweep0_wa10k.npz')
with open('params/ddjd_csweep0_wa10k.json') as param_file:
params_4 = json.load(param_file)
In [123]:
pdist_params = params_4['pdist_params']
binrange = pdist_params['binrange']
nbins = pdist_params['nbins']
bin_lefts = np.linspace(*binrange, num=nbins, endpoint=False)
bin_ctrs = bin_lefts + (bin_lefts[1] - bin_lefts[0]) / 2.0
In [124]:
rxn_params = params_4['rxn_params']
A_val = rxn_params['k_plus']
B_val = rxn_params['k_minus']
C_vals = ddj_4['C_vals']
j_pdists = ddj_4['j_pdists']
tau = rxn_params['tau_delay']
In [125]:
fig = figure()
ax = fig.add_subplot(111, frameon=False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel('$n(t)$')
ax.set_ylabel('$n(t - \\tau)$')
ax.set_aspect(0.75)
ax.set_position((0.05, 0.05, 0.9, 0.9))
tiles = ImageGrid(fig, 111, nrows_ncols=(3, 4), axes_pad=0.1)
for gridx in range(12):
tiles[gridx].pcolormesh(bin_lefts, bin_lefts, j_pdists[gridx,...].T, cmap='hot')
fig.suptitle('Parameter sweep, $\\Omega A=10\,000$, $B=3$, $C \in [1, 5]$')
Out[125]:
$A\Omega = 100$, $B = 4.5$, $\tau = 1.0$
In [170]:
ddj1_1 = np.load('output/ddjd_csweep1_wa100.npz')
with open('params/ddjd_csweep1_wa100.json') as param_file:
params1_1 = json.load(param_file)
In [171]:
pdist_params = params1_1['pdist_params']
binrange = pdist_params['binrange']
nbins = pdist_params['nbins']
bin_lefts = np.linspace(*binrange, num=nbins, endpoint=False)
bin_ctrs = bin_lefts + (bin_lefts[1] - bin_lefts[0]) / 2.0
In [172]:
rxn_params = params1_1['rxn_params']
A_val = rxn_params['k_plus']
B_val = rxn_params['k_minus']
C_vals = ddj1_1['C_vals']
j_pdists = ddj1_1['j_pdists']
tau = rxn_params['tau_delay']
In [279]:
C_vals
Out[279]:
In [173]:
fig = figure()
ax = fig.add_subplot(111, frameon=False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel('$n(t)$')
ax.set_ylabel('$n(t - \\tau)$')
ax.set_aspect(0.4)
ax.set_position((0.05, 0.05, 0.9, 0.9))
tiles = ImageGrid(fig, 111, nrows_ncols=(2, 5), axes_pad=0.1)
for gridx in range(10):
tiles[gridx].pcolormesh(bin_lefts, bin_lefts, j_pdists[gridx + 1,...].T, cmap='hot')
fig.suptitle('Parameter sweep, $\\Omega A=100$, $B=4.5$, $\\tau = 1.0$, $C \in [4, 5.8]$')
Out[173]:
$A\Omega = 400$, $B = 4.5$, $\tau = 1.0$
In [158]:
ddj1_4 = np.load('output/ddjd_csweep1_wa400.npz')
with open('params/ddjd_csweep1_wa400.json') as param_file:
params1_4 = json.load(param_file)
In [159]:
pdist_params = params1_4['pdist_params']
binrange = pdist_params['binrange']
nbins = pdist_params['nbins']
bin_lefts = np.linspace(*binrange, num=nbins, endpoint=False)
bin_ctrs = bin_lefts + (bin_lefts[1] - bin_lefts[0]) / 2.0
In [160]:
rxn_params = params1_4['rxn_params']
A_val = rxn_params['k_plus']
B_val = rxn_params['k_minus']
C_vals = ddj1_4['C_vals']
j_pdists = ddj1_4['j_pdists']
tau = rxn_params['tau_delay']
In [161]:
fig = figure()
ax = fig.add_subplot(111, frameon=False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel('$n(t)$')
ax.set_ylabel('$n(t - \\tau)$')
ax.set_aspect(0.4)
ax.set_position((0.05, 0.05, 0.9, 0.9))
tiles = ImageGrid(fig, 111, nrows_ncols=(2, 5), axes_pad=0.1)
for gridx in range(10):
tiles[gridx].pcolormesh(bin_lefts, bin_lefts, j_pdists[gridx,...].T, cmap='hot')
fig.suptitle('Parameter sweep, $\\Omega A=400$, $B=4.5$, $\\tau = 1.0$, $C \in [4, 5.8]$')
Out[161]:
In [174]:
margdist_n = np.sum(j_pdists, axis=2)
# The division here takes care of any normalization factors on the joint distribution
cond_avg = np.sum(j_pdists * bin_ctrs[np.newaxis,np.newaxis,:] / margdist_n[:,:,np.newaxis], axis=2)
In [147]:
cond_avg.shape
Out[147]:
The Langevin limit: a straight line
In [175]:
C_idx = 8
C_val = C_vals[C_idx]
omega_h = np.sqrt(C_val**2 - B_val**2)
ca_slope = (omega_h * np.cos(omega_h * tau) - B_val * np.sin(omega_h * tau)) / (omega_h + C_val * np.sin(omega_h * tau))
n_star = A_val / (B_val + C_val)
In [176]:
omega_h
Out[176]:
In [168]:
fig = figure()
plot((bin_ctrs - n_star) / np.sqrt(A_val), (cond_avg[C_idx,...] - n_star) / np.sqrt(A_val), figure=fig)
xi_grid = np.linspace(-8, 8, 100)
plot(xi_grid, ca_slope * xi_grid)
print(C_val)
In [177]:
fig = figure()
ax = fig.gca()
im = ax.pcolormesh(bin_lefts, bin_lefts, j_pdists[C_idx,...].T, cmap='hot')
ax.plot(bin_ctrs, cond_avg[C_idx,...], color='green', linewidth=2.0)
ax.plot((xi_grid * np.sqrt(A_val)) + n_star, (ca_slope * xi_grid * np.sqrt(A_val)) + n_star, color='blue', linewidth=2.0)
ax.set_ylim(min(bin_lefts), max(bin_lefts))
ax.set_xlim(min(bin_lefts), max(bin_lefts))
Out[177]:
In [157]:
plot(bin_ctrs, cond_avg.T)
legend(['$C = {:.2g}$'.format(C_val) for C_val in C_vals[:12]])
Out[157]:
For the delayed-deg system below the Hopf bifurcation:
Analytical distribution is $P(\xi) = \sqrt{\frac{b + c}{2\pi (1 + c\tau)}}\exp\left(- \frac{b + c}{2 (1 + c\tau)}\xi^2 \right)$ (as a function of fluctuation $\xi$ from the mean of $x^\star = \frac{a}{b+c}$)
$A\Omega = 100$, $\tau = 20$
In [91]:
ddwe1_res = np.load('output/ddwe_gauss_t20_wa100_res.npz')
ddwe1_nores = np.load('output/ddwe_gauss_t20_wa100_nores.npz')
In [92]:
param_fname = str(ddwe1_res['param_fname'])
with open(param_fname, 'r') as param_file:
pdwe_params = json.load(param_file)
rxn_list = pdwe_params['rxn_list']
prmA = rxn_list[0]['propensity_const']
prmB = rxn_list[1]['propensity_const']
prmC = rxn_list[2]['propensity_const']
prmTau = rxn_list[2]['delay']
In [93]:
pdists_res = ddwe1_res['prob_dists']
pdists_nores = ddwe1_nores['prob_dists']
bincts_res = ddwe1_res['bin_counts']
bincts_nores = ddwe1_nores['bin_counts']
bin_xs = ddwe1_res['bin_xs']
nens = pdists_res.shape[0]
$A\Omega = 200$
In [477]:
ddwe2_res = np.load('output/ddwe_gauss_wa200_freqres.npz')
ddwe2_nores = np.load('output/ddwe_gauss_wa200_freqnores.npz')
In [478]:
rxn_params = ddwe2_res['rxn_params'][()]
print(rxn_params)
prmA = rxn_params['k_plus']
prmB = rxn_params['k_minus']
prmC = rxn_params['k_delayed']
prmTau = rxn_params['tau_delay']
In [479]:
pdists_res = ddwe2_res['prob_dists']
pdists_nores = ddwe2_nores['prob_dists']
bincts_res = ddwe2_res['bin_counts']
bincts_nores = ddwe2_nores['bin_counts']
bin_xs = ddwe2_res['bin_xs']
nens = pdists_res.shape[0]
$A\Omega = 400$
In [443]:
ddwe4_res = np.load('output/ddwe_gauss_wa400_freqres.npz')
ddwe4_nores = np.load('output/ddwe_gauss_wa400_freqnores.npz')
In [444]:
rxn_params = ddwe4_res['rxn_params'][()]
print(rxn_params)
prmA = rxn_params['k_plus']
prmB = rxn_params['k_minus']
prmC = rxn_params['k_delayed']
prmTau = rxn_params['tau_delay']
In [445]:
pdists_res = ddwe4_res['prob_dists']
pdists_nores = ddwe4_nores['prob_dists']
bincts_res = ddwe4_res['bin_counts']
bincts_nores = ddwe4_nores['bin_counts']
bin_xs = ddwe4_res['bin_xs']
nens = pdists_res.shape[0]
$A\Omega = 800$
In [54]:
ddwe8_res = np.load('output/ddwe_gauss_wa800_freqres.npz')
ddwe8_nores = np.load('output/ddwe_gauss_wa800_freqnores.npz')
In [55]:
rxn_params = ddwe8_res['rxn_params'][()]
print(rxn_params)
prmA = rxn_params['k_plus']
prmB = rxn_params['k_minus']
prmC = rxn_params['k_delayed']
prmTau = rxn_params['tau_delay']
In [491]:
pdists_res = ddwe8_res['prob_dists']
pdists_nores = ddwe8_nores['prob_dists']
bincts_res = ddwe8_res['bin_counts']
bincts_nores = ddwe8_nores['bin_counts']
bin_xs = ddwe8_res['bin_xs']
nens = pdists_res.shape[0]
$A\Omega = 100$, $\tau = 1$
In [98]:
ddwe11_res = np.load('output/ddwe_gauss_t1_wa100_res.npz')
ddwe11_nores = np.load('output/ddwe_gauss_t1_wa100_nores.npz')
In [99]:
param_fname = str(ddwe11_res['param_fname'])
with open(param_fname, 'r') as param_file:
pdwe_params = json.load(param_file)
rxn_list = pdwe_params['rxn_list']
prmA = rxn_list[0]['propensity_const']
prmB = rxn_list[1]['propensity_const']
prmC = rxn_list[2]['propensity_const']
prmTau = rxn_list[2]['delay']
In [100]:
pdists_res = ddwe11_res['prob_dists']
pdists_nores = ddwe11_nores['prob_dists']
bincts_res = ddwe11_res['bin_counts']
bincts_nores = ddwe11_nores['bin_counts']
bin_xs = ddwe11_res['bin_xs']
nens = pdists_res.shape[0]
$A\Omega = 400$, $\tau = 1$
In [28]:
ddwe14_res = np.load('output/ddwe_gauss_t1_wa400_res.npz')
ddwe14_nores = np.load('output/ddwe_gauss_t1_wa400_nores.npz')
In [29]:
param_fname = str(ddwe14_res['param_fname'])
with open(param_fname, 'r') as param_file:
pdwe_params = json.load(param_file)
rxn_list = pdwe_params['rxn_list']
prmA = rxn_list[0]['propensity_const']
prmB = rxn_list[1]['propensity_const']
prmC = rxn_list[2]['propensity_const']
prmTau = rxn_list[2]['delay']
In [30]:
pdists_res = ddwe14_res['prob_dists']
pdists_nores = ddwe14_nores['prob_dists']
bincts_res = ddwe14_res['bin_counts']
bincts_nores = ddwe14_nores['bin_counts']
bin_xs = ddwe14_res['bin_xs']
nens = pdists_res.shape[0]
$A\Omega = 800$, $\tau = 1$
In [9]:
ddwe18_res = np.load('output/ddwe_gauss_t1_wa800_res.npz')
ddwe18_nores = np.load('output/ddwe_gauss_t1_wa800_nores.npz')
In [10]:
param_fname = str(ddwe18_res['param_fname'])
with open(param_fname, 'r') as param_file:
pdwe_params = json.load(param_file)
rxn_list = pdwe_params['rxn_list']
prmA = rxn_list[0]['propensity_const']
prmB = rxn_list[1]['propensity_const']
prmC = rxn_list[2]['propensity_const']
prmTau = rxn_list[2]['delay']
In [11]:
pdists_res = ddwe18_res['prob_dists']
pdists_nores = ddwe18_nores['prob_dists']
bincts_res = ddwe18_res['bin_counts']
bincts_nores = ddwe18_nores['bin_counts']
bin_xs = ddwe18_res['bin_xs']
nens = pdists_res.shape[0]
In [12]:
np.sum(bincts_nores, axis=1)
Out[12]:
$A\Omega = 100$
In [4]:
pdwe1_res = np.load('output/pdwe_wa100_res.npz')
pdwe1_nores = np.load('output/pdwe_wa100_nores.npz')
In [5]:
param_fname = str(pdwe1_res['param_fname'])
with open(param_fname, 'r') as param_file:
pdwe_params = json.load(param_file)
rxn_list = pdwe_params['rxn_list']
prmA = rxn_list[0]['propensity_const']
prmB = rxn_list[1]['propensity_const']
In [6]:
pdists_res = pdwe1_res['prob_dists']
pdists_nores = pdwe1_nores['prob_dists']
bincts_res = pdwe1_res['bin_counts']
bincts_nores = pdwe1_nores['bin_counts']
bin_xs = pdwe1_res['bin_xs']
nens = pdists_res.shape[0]
In [8]:
np.sum(bincts_nores, axis=1)
Out[8]:
$A\Omega = 400$
In [24]:
pdwe4_res = np.load('output/pdwe_wa400_res.npz')
pdwe4_nores = np.load('output/pdwe_wa400_nores.npz')
In [25]:
param_fname = str(pdwe4_res['param_fname'])
with open(param_fname, 'r') as param_file:
pdwe_params = json.load(param_file)
rxn_list = pdwe_params['rxn_list']
prmA = rxn_list[0]['propensity_const']
prmB = rxn_list[1]['propensity_const']
In [26]:
pdists_res = pdwe4_res['prob_dists']
pdists_nores = pdwe4_nores['prob_dists']
bincts_res = pdwe4_res['bin_counts']
bincts_nores = pdwe4_nores['bin_counts']
bin_xs = pdwe4_res['bin_xs']
nens = pdists_res.shape[0]
Using $P(\xi) = \sqrt{\frac{B + C}{2\pi (1 + C\tau)}}\exp\left(- \frac{B + C}{2 (1 + C\tau)}\xi^2 \right)$ (as a function of normalized fluctuation $\xi$ from the mean of $n^\star = \frac{A}{B+C}$; $\xi = \frac{n - n^\star}{\sqrt{A}}$)
In [181]:
an_normfact = np.sqrt((prmB + prmC) / (2*np.pi * (1 + prmC*prmTau)))
an_sigma = np.sqrt((1 + prmC*prmTau) / (prmB + prmC))
xi_max = 3 * an_sigma
xi_grid = np.linspace(-1*xi_max, xi_max, 1000)
an_dist = an_normfact * np.exp((prmB + prmC) * -1 * xi_grid**2 / (2 * (1 + prmC*prmTau)))
In [182]:
n_star = prmA / (prmB + prmC)
bin_xis = (bin_xs - n_star) / np.sqrt(prmA)
pdist_scale = np.sqrt(prmA)
discrete_shift = 0.5 * (bin_xs[1] - bin_xs[0] - 1.0)
bin_xis_shifted = (bin_xs + discrete_shift - n_star) / pdist_scale
an_coarse = np.sqrt(prmB / (2 * np.pi)) * np.exp(-1.0 * prmB * bin_xis_shifted**2 / 2.0)
In [89]:
bin_xis_shifted
Out[89]:
In [64]:
an_sigma * np.sqrt(prmA)
Out[64]:
In [65]:
an_sigma_noct = np.sqrt(1.0 / (prmB + prmC))
an_sigma_noct * np.sqrt(prmA)
Out[65]:
In [66]:
np.trapz(an_dist, x=xi_grid)
Out[66]:
$P(n) = \sqrt{\frac{B}{2\pi A}}\exp\left(-\frac{B}{2A} \left(n - \frac{A}{B}\right)^2\right)$ (again, using $\Omega = 1$ and writing everything in terms of $A$).
In [12]:
n_grid = np.linspace(min(bin_xs), max(bin_xs), 1000)
pdist_analytic = sqrt(prmB / (2.0 * pi * prmA)) * exp(-1.0 * prmB / (2.0 * prmA) * (n_grid - prmA * 1.0 / prmB)**2)
In [13]:
np.trapz(pdist_analytic, n_grid)
Out[13]:
Or, in terms of $\xi = A^{-1/2}\left(n - \frac{A}{B}\right)$: $P(\xi) = \sqrt{\frac{B}{2\pi}}\exp\left(-\frac{B\xi^2}{2}\right)$:
In [66]:
an_sigma = 1.0 / np.sqrt(prmB)
xi_max = 4 * an_sigma
xi_grid = np.linspace(-1*xi_max, xi_max, 1000)
an_dist = np.sqrt(prmB / (2 * np.pi)) * np.exp(-1.0 * prmB * xi_grid**2 / 2.0)
pdist_scale = np.sqrt(prmA)
n_star = prmA / prmB
bin_xis = (bin_xs - n_star) / pdist_scale
discrete_shift = 0.5 * (bin_xs[1] - bin_xs[0] - 1.0)
bin_xis_shifted = (bin_xs + discrete_shift - n_star) / pdist_scale
an_coarse = np.sqrt(prmB / (2 * np.pi)) * np.exp(-1.0 * prmB * bin_xis_shifted**2 / 2.0)
In [48]:
an_sigma * np.sqrt(prmA)
Out[48]:
In [49]:
n_star - 45, n_star + 45
Out[49]:
In [27]:
pd_mean_res = np.mean(pdists_res, axis=0)
pd_mean_nores = np.mean(pdists_nores, axis=0)
pd_sem_res = np.std(pdists_res, axis=0) / np.sqrt(nens)
pd_sem_nores = np.std(pdists_nores, axis=0) / np.sqrt(nens)
In [68]:
fig = figure()
ax = fig.gca()
im = ax.pcolormesh(bin_xs, array(range(nens)), pdists_res, cmap='hot')
ax.set_xlabel('Number of X')
ax.set_ylabel('Ensemble Index')
ax.set_title('Probability Distribution, Resampling')
fig.colorbar(im)
Out[68]:
In [13]:
fig = figure()
ax = fig.gca()
im = ax.pcolormesh(bin_xs, array(range(nens)), bincts_res, cmap='hot')
ax.set_ylabel('Ensemble Index')
ax.set_xlabel('Number of X')
ax.set_title('Bin Counts, Resampling')
fig.colorbar(im)
Out[13]:
In [40]:
fig = figure()
ax = fig.gca()
im = ax.pcolormesh(bin_xs, array(range(nens)), pdists_nores, cmap='hot')
ax.set_xlabel('Number of X')
ax.set_ylabel('Ensemble Index')
ax.set_title('Probability Distribution, No Resampling')
fig.colorbar(im)
Out[40]:
In [14]:
fig = figure()
ax = fig.gca()
im = ax.pcolormesh(bin_xs, array(range(nens)), bincts_nores, cmap='hot')
ax.set_ylabel('Ensemble Index')
ax.set_xlabel('Number of X')
ax.set_title('Bin Counts, No Resampling')
fig.colorbar(im)
Out[14]:
In [39]:
pd_sigmas = dict()
In [71]:
pd_mean_ctr = np.sum(bin_xs * pd_mean_res)
pd_mean_sigma = np.sqrt(np.sum(bin_xs**2 * pd_mean_res) - pd_mean_ctr**2)
print(pd_mean_ctr)
print(prmA / (prmB + prmC))
print(pd_mean_sigma)
In [72]:
pd_sigmas[prmA] = pd_mean_sigma
In [73]:
pd_sigmas
Out[73]:
In [81]:
A_sigmas = array(list(pd_sigmas.items()))
plot(log(A_sigmas[:,0]), log(A_sigmas[:,1]), 'o', label='Simulation')
A_1 = 100
sigma_1 = pd_sigmas[100]
A_grid = np.logspace(0.9, 1.5, 1000, base=A_1)
plot(log(A_grid), log(sigma_1) + 0.5 * (log(A_grid) - log(A_1)), label=r'$\sigma \propto \sqrt{A}$')
legend(loc='best')
xlabel(r'$\ln(A)$')
ylabel(r'$\ln(\sigma)$')
Out[81]:
In [76]:
log(800) / log(100)
Out[76]:
In [30]:
fig = figure()
ax = fig.gca()
bar_width = bin_xs[1] - bin_xs[0]
ax.bar(bin_xs, pd_mean_res, yerr=pd_sem_res, width=bar_width, color='blue', alpha=0.6, ecolor='blue')
ax.bar(bin_xs, pd_mean_nores, yerr=pd_sem_nores, width=bar_width, color='red', alpha=0.6, ecolor='red')
ax.set_ylim(bottom=0.0)
ax.legend(('With Resampling', 'No Resampling'))
ax.set_xlabel('Number $n$')
ax.set_ylabel('Probability $P(n)$')
#ax.set_title(r'Delayed-degradation system, $\Omega A = {:g}$, $B=3$, $C=1$, $\tau=20$'.format(prmA))
Out[30]:
In terms of $\xi = \frac{n - n^\star}{\sqrt{A}}$:
In [29]:
fig = figure()
ax = fig.gca()
bar_width = bin_xis[1] - bin_xis[0]
bar_scale = 1.0 / bar_width
ax.bar(bin_xis_shifted - 0.5*bar_width, pd_mean_res * bar_scale, yerr=pd_sem_res * bar_scale,
width=bar_width, color='blue', alpha=0.6, ecolor='blue', label='With Resampling')
ax.bar(bin_xis_shifted - 0.5*bar_width, pd_mean_nores * bar_scale, yerr=pd_sem_nores * bar_scale,
width=bar_width, color='red', alpha=0.6, ecolor='red', label='No Resampling')
ax.set_ylim(bottom=0.0)
ax.set_xlim(-3.0, 3.0)
ax.plot(xi_grid, an_dist, color='black', label='Analytical')
ax.legend(loc='best')
ax.set_xlabel(r'Fluctuation $\xi = A^{-1/2}(n - n^\star)$')
ax.set_ylabel('Probability $P(n)$')
#ax.set_title(r'Delayed-degradation system, $A = {:g}$, $B={:g}$, $C={:g}$, $\tau={:g}$'.format(prmA, prmB, prmC, prmTau))
ax.set_title(r'Production-Degradation system, $A = {:g}$, $B={:g}$'.format(prmA, prmB))
Only without resampling:
In [52]:
fig = figure()
ax = fig.gca()
bar_width = bin_xis[1] - bin_xis[0]
bar_scale = 1.0 / bar_width
#ax.bar(bin_xis_shifted - 0.5*bar_width, pd_mean_res * bar_scale, yerr=pd_sem_res * bar_scale,
# width=bar_width, color='blue', alpha=0.6, ecolor='blue', label='With Resampling')
ax.bar(bin_xis_shifted - 0.5*bar_width, pd_mean_nores * bar_scale, yerr=pd_sem_nores * bar_scale,
width=bar_width, color='red', alpha=0.6, ecolor='red', label='No Resampling')
ax.set_ylim(bottom=0.0)
ax.plot(xi_grid, an_dist, color='black', label='Analytical')
ax.legend(loc='best')
ax.set_xlabel(r'Fluctuation $\xi = A^{-1/2}(n - n^\star)$')
ax.set_ylabel('Probability $P(n)$')
ax.set_title(r'Delayed-degradation system, $A = {:g}$, $B={:g}$, $C={:g}$, $\tau={:g}$'.format(prmA, prmB, prmC, prmTau))
#ax.set_title(r'Production-Degradation system, $A = {:g}$, $B={:g}$'.format(prmA, prmB))
Out[52]:
And now a plot of the deviations from analytical
In [70]:
fig = figure()
ax = fig.gca()
bar_width = bin_xis[1] - bin_xis[0]
bar_scale = 1.0 / bar_width
devs_res = (pd_mean_res * bar_scale - an_coarse) / (pd_sem_res * bar_scale)
devs_nores = (pd_mean_nores * bar_scale - an_coarse) / (pd_sem_nores * bar_scale)
keep_range = np.abs(bin_xis) <= 2.0
ax.bar(bin_xis[keep_range], devs_res[keep_range], width=bar_width, color='blue', alpha=0.6, label='With Resampling')
ax.bar(bin_xis[keep_range], devs_nores[keep_range], width=bar_width, color='red', alpha=0.6, label='No Resampling')
ax.grid()
ax.axhline(0.0, color='black', linewidth=1.5)
ax.legend(loc='best')
ax.set_xlabel(r'Fluctuation $\xi = A^{-1/2}(n - n^\star)$')
ax.set_ylabel(r'Normalized Difference $(P - P_\mathrm{an})/\sigma_\mathrm{est}$')
ax.set_title('Difference from Analytical Distribution')
Out[70]:
In [71]:
fig = figure()
ax = fig.gca()
bar_width = bin_xis[1] - bin_xis[0]
bar_scale = 1.0 / bar_width
devs_rescmp = (pd_mean_res - pd_mean_nores) / sqrt(pd_sem_res**2 + pd_sem_nores**2)
keep_range = np.abs(bin_xis) <= 2.0
ax.bar(bin_xis[keep_range], devs_rescmp[keep_range], width=bar_width, color='blue')
ax.grid()
ax.axhline(0.0, color='black', linewidth=1.5)
ax.legend(loc='best')
ax.set_xlabel(r'Fluctuation $\xi = A^{-1/2}(n - n^\star)$')
ax.set_ylabel(r'Normalized Difference $(P_\mathrm{res} - P_\mathrm{non})/\sqrt{\sigma_\mathrm{res}^2 + \sigma_\mathrm{non}^2}$')
ax.set_title('With vs. Without Resampling')
Out[71]:
In [ ]: