Skip to content

Instantly share code, notes, and snippets.

@naught101
Created August 26, 2015 06:42
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save naught101/d3e28b73b87d6d69c0c4 to your computer and use it in GitHub Desktop.
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