Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Toy rank histograms
# -*- coding: utf-8 -*-
"""
Created on Wed Aug 26 12:12:37 2015
@author: naught101
"""
import numpy as np
import matplotlib.pyplot as pl
year = 30
nyears = 9
nmodels = 10
trend=0.5
ensemble = np.random.uniform(high=trend*nyears + 1, size=(year*nyears, nmodels))
obs = np.random.uniform(size=(year*nyears)) + \
np.linspace(0, trend*nyears, num=year*nyears)
obs = np.reshape(obs, (year*nyears,1))
combined = np.concatenate((obs, ensemble), axis=1)
obs_ranks = np.apply_along_axis(lambda x: x.argsort().argsort()[0], 1, combined)
# timeseries
pl.subplot(5,3,(1,3))
pl.plot(ensemble[1:1000,], alpha=0.1)
pl.plot(obs[1:1000,])
pl.title("{0} models; {1} 'years' of {2} timesteps. Obs with reduced variance + trend".format(
nmodels, nyears, year))
# full RH
pl.subplot(5,3,4)
pl.hist(obs_ranks)
pl.title("full rank histogram")
# Overplotted
for i in range(nyears):
pl.subplot(5,3,5)
pl.hist(obs_ranks[(i*year):((i+1)*year),], range=(0,nmodels), alpha=0.1)
pl.title("annual histograms overplotted")
# boxplots
ranks = np.full((nyears,nmodels), np.nan)
for i in range(nyears):
ranks[i,] = pl.hist(obs_ranks[(i*year):((i+1)*year),], range=(0,nmodels), alpha=0.1)[0]
pl.subplot(5,3,6)
pl.boxplot(ranks, positions=list(range(nmodels)))
pl.plot(ranks.T, alpha=0.5)
pl.title("annual histograms averaged")
# individual years
for i in range(nyears):
pl.subplot(5,3,i+7)
pl.hist(obs_ranks[(i*year):((i+1)*year),], range=(0,nmodels))
pl.title('year ' + str(i))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment