<jupyter_start><jupyter_text>Time series analysis<jupyter_code>lon = -129
lat= -55

p0_ = p0.sel(longitude=lon, latitude=lat, method='nearest')
mu0_ = mu0.sel(longitude=lon, latitude=lat, method='nearest')
mu1_ = mu1.sel(longitude=lon, latitude=lat, method='nearest')
beta0_ = beta0.sel(longitude=lon, latitude=lat, method='nearest')
beta1_ = beta1.sel(longitude=lon, latitude=lat, method='nearest')
true_ = true.sel(longitude=lon, latitude=lat, method='nearest')

plt.figure()
plt.plot(p0_)
plt.plot(true_, '-x')
plt.plot(mu0_, 'r')
plt.plot(mu1_, 'k')
plt.plot(mu0_ -1.96 * 1 / beta0_, '--r')
plt.plot(mu0_ + 1.96 * 1 / beta0_, '--r')
plt.plot(mu1_ - 1.96 * 1 / beta1_, '--k')
plt.plot(mu1_ + 1.96 * 1 / beta1_, '--k')
time = 796
n_samples = 1000000
sel = np.random.rand(n_samples) > float(p0_.isel(time=time))
sel = np.arange(n_samples) + sel * n_samples
samples0 = np.random.randn(n_samples) / float(beta0_.isel(time=time)) + float(mu0_.isel(time=time))
samples1 = np.random.randn(n_samples) / float(beta1_.isel(time=time)) + float(mu1_.isel(time=time))
samples = np.concatenate((samples0, samples1))
final_samples = samples[sel]
plt.figure()
print(float(true_.isel(time=time)))
_ = plt.hist(final_samples, bins=200)<jupyter_output><empty_output><jupyter_text>Global distribution of true and simulated forcing We assess the goodness of fit  Likelihood map<jupyter_code>from scipy.stats import norm
def pdf(x, p0, mu0, beta0, mu1, beta1):
    return p0 * norm.pdf((x - mu0) * beta0) + (1 - p0) * norm.pdf((x - mu1) * beta1)

lkh = pdf(true, p0, mu0, beta0, mu1, beta1)
plt.figure()
plt.imshow(lkh.mean(dim='time'), origin='lower', vmin=0, vmax=0.5)<jupyter_output><empty_output><jupyter_text>Mean likelihood<jupyter_code>lat = slice(-40, 40)
with ProgressBar():
    mean_lkh = apply_complete_mask(lkh.sel(latitude=lat)).mean().compute()
mean_lkh<jupyter_output><empty_output><jupyter_text>Goodness of fit<jupyter_code>from scipy.stats import norm
def my_transform(x ,p0, mu0, beta0, mu1, beta1):
    cdf = lambda x: p0 * norm.cdf((x - mu0) * beta0) + (1 - p0) * norm.cdf((x - mu1) * beta1)
    return cdf(x)

v = my_transform(true, p0, mu0, beta0, mu1, beta1)
v = v.sel(latitude=slice(-40, 40))
v = apply_complete_mask(v)
v
plt.figure()
_ = plt.hist(v.values.flatten(), bins=200, density=True)
v.shape
quantiles = np.exp(np.linspace(-10, 10, 100)) / (1 + np.exp(np.linspace(-10, 10, 100)))

q = np.nanquantile(v, quantiles)
plt.figure()
plt.plot(quantiles, q, 'x')
plt.plot(quantiles, quantiles)
norm.ppf(0.001)
from scipy.stats import norm
s = norm.ppf(q)
plt.figure()
plt.plot(norm.ppf(quantiles),  s, 'x')
plt.plot(norm.ppf(quantiles),  norm.ppf(quantiles))
plt.axis([None, None, -18, None])<jupyter_output><empty_output><jupyter_text>Distribution of simulated vs true forcing<jupyter_code>saved_true = true
def apply_complete_mask(array):
    mask = uv_plotter.borders
    mask2 = uv_plotter.mask
    mask = mask.interp({k: array.coords[k] for k in ['longitude', 'latitude']})
    mask2 = mask2.interp({k: array.coords[k] for k in ['longitude', 'latitude']})
    array = array.where(np.isnan(mask) & (~np.isnan(mask2)))
    array = array.sel(latitude=slice(pred['latitude'][0], pred['latitude'][-1]))
    return array

true = apply_complete_mask(saved_true)

epsilon = np.random.randn(*true.shape)
epsilon2 = np.random.randn(*true.shape)
bernouilli = np.random.rand(*true.shape) > p0
simulated = bernouilli * (mu0 + epsilon / beta0) + (1 - bernouilli) * (mu1 + epsilon / beta1)

simulated = apply_complete_mask(simulated)

plt.figure()
_ = plt.hist(true.values.ravel(), log=True, bins=np.arange(-20, 21, 0.5), density=True)
_ = plt.hist(simulated.values.ravel(), log=True, bins=np.arange(-20, 21, 0.5), density=True)

quantiles = np.exp(np.linspace(-3, 3, 100)) / (1 + np.exp(np.linspace(-3, 3, 100)))
quantiles
q_true = np.nanquantile(true.values.ravel(), quantiles)
q_simu = np.nanquantile(simulated.values.ravel(), quantiles)
plt.figure()
plt.plot(q_true, q_simu, 'x')
plt.plot(q_true, q_true)
plt.figure()
plt.plot(q_true, q_simu, 'x')
plt.plot(q_true, q_true, 'o')
simulated.shape<jupyter_output><empty_output><jupyter_text>Quantiles analysis<jupyter_code>import matplotlib
import cmocean
mean_pred = p0 * mu0 + p1 * mu1
plt.figure()
plt.imshow(np.abs((mean_pred.mean(dim='time') - true.mean(dim='time'))) / true.std(dim='time'), vmin=0.01, vmax=1, norm=matplotlib.colors.LogNorm()
, origin='lower', cmap=cmocean.cm.delta)
plt.colorbar()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17950435365036838, 'fraction_numerical': 0.03951774949765573, 'mean_word_length': 5.676602086438152, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 27, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '443663', 'license': 'MIT', 'n_tokens_mistral': 1951, 'n_tokens_neox': 1833, 'n_words': 336, 'path': 'Untitled1.ipynb', 'repo_name': 'm2lines/subgrid', 'size': '4479'}