Toy rank histograms

toy_rank_histograms.py

# -*- 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))