add diagnostic.jl

src/DataAssim.jl

@@ -18,6 +18,7 @@ include("ensemble.jl")
 include("fourDVar.jl")
 include("KalmanFilter.jl")
 include("TwinExperiment.jl")
+include("diagnostic.jl")
 
 export FreeRun, fourDVar, TwinExperiment, LinShallowWater1DModel, KalmanFilter
 export AbstractModel, ModelFun

src/diagnostic.jl

@@ -0,0 +1,47 @@
+
+
+function _talagrand_diagram!(nbins,buffer,x,y)
+    @assert size(x,1) == size(y,1)
+    @assert size(x,2) == size(buffer,1)
+    @assert size(nbins,1) == size(buffer,1) + 1
+
+    @inbounds for i = 1:size(x,1)
+        buffer .= @view x[i,:]
+        sort!(buffer)
+        ibin = count(<(y[i]), buffer) + 1
+        nbins[ibin] += 1
+    end
+
+    return nbins
+end
+
+_talagrand_diagram(x,y) = _talagrand_diagram!(zeros(Int,size(x,2)+1),x[1,:],x,y)
+
+"""
+    freq = talagrand_diagram(x,y)
+
+
+Compute the frequency of each bins for a Talagrand Diagram where
+`x` is the ensemble (array of the size sample × ensemble size) and `y` are
+the obervations (sample).
+
+## Example
+
+```julia
+using PyPlot
+Nens = 100
+Nsample = 10000
+# Ensemble
+x = randn(Nsample,Nens)
+# Observation
+y = randn(Nsample)
+
+freq = talagrand_diagram(x,y)
+plot(freq)
+plot(ones(size(freq)) / size(freq,1))
+ylim(0,ylim()[2])
+xlabel("bins")
+ylabel("frequency")
+```
+"""
+talagrand_diagram(x,y) = _talagrand_diagram(x,y) / size(x,1)

test/runtests.jl

@@ -5,3 +5,4 @@ using Statistics
 
 include("test_ensemble.jl")
 include("test_4dvar_assim.jl")
+include("test_diagnostic.jl")

test/test_diagnostic.jl

@@ -0,0 +1,16 @@
+using Random
+using Test
+using DataAssim: talagrand_diagram
+
+# Example from
+# ABC of ensemble forecasting
+# by Peter Houtekamer, ARMA
+# https://collaboration.cmc.ec.gc.ca/cmc/cmoi/product_guide/docs/lib/ens_en.pdf
+# https://web.archive.org/web/20231118202650/https://collaboration.cmc.ec.gc.ca/cmc/cmoi/product_guide/docs/lib/ens_en.pdf
+
+x = [1.5;; 2.3;; 0.8;; 4.1;; 0.3]
+y = [2.95]
+
+
+freq = talagrand_diagram(x,y)
+@test nbins == [0,0,0,0,1,0]