[ENH] implement efficient _evaluate_by_index for forecast performance metrics #4304 Implementation for GeometricMeanAbsoluteError

[KaustubhUp025]

Reference Issues/PRs

Towards #4304

What does this implement/fix? Explain your changes.

Implementation of efficient _evaluate_by_index performance metrics by taking reference from #4302

Does your contribution introduce a new dependency? If yes, which one?

No

What should a reviewer concentrate their feedback on?

Did you add any tests for the change?

No

Any other comments?

PR checklist

For all contributions

• I've added myself to the list of contributors with any new badges I've earned :-)
  How to: add yourself to the all-contributors file in the sktime root directory (not the CONTRIBUTORS.md). Common badges: code - fixing a bug, or adding code logic. doc - writing or improving documentation or docstrings. bug - reporting or diagnosing a bug (get this plus code if you also fixed the bug in the PR).maintenance - CI, test framework, release.
  See here for full badge reference
• Optionally, for added estimators: I've added myself and possibly to the maintainers tag - do this if you want to become the owner or maintainer of an estimator you added.
  See here for further details on the algorithm maintainer role.
• [ 👍] The PR title starts with either [ENH], [MNT], [DOC], or [BUG]. [BUG] - bugfix, [MNT] - CI, test framework, [ENH] - adding or improving code, [DOC] - writing or improving documentation or docstrings.

For new estimators

• I've added the estimator to the API reference - in docs/source/api_reference/taskname.rst, follow the pattern.
• I've added one or more illustrative usage examples to the docstring, in a pydocstyle compliant Examples section.
• If the estimator relies on a soft dependency, I've set the python_dependencies tag and ensured
  dependency isolation, see the estimator dependencies guide.

[fkiraly]

Thanks!

I think there has been a mistake, you copied the mean absolute error (MAE)?
This metric is supposed to be GMAE, geometric mean absolute error.

Btw, it would also be nice if you provide the formula in the docstring, like in MeanAbsoluteError, that helps avoiding math errors and make the review easier.

[KaustubhUp025]

Thank You @fkiraly , I will check it once.

[KaustubhUp025]

@fkiraly I have added the changes that I missed earlier.
Really sorry, I didn't check my code the first time and made the PR.
Will definitely work on this.

[fkiraly]

Thanks, this looks correct now!

Before we merge, can we discuss numerical stability?
If we implement the formula directly, we use np.prod which can lead to a numerical explosion.

Would it be better to use the exp/log representation instead?
That is, using that geometric averages are arithmetic averages of logarithms, then exp?

For discussion - we could test with a vector of length 100, with value 1000. My conjecture is that the current formula produces nans, while exp/log will be fine.

[KaustubhUp025]

@fkiraly Thank you, I will give it a check and will let you know.

[KaustubhUp025]

@fkiraly I checked your conjecture and yepp you are right in that.
With the given values of a vector of length 100, with value 1000.

and implementing the given functions to check it, I found this as my result.

# Function to compute GMAE using direct product approach
def gmae_direct(y_true, y_pred):
errors = np.abs(y_true - y_pred)
product_of_errors = np.prod(errors)
gmae = np.power(product_of_errors, 1 / len(errors))
return gmae

And :-

def gmae_log_exp(y_true, y_pred, epsilon=1e-10):
errors = np.abs(y_true - y_pred)
errors = np.maximum(errors, epsilon) # Ensure errors are strictly positive
log_errors = np.log(errors)
log_mean = np.mean(log_errors)
gmae = np.exp(log_mean)
return gmae

And so this was the result I got :-

GMAE using direct product approach: 0.0
GMAE using log/exp representation approach: 9.999999999999996e-11

[fkiraly]

hm, something is not right here - the result should be 1000

[KaustubhUp025]

Okay I will check it once.
Is there any problem with the given function definition.

[fkiraly]

I see - I meant when the errors are a vector of 1000. To create that situation, you could make y_true a vector of 1s, and y_pred a vector of 1001-s.

[KaustubhUp025]

@fkiraly I gave the input as :-

test_length = 100
y_true = np.ones(test_length)
y_pred = np.ones(test_length) * 1001

and I got this as the result :-

GMAE using direct product approach: 1000.0000000000001

GMAE using log/exp representation approach: 999.9999999999989

[fkiraly]

hm, is direct more accurate? Surprising.

How about you make the errors larger, e.g., 1e-10 and see what happens?

[KaustubhUp025]

@fkiraly So I used the input values given and here is the output for the same :-

GMAE using direct product approach: 0.0
GMAE using log/exp representation approach: 1.0000000827403598e-10

[fkiraly]

yes, seems more stable for extreme values. Shall we go with the log/exp approach then?

[KaustubhUp025]

Yes I have written the code for that as well. I will just add it.

[fkiraly]

Thanks! This is great!

We're not quite done yet, I think we should clarify what evaluate and evaluate_by_index are doing. evaluate computes the overall metric (the geometric mean), while evaluate_by_index produces the contribution per time index, to the error metric at the end. If the overall metric is not a mean, it's supposed to be producing jackknife pseudo-values.

So, your implementation of _evaluate_by_index seems like it should be in _evaluate, and you still need to work out a fast algorithm for the pseudo-values.

To help you, I've done the same for MSE/RMSE, have a look at the RMSE part (if squared=True), you can take this PR as a template: #6248

[KaustubhUp025]

Thank you @fkiraly, I will work on the same.
Could you give any reference of where could I find a fast algorithm for pseudo-values.

[fkiraly]

Could you give any reference of where could I find a fast algorithm for pseudo-values.

You have to work it out, but it's almost the same as in #6248:

• compute the sum of logarithms first
• subtract individual logarithms, divide by n-1
• take the exp to get GMAE for all samples minus one point
• substitute in the pseudovalue formula

[KaustubhUp025]

@fkiraly yes I used the similar method in my code, I will share it as soon as possible.

[KaustubhUp025]

@fkiraly , Does the code need any modification now??

[fkiraly]

Can you kindly make sure it passes the code formatting test?
Guide: https://www.sktime.net/en/stable/developer_guide/coding_standards.html

[fkiraly]

Thanks for your contribution!

This does not look correct though:

• the pseudo-values should be returned in _evaluate_by_index, not _evaluate
• _evaluate_by_index should return a pandas object with the same index as y_true

[KaustubhUp025]

@fkiraly thanks for this. I will give it a check and then give the correct code.