JuliaExtremes / IDFCurves.jl

MIT License
0 stars 0 forks source link

Conversion of the hessian into a PDMat object may cause an error #163

Open AugustePaoli99 opened 5 days ago

AugustePaoli99 commented 5 days ago

L'appel à la fonction hessian(pd::DependentScalingModel, data::IDFdata) peut causer une erreur dans certains cas. Cette erreur a lieu lors de l'exécution de PDMat(Symmetric(H))

L'erreur est la suivante : PosDefException: matrix is not positive definite; Cholesky factorization failed. C'est dû au fait que la matrice H (calculée comme la hessienne de la log-vraisemblance) n'est pas inversible (erreur numérique probablement). On peut recréer la même erreur en exécutant :

H = [0 0; 0 0]
PDMat(Symmetric(H))
AugustePaoli99 commented 5 days ago

Dans les faits, cette erreur finit par arriver lorsqu'on applique le test à un grand nombre de données de simulation. Proposition de "quick fix" : dans la fonction scalingtest(pd_type::Type{<:MarginalScalingModel}, data::IDFdata, tag_out::String, q::Integer), je try / catch au moment de l'appel à hessian(fitted_model, train_data). Si une erreur est renvoyée, j'en déduis qu'on a une matrice d'information de Fisher singulière, donc pas assez d'information pour rejeter : je renvoie alors une valeur-p égale à 1.

AugustePaoli99 commented 4 days ago

Mon quick fix n'a pas vraiment réglé le problème. Certes il n'y a maintenant plus de bug, mais il y a énormément d'instances dans l'étude de simulation pour lesquelles la situation ci-dessus se produit, et la valeur-p renvoyée est alors 1.

Précisément, lorsque les données sont générées par le modèle hybride et que l'exposant d'échelle $\alpha$ est plus grand pour les petits durées (cassure convexe), et que l'on estime un modèle d'échelle général sur ces données, la valeur estimée de $\delta$ est égale précisément à $0$. Lorsqu'on calcule ensuite la hessienne, la dernière ligne et la dernière colonne sont alors nulles, d'où le problème lors de la conversion en PDMat.

Dans IDF.jl, j'avais réglé le problème de manière peu élégante en forçant la valeur retournée de $\delta$ à être supérieure à 1/100 de minute. Il y a forcément moyen de faire mieux mais je ne sais pas encore comment. Peut être utiliser le fait que $\delta$ soit égal à $0$ pour tranformer le modèle GeneralScaling en modèle SimpleScaling ?

@jojal5, vois-tu une façon de prendre en compte ce problème ? Tant qu'il n'est pas réglé, impossible de reproduire les résultats de l'article avec IDFCurves.jl...

jojal5 commented 4 days ago

Lorsque delta = 0, il faudrait effectivement calculer la hessienne avec le simple scaling. Peut-être même renvoyer le un modèle simple scaling lors de l'estimation si delta hat = 0. Un avertissement pourrait être renvoyer avec @warn("Delta hat is 0, so a fitted simple scaling model is returned.")