MBCn - Why an univariate adjustment before applying NpdfTransform?

[sylvainmarchi]

Setup Information

Xclim version: 0.52.0
Python version: 3.10.14
Operating System: Linux

Context

I would like to apply the MBCn algorithm to Euro-CORDEX data. I followed this helpful tutorial:
https://xclim.readthedocs.io/en/stable/notebooks/sdba.html#Fourth-example-:-Multivariate-bias-adjustment-with-multiple-steps-(Cannon,-2018)

I successfully reproduced the example. However, I’m curious about why a univariate adjustment is performed before applying the N-dimensional probability density function transform. I understand that the univariate adjustments are used afterward, but in Cannon (2018), isn’t the NpdfTransform applied to the standardized raw data?

Perhaps it doesn’t make much difference. I tested using the raw standardized stacked data, and it produces similar (by eye, I haven't investigated much) results. I thought the reason might be related to the convergence speed of the algorithm. However, the E-score converges similarly in both cases.

Thanks for sharing your thoughts on this and for taking the time to develop this useful package.

Sylvain

Steps To Reproduce

$ pip install foo --bar

[coxipi]

Hi Sylvain,

You are right, in (Cannon, 2018), we first have:

1. Npdf transform on raw standarized data
2. Univariate QDM on raw data
3. Reorder QDM-adjusted (output of 2.) using the Npdf transformed data (output of 1.)

In the xclim notebook you point out, we had:

1. Univariate QDM on raw data
2. Npdf transform on QDM-adjusted and standarized data
3. Reorder QDM-adjusted (output of 1.) using the Npdf transformed data (output of 2.)

I also concluded there was not much difference, but in newer version of xclim, we follow the exact methology proposed in (Cannon, 2018).

Both the standardization and univariate QDM are now performed in a single class MBCn in recent versions of xclim, you can see recent notebooks for more information (I found out the first two cells in the MBCn example in the latest sdba.ipynb don't serve any purpose, this might be a bit confusing, I will fix this for future releases)

Let me know if you have more questions.

Cheers,
Éric

[bweeding]

Hi Éric, I also have a question about the change. The new example code reads:

ADJ = sdba.MBCn.train(
    ref,
    hist,
    base_kws={"nquantiles": 20, "group": "time"},
    adj_kws={"interp": "nearest", "extrapolation": "constant"},
    n_iter=20,  # perform 20 iteration
    n_escore=1000,  # only send 1000 points to the escore metric
)

scenh, scens = (
    ADJ.adjust(
        sim=ds,
        ref=ref,
        hist=hist,
        base=[sdba.QuantileDeltaMapping](https://xclim.readthedocs.io/en/stable/apidoc/xclim.sdba.html#xclim.sdba.adjustment.QuantileDeltaMapping),
        base_kws_vars={
            "pr": {
                "kind": "*",
                "jitter_under_thresh_value": "0.01 mm d-1",
                "adapt_freq_thresh": "0.1 mm d-1",
            },
            "tasmax": {"kind": "+"},
        },
        adj_kws={"interp": "linear", "extrapolation": "constant"},
    )
    for ds in (hist, sim)
)

I'm confused about the nquantiles and group parameter here... In the original xclim example we could perform univariate adjustments using different groups on each variable, and then apply the N-pdft part of the MBCn process with another grouping. Does the new way the example is written mean that both the N-pdft part and the univariate parts are performed with the same grouping?

Cheers
Ben

[coxipi]

Hi Ben,

Yes, you are correct. The previous implementation gave more flexibility, since you would simply compute a univariate adjustment(let's say QDM) and a multivariate NPDf adjustment separately, and combine them outside of a class.

I didn't anticipate that mixing groups could be something of interest. Is there a motivation for having a different grouping, is there some work or insight you could share that highlights how this implementation could be preferred? Or is it more of a performance issue?

The motivation to perform both procedures at once was because of the rolling window. With Grouper("time.dayofyear", window=31) for instance, you train 31 days * N years, and just keep the central day in the window. Then rinse and repeat for each day. This is true for NPDf rotation and univariate QDM. But I feel that if you reorder a sample of 31 days.N years with the NPDf, you should not drop things yet. You should then perform your univariate adjustment, also on the 31 days* N years. Then, reorder this adjusted sample of 31 days * N years. Then, finally, you can keep the central day. In the old way this was configured, this was not possible.

Best,
Éric

[bweeding]

Ah, rolling windows make everything fun don't they!

I used the more flexible version to perform a "cascading" bias correction. I'll be interested to experiment with the new example code.

The paper and code showing my implementation are:
paper
code

I'd be very interested to hear what you think if you get a moment to have a look.

Cheers
Ben