Different results obtained with CRAN/verification package for CRPS decomposition

[vinfort]

I've compared the result of the crps decomposition to those obtained with the CRAN verification package and the numbers do not match.

For some test data, here is what I obtain for CRPS, CRPS_Reli and CRPS_Pot:

with CRAN/verification package: [3.1440921604902003, 0.515740759959203, 2.628351400530997]

with ensverif package, using HEAD of main branch: [3.1440921604901995, 1.3397283115077352, 1.8043638489824645]

with ensverif package, after applying the patch proposed in issue 1 to fix a problem with other test data: [3.1440921604901995, 1.2974725396334514, 1.8466196208567482]

The total CRPS matches, but not the decomposition. The correction coming from issue 1 is minor compared to the differences with the CRAN package.

I would advise against using this function in ensverif until the issue is resolved.

[vinfort]

Here is the test data: Forecasts:

,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20
2023-09-21,41.527937235682984,47.564373423486124,51.09494434858676,53.811769699237814,56.12117823457851,58.19342667113706,60.11936210488422,61.95489180872106,63.74019718874394,65.5070689665543,67.2838812090359,69.09887281866106,70.98359073251717,72.97689147207892,75.13175091902444,77.52647678634591,80.29357968130151,83.69040874535996,88.35297109310373,97.06390687297073
2023-09-22,57.72709464205502,65.26648544366056,69.65510484601847,73.02282973835766,75.87966393397002,78.43881760946519,80.81381942312693,83.07438807308336,85.27046567136007,87.44141548064378,89.62221720457707,91.84753052142645,94.15592046942182,96.59473365272939,99.2283646687571,102.15183296182752,105.52574767451509,109.6617207223748,115.32923126392387,125.89014372875299
2023-09-23,19.321727203246073,23.79595686164574,26.482912964631083,28.58148352100599,30.384847692182156,32.0173089250475,33.54597181369877,35.01266320430217,36.447976845581444,37.87661317959488,39.321125535521,40.80447173340417,42.3528126528601,43.99888232463375,45.7878206112342,47.786921842423226,50.11066206234185,52.98236510799713,56.956360972987554,64.47134651667623
2023-09-24,5.088600383550843,7.5019399347084486,9.044505399415442,10.289430345205314,11.384143281844716,12.393176103481567,13.352418589938887,14.284943439894684,15.208333049947251,16.137424809878844,17.086411244347467,18.070374716106617,19.107125974159516,20.21955976443203,21.439847625159363,22.816605062602733,24.43323346984912,26.453578088751893,29.286991243127595,34.75008169348809
2023-09-25,1.8323848475237756,3.35192054851929,4.392918192790584,5.2612654230500135,6.0417775375151805,6.773205033043798,7.477917083187808,8.170823889447673,8.863804079281792,9.567332309926273,10.29186104284194,11.048905304021279,11.852446379494202,12.720817208787889,13.680125136906018,14.770192144775372,16.05972509006288,17.68433408509734,19.984283239688597,24.478225379252656
2023-09-26,0.8117339057144826,1.861769846096195,2.630806555333289,3.291054455550672,3.8954117086985334,4.4693153130584395,5.028060064796095,5.582221886806186,6.1405849302581945,6.711194276137087,7.30235315087387,7.92345691110563,8.586138372458912,9.30585721506104,10.104822917880195,11.017116497885711,12.101754255674516,13.475564129873842,15.432460614456588,19.288971546730508
2023-09-27,0.5309839508508899,1.39752495553691,2.056669897675647,2.6312240819020096,3.1620244546854246,3.66940170803442,4.165897478737873,4.660374054553755,5.160364429185901,5.672905663508547,6.2053827137273165,6.766258303079494,7.366104543549275,8.019055911062326,8.745496405737468,9.576787933493629,10.567323306530367,11.824914560466192,13.621098169998113,17.174064906828438
2023-09-28,0.45543927131670076,1.2513931407080734,1.863803677313699,2.3999904617833483,2.8966611558971467,3.372302405285338,3.838411349359157,4.303166567620608,4.773568091730482,5.256192807169292,5.75797548109955,6.286891007287524,6.852926560430661,7.469456369550556,8.155788133946483,8.941648761155102,9.87861716242262,11.068962706499665,12.770332272344367,16.13913805285335
2023-09-29,0.7970557244911081,1.7546245190542569,2.447769725793726,3.039994176457527,3.5804622118697416,4.0925903758629145,4.590353099539331,5.083350191283921,5.5794977205477405,6.0859998793016805,6.61025141146697,7.160584311064936,7.747283507953984,8.383989296117146,9.090271349652221,9.896132686400145,10.853499340539878,12.065117982938409,13.789342933771128,17.183040225066282
2023-09-30,0.8175381951646838,1.7514968140527125,2.4222641199449733,2.9934872451196863,3.5137198189086654,4.005943295696447,4.483804732687924,4.9566382482766755,5.432103527307439,5.917141366247464,6.418849356311661,6.945201204470851,7.506018742428262,8.114309318748656,8.788717061826281,9.557807703407384,10.471000847550933,11.626053080297485,13.268687981453315,16.49890266097542
2023-10-01,0.6758642366168103,1.5281732333399143,2.149846447521487,2.6827006328982943,3.1699455011962403,3.6322958581040643,4.082174547775706,4.528151923679227,4.977328661475443,5.436193841481194,5.911432196069014,6.410596743255801,6.94302952873512,7.521137175160717,8.162735097505303,8.8951530998696,9.76570772545774,10.868052271392754,12.437721255548512,15.529894842961744
2023-10-02,0.5736461231423647,1.3775074231695792,1.9738420807871786,2.488521023215921,2.9611499463464983,3.4109965470656136,3.849744902872587,4.285532997816089,4.725173715317313,5.17494964201856,5.641382926164623,6.1318866450388585,6.655668490436279,7.224991482297908,7.857497205806315,8.58028268849988,9.440296386013904,10.530519659314631,12.084929766703405,15.152465815098905
2023-10-03,0.4802547480742554,1.2305415470843062,1.7971241590099316,2.2896068051704646,2.7438084189187144,3.177442856859658,3.6013825246446225,4.023278045526499,4.4496040455362,4.886386631380106,5.339930493003715,5.817443857996806,6.3279171912186385,6.883357503107837,7.501068817614201,8.207661932796292,9.049276696860058,10.117346417926251,11.642074771110128,14.656203178390188
2023-10-04,0.4255560171405164,1.1413329601519024,1.6885376936235958,2.1664217060438893,2.60847947584695,3.031375713488312,3.445465686270562,3.8580589350762846,4.275449086940716,4.703481569520264,5.148314989349604,5.617017406033472,6.118434447965855,6.66450188852462,7.271975271674106,7.967388515158372,8.79617029081174,9.84896649834713,11.353302083752999,14.329256615275163
2023-10-05,0.3686361481483196,1.0506865731290054,1.5805077221040753,2.0460809044210038,2.4782909482976754,2.8928355800593906,3.299558447417576,3.7054556553701645,4.116628024018337,4.538779190540858,4.977960943515847,5.441152830073025,5.937114679509892,6.477679803619778,7.079552616498884,7.769120310701608,8.59156165509469,9.637258110523799,11.132518452128267,14.09568020625373
2023-10-06,0.1829621757808343,0.7169288517230669,1.1652727234233808,1.5704252901252609,1.9526969496415703,2.3234772607174983,2.6903561143593286,3.0590259456070106,3.434617397805488,3.822151759675052,4.2270968956030845,4.655888646850785,5.116715379320089,5.6206566152341075,6.183812406780934,6.831080474505732,7.605774157889614,8.593955633454827,10.012834984503863,12.83991081182086
2023-10-07,0.1488798252999109,0.6468084949266091,1.0747106481273767,1.464413869426849,1.833729882695847,2.1930147150352974,2.549313879838747,2.9079894757798974,3.273939624398624,3.652006148914766,4.047501345826545,4.466710866908623,4.9176599603661915,5.41123187920314,5.963263758549886,6.5982697841598785,7.358919774625418,8.330035257104157,9.72578152732363,12.51046155501188
2023-10-08,0.12197180217158345,0.5878602684622436,0.9972988055644855,1.372915283197073,1.7303339071732375,2.0789995540725075,2.4254743072034697,2.7748250149727114,3.1317399047361767,3.5008965145524926,3.887462498824844,4.297581572568089,4.739122553111325,5.222778019910515,5.764126978213861,6.387303131168234,7.134338270466061,8.088817812935586,9.461853102857955,12.20446217563911
2023-10-09,0.42262574843659195,1.1523655842855949,1.7125871223844955,2.202610458559941,2.656156064098848,3.090429360016428,3.5156234152856713,3.9394240016502033,4.368322238392603,4.808226346900435,5.265477820825203,5.747322504363361,6.262846631710388,6.824214548678441,7.449012013333949,8.164208178852203,9.016618734731821,10.0991303798744,11.6458515219938,14.706304366693807
2023-10-10,35.61931489773349,42.07035212449426,45.88638053601694,48.841690085082014,51.36573557455371,53.639311115206205,55.7593607512528,57.78587262449331,59.76228386949001,61.72326377138145,63.700065658323375,65.72410709059731,67.83079969933861,70.0640727614956,72.48412993049756,75.18031412581702,78.30417293003293,82.15063732146368,87.45020273317267,97.40615404502351
2023-10-11,19.468255653588894,24.299909241406457,27.215534741366618,29.498673718848767,31.464370317038792,33.24648266465108,34.91743079077736,36.5224599560739,38.09476988229679,39.66126099506868,41.246591852628356,42.87595832098041,44.578164963894274,46.38934483982698,48.359414739482304,50.562888232434624,53.12662738158901,56.29828991542093,60.693072081011515,69.0192923030043
2023-10-12,6.269984773440503,9.128365558075211,10.947760422866535,12.412973611104244,13.699530674143391,14.884071415308227,16.009129231945035,17.10198980312654,18.18338926167894,19.270776635190796,20.380794907731165,21.531085728095427,22.742439994977925,24.04154355161633,25.46585694338322,27.07194668019062,28.956813337764608,31.310948030207147,34.61016053195082,40.964670522542875
2023-10-13,2.254450287835515,4.07792311260215,5.322852632305163,6.359656013129408,7.290638038373466,8.16241655022003,9.00184715292623,9.826798839473074,10.651477525813556,11.488382765612055,12.349963630723096,13.249913172503135,14.20483830972706,15.236496764638243,16.375854916775356,17.67012176280754,19.200752013524735,21.128477892487115,23.856539559021414,29.183985207727343
2023-10-14,3.0147963894767287,5.077524819810424,6.454591674334768,7.589622339484081,8.601930495883485,9.545091230699938,10.44958508238284,11.335455172550356,12.218409861331283,13.112081937396189,14.029875307798067,14.986376407696046,15.99913302456759,17.091004725356157,18.294402482929172,19.658607341522167,21.26851593319832,23.291430369842296,26.146564432258025,31.701306901344697
2023-10-15,2.9104211405641465,4.9393954641228905,6.2972226174011325,7.417677268782364,8.417732893050998,9.350001063665236,10.244456172266215,11.120830524697105,11.994612230066362,12.879264749768337,13.788044616429415,14.735395165107148,15.738707267874997,16.82065101161541,18.013384395193977,19.365825132252773,20.962236519025343,22.9687093658928,25.80152107395988,31.315207445903443
2023-10-16,1.7849799556311488,3.3814578083401026,4.48696994028258,5.413567278942445,6.249021162485882,7.033726646028929,7.791146449027226,8.537010399047121,9.28393567065179,10.04311866745902,10.825797309087791,11.644410745898796,12.514115962484306,13.454836371782221,14.49498985658175,15.677968687002695,17.078698017572528,18.84514082797621,21.348757959453984,26.248351118191085
2023-10-17,1.3903674699687476,2.8262359170717457,3.841787331715578,4.700918665937827,5.480130233327845,6.215181930197629,6.927106630904766,7.630165840036392,8.335955000651923,9.054889288698252,9.797538216604206,10.575704486268366,11.40386407907963,12.301266919264473,13.294846674197206,14.426837937488104,15.769436540425167,17.46561674200021,19.874604314270595,24.60256975937972
2023-10-18,1.1839744212756167,2.524710154327395,3.4863866863904422,4.304844128863584,5.049975541873573,5.7548074401743445,6.438934288069463,7.115748873895431,7.796232856117727,8.490329071201916,9.208198251485998,9.961249500047746,10.763531425797012,11.633644228993397,12.598243510187746,13.698161679901409,15.004046435109217,16.65563371309994,19.004230030118165,23.62166344406455
2023-10-19,0.8199962400122132,1.9789637417410573,2.839991418677009,3.583506395875156,4.266510622778499,4.9167514474439855,5.551070241923478,6.18120714200644,6.816999079684354,7.467523698993359,8.14221063476344,8.851782742785778,9.609563351265933,10.433300756695463,11.348531037294116,12.394479710166026,13.639117196011297,15.217063945147618,17.467186488244817,21.907983319259888
2023-10-20,0.8320026734818683,1.9990641359484804,2.865040493359555,3.6124579475387195,4.2988392120399705,4.952154004956935,5.589364696235127,6.22228771796258,6.8608163394277035,7.514074304818392,8.191533944066336,8.903962188188782,9.664733158076995,10.491659278631616,11.410365589881389,12.460211127765486,13.709393619402816,15.292978018794091,17.550938227898996,22.006653316077088
2023-10-21,1.3213224600219549,2.297605413788693,2.956214516186408,3.5015272963964903,3.9893874635977222,4.444977522434944,4.882705776682827,5.312095805641869,5.74065941608923,6.1749576028599815,6.621479538187183,7.087321849650314,7.581053714099549,8.113868057570034,8.701663572669975,9.368647672384315,10.156541607758815,11.147620154220045,12.548486990174945,15.277793852705466
2023-10-22,1.3606961168750171,2.285539793472162,2.902651653928789,3.4109938741023065,3.8642566205918425,4.286477395844271,4.691324531490815,5.087782077383905,5.482889306460354,5.882750959051172,6.293366268339506,6.721260879922948,7.174283643777686,7.662733608615843,8.200870584250739,8.81094720110965,9.530827193328127,10.435353454787716,11.712103447320438,14.194869812268994
2023-10-23,1.411697256541006,2.355935170313856,2.9843853507243296,3.5016928271979304,3.9626708121311656,4.391875437798038,4.803269527264924,5.205975593956691,5.607207721793531,6.013164555829712,6.429941512401856,6.864159514188693,7.323783869888669,7.8192808576402015,8.365031374838885,8.98358962097493,9.713400407150859,10.630143597491235,11.923533235487865,14.438699049814486
2023-10-24,2.852319007435368,4.136613303961507,4.953300938204412,5.61044223561167,6.187200145582611,6.718042149882603,7.222084286861536,7.71158218126856,8.195841937664488,8.682687431098811,9.17957452597105,9.694400289666854,10.236465747304557,10.817703543940759,11.454858749513997,12.173210427022473,13.016103988361214,14.068651011158392,15.543736217333711,18.38300156868037
2023-10-25,2.572277583578832,3.7979652721004746,4.582091652598515,5.214947940159319,5.771543530039918,6.284646394061021,6.77248539824758,7.2467825871272415,7.716473708517397,8.189102068051097,8.671885663447513,9.172497434602295,9.700000605692685,10.266047572125242,10.887014215403124,11.587648789568538,12.41041020131823,13.438714256702141,14.881312035755373,17.66214194624914
2023-10-26,1.610826735371835,2.6062234965228286,3.2624747697960297,3.7999878111149172,4.277449193074135,4.720945687909949,5.145217263179053,5.559890205237559,5.972446894621793,6.389329131325186,6.816821556817731,7.261718335609284,7.732151484954132,8.238676566044445,8.79622416631242,9.427460295678088,10.171383348425975,11.104785737176751,12.420241004844785,14.972512321148633
2023-10-27,1.0203897895902951,1.8364296967672034,2.392891624741474,2.8558980727056977,3.2714668488117082,3.660478756109655,4.0349646845721265,4.402881662669429,4.770617673815294,5.143740681311108,5.527803593462894,5.928909825324162,6.354459490434196,6.814219973706927,7.321749918329221,7.898278626894503,8.580084184718476,9.438491639762901,10.653422111142056,13.024641961176602
2023-10-28,0.9659668594147975,1.762084100732123,2.3069509166770974,2.7612833527971223,3.169583750580951,3.5521361929337845,3.920664896509483,4.282949292224084,4.645244685575242,5.013020100557687,5.391740560760918,5.787423487505518,6.207378505744378,6.661272250499082,7.162484576232532,7.732031770821429,8.405694698143037,9.25454547434893,10.456264329322831,12.803274827609364
2023-10-29,3.01378875628908,4.333163797169246,5.169716574181964,5.841846343071081,6.431160666420517,6.973135515653065,7.487416488338682,7.986578974629123,8.480156164649719,8.976145924569423,9.482154604618797,10.006224040538301,10.557812810895046,11.149132635428579,11.796911222817164,12.527058210428383,13.38345339184682,14.452401377566002,15.94969122465619,18.829573850842475
2023-10-30,4.223889062814202,5.773528411809085,6.739257300801532,7.508316093183707,8.178479238551159,8.791866895997865,9.371609852258501,9.932387527716774,10.485200314922654,11.03916890936001,11.602861619169204,12.18523728290376,12.796737364350628,13.45065385217946,14.165539930800863,14.969298499440354,15.909657122746712,17.080150493514523,18.714319771217983,21.84264540829748
2023-10-31,2.3569397114292063,3.5411697397343587,4.303189287313782,4.919988020346216,5.463532490766365,5.965363225850026,6.443077867283195,6.908026192650701,7.368892972983997,7.8330380382384766,8.307531205691001,8.799913641063815,9.319117194808866,9.876648139823024,10.48869876768322,11.17976269242725,11.99188837503487,13.007724727313565,14.434191556496451,17.18767867722979

Observations;

,Obs
2023-09-21,116.14925
2023-09-22,77.34921250000001
2023-09-23,33.762283333333336
2023-09-24,16.902504166666667
2023-09-25,10.817033333333333
2023-09-26,7.760608333333334
2023-09-27,6.157825
2023-09-28,5.104141666666666
2023-09-29,4.415979166666666
2023-09-30,3.90935
2023-10-01,3.6289625
2023-10-02,3.3490458333333333
2023-10-03,3.048633333333333
2023-10-04,2.9210499999999997
2023-10-05,2.7403333333333335
2023-10-06,2.5623125
2023-10-07,2.3800208333333335
2023-10-08,2.2624166666666667
2023-10-09,31.1369375
2023-10-10,58.83715
2023-10-11,34.48285833333333
2023-10-12,16.837333333333333
2023-10-13,12.270874999999998
2023-10-14,13.359291666666666
2023-10-15,11.410545833333332
2023-10-16,9.407925
2023-10-17,8.507729166666666
2023-10-18,7.5906875000000005
2023-10-19,6.358270833333333
2023-10-20,5.342191666666667
2023-10-21,4.737183333333333
2023-10-22,4.4277375
2023-10-23,6.032966666666667
2023-10-24,9.236491666666668
2023-10-25,7.1585624999999995
2023-10-26,5.490570833333333
2023-10-27,4.976475
2023-10-28,6.67305
2023-10-29,12.142620833333334
2023-10-30,11.530054166666666
2023-10-31,8.564754166666667

[vinfort]

I gave up on trying to fix the code, and instead translated the code found in the CRAN/verification package, with (a lot of) help from chatGPT3.5. Here is the new code:

def crps_hersbach_decomposition(eps, obs):
    """
    This function decomposes the CRPS into reliability and "potential"
    components according to Hersbach (2000). The potential CRPS
    represents the best possible CRPS value that could be achieved, if forecasts
    were perfectly reliable.

    Parameters
    ----------
    eps : Ensemble forecasts or ensemble simulations. It must be a T x M matrix,
        with T the time steps and M the members.
    obs : A vector of corresponding observations to match the forecasts. It must
        be a T x 1 vector

    Returns
    -------
    crps_tot : The total CRPS (reliability + potential)
    reliability_component : The reliability component of the CRPS according to Hersbach (2000)
    potential_component : The potential component of the CRPS according to Hersbach (2000)

    Hersbach, H., 2000. Decomposition of the continuous ranked probability score
    for ensemble prediction systems. Weather Forecast. 15, 550?570.

    Author:
        Ronald Frenette, Severe Weather Lab, Quebec region, Jun 2009
        Vincent Fortin, ECCC and chatGPT3.5 for the python version, Nov 2023
    """
    n_member = eps.shape[1]
    n_obs = len(obs)
    alpha = np.zeros((n_obs, n_member + 1))
    beta = np.zeros((n_obs, n_member + 1))
    heaviside_0 = np.zeros(n_obs)
    heaviside_n = np.zeros(n_obs)

    prev = np.sort(eps, axis=1)

    # Calculate alpha and beta of observation
    # heaviside for the two outliers

    # 1) Beta and alpha for Outliers
    index = np.where(obs < prev[:, 0])
    beta[index, 0] = prev[index, 0] - obs[index]
    index = np.where(obs > prev[:, n_member - 1])
    alpha[index, n_member] = obs[index] - prev[index, n_member - 1]

    # 2) Heavisides for Outliers
    index = np.where(obs <= prev[:, 0])
    heaviside_0[index] = 1
    index = np.where(obs <= prev[:, n_member - 1])
    heaviside_n[index] = 1

    # 3) Non-outlier
    for i in range(n_member - 1):
        index = np.where(obs > prev[:, i + 1])
        alpha[index, i + 1] = prev[index, i + 1] - prev[index, i]
        index = np.where(obs < prev[:, i])
        beta[index, i + 1] = prev[index, i + 1] - prev[index, i]
        index = np.where((prev[:, i + 1] > obs) & (obs > prev[:, i]))
        alpha[index, i + 1] = obs[index] - prev[index, i]
        beta[index, i + 1] = prev[index, i + 1] - obs[index]

    # Compute the components of CRPS from these quantities
    crps_tot, reliability_component, potential_component = crps_from_alpha_beta(alpha, beta, heaviside_0, heaviside_n)

    return crps_tot, reliability_component, potential_component

def crps_from_alpha_beta(alpha, beta, heaviside_0, heaviside_n):
    """
    Calculate CRPS from alpha, beta, heavisides.

    Parameters
    ----------
        alpha (numpy.ndarray): Alpha array from crps_decomposition.
        beta (numpy.ndarray): Beta array from crps_decomposition.
        heaviside_0 (numpy.ndarray): Heaviside array for first outliers from crps_decomposition.
        heaviside_n (numpy.ndarray): Heaviside array for last outliers from crps_decomposition.

    Returns
    -------
    crps_tot : The total CRPS (reliability + potential)
    reliability_component : The reliability component of the CRPS according to Hersbach (2000)
    potential_component : The potential component of the CRPS according to Hersbach (2000)

    Author:
        Ronald Frenette, Severe Weather Lab, Quebec region, Jun 2009
        Vincent Fortin, ECCC and chatGPT3.5 for the python version, Nov 2023
    """
    n_member = alpha.shape[1] - 1
    reli = 0
    crps_pot = 0

    for i in range(n_member + 1):
        index = i

        mean_oi = 0
        mean_gi = 0

        # Outlier
        if i == 0:
            mean_beta = np.mean(beta[:, index])
            mean_oi = np.mean(heaviside_0)
            if mean_oi != 0:
                mean_gi = mean_beta / mean_oi

        if i == n_member:
            mean_oi = np.mean(heaviside_n)
            mean_alpha = np.mean(alpha[:, index])
            if mean_oi != 1:
                mean_gi = mean_alpha / (1 - mean_oi)

        # Non-outliers
        if 0 < i < n_member:
            mean_beta = np.mean(beta[:, index])
            mean_alpha = np.mean(alpha[:, index])
            mean_oi = mean_beta / (mean_alpha + mean_beta)
            mean_gi = mean_alpha + mean_beta

        pi = i / n_member

        reli += mean_gi * (mean_oi - pi) * (mean_oi - pi)
        crps_pot += mean_gi * mean_oi * (1.0 - mean_oi)

    crps = reli + crps_pot

    return crps, reli, crps_pot