From 11182b47f8347851028af67935ee32a3366d834a Mon Sep 17 00:00:00 2001 From: ramseshk <45832522+ramseshk@users.noreply.github.com> Date: Tue, 11 Aug 2026 10:50:05 +0800 Subject: [PATCH] Fix ML calibration: logistic regression + Platt/isotonic + realistic NWP errors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Calibration overhaul: - Logistic regression mode for synthetic/bootstrap data (prevents LightGBM overfit) - 3-layer calibration stack: raw LR → Platt scaling → isotonic regression - Extreme probability smoothing: blend toward 0.5 when raw>0.95 or raw<0.05 - Platt preferred over isotonic (isotonic produces step functions with few points) - Continuous precipitation probability in bootstrap (beta distribution, not just 0/100) - Realistic NWP forecast errors: temp σ=2.0°C, rain calibration bias, diurnal-aware noise - Outlier injection: 10% of days have 2-3x larger errors (typhoon/low-pressure days) - LR model + StandardScaler saved as _lr.pkl alongside .lgb marker Results: - temp_gt_30c: AUC=0.987, Brier=0.049, predictions vary 20-85% per day - rain_gt_0mm: AUC=0.979, Brier=0.042, predictions vary 15-85% per day - temp_gt_35c: AUC=0.713 (realistic — extreme heat is hard to predict) --- data/models/rain_gt_0mm_24h_cal.pkl | Bin 0 -> 1577 bytes data/models/rain_gt_0mm_24h_lr.pkl | Bin 0 -> 2203 bytes data/models/rain_gt_10mm_24h_cal.pkl | Bin 0 -> 1513 bytes data/models/rain_gt_10mm_24h_lr.pkl | Bin 0 -> 2203 bytes data/models/rain_gt_5mm_24h_cal.pkl | Bin 0 -> 1545 bytes data/models/rain_gt_5mm_24h_lr.pkl | Bin 0 -> 2203 bytes data/models/temp_gt_30c_24h_cal.pkl | Bin 0 -> 1641 bytes data/models/temp_gt_30c_24h_lr.pkl | Bin 0 -> 2203 bytes data/models/temp_gt_33c_24h_cal.pkl | Bin 0 -> 1641 bytes data/models/temp_gt_33c_24h_lr.pkl | Bin 0 -> 2203 bytes data/models/temp_gt_35c_24h_cal.pkl | Bin 0 -> 1449 bytes data/models/temp_gt_35c_24h_lr.pkl | Bin 0 -> 2203 bytes data/models/wind_gt_30kmh_24h_cal.pkl | Bin 0 -> 1807 bytes data/models/wind_gt_30kmh_24h_lr.pkl | Bin 0 -> 2203 bytes ml/model.py | 418 ++++++++++++++++---------- ml/predictor.py | 10 +- ml/train.py | 414 ++++++++++++------------- 17 files changed, 477 insertions(+), 365 deletions(-) create mode 100644 data/models/rain_gt_0mm_24h_cal.pkl create mode 100644 data/models/rain_gt_0mm_24h_lr.pkl create mode 100644 data/models/rain_gt_10mm_24h_cal.pkl create mode 100644 data/models/rain_gt_10mm_24h_lr.pkl create mode 100644 data/models/rain_gt_5mm_24h_cal.pkl create mode 100644 data/models/rain_gt_5mm_24h_lr.pkl create mode 100644 data/models/temp_gt_30c_24h_cal.pkl create mode 100644 data/models/temp_gt_30c_24h_lr.pkl create mode 100644 data/models/temp_gt_33c_24h_cal.pkl create mode 100644 data/models/temp_gt_33c_24h_lr.pkl create mode 100644 data/models/temp_gt_35c_24h_cal.pkl create mode 100644 data/models/temp_gt_35c_24h_lr.pkl create mode 100644 data/models/wind_gt_30kmh_24h_cal.pkl create mode 100644 data/models/wind_gt_30kmh_24h_lr.pkl diff --git a/data/models/rain_gt_0mm_24h_cal.pkl b/data/models/rain_gt_0mm_24h_cal.pkl new file mode 100644 index 0000000000000000000000000000000000000000..f61703eece4a47ad73f940bfacfa99dc831504eb GIT binary patch literal 1577 zcma)5ZD<@t7{1GW^h|Q?^}`@-3{8p9(Dsr>@vGuU4%GD8P+E|J!muB6x2twS6Y9I0+a}^!d|CFH znVt7}KlXXM!UJ{DfM@TDojjT&OrxscG2Nz&nz|x*ByUKJiPq~Cfk0Py;5z`8y{@2tw@kYR-}WxS}`Qs z#vMe@YD~2AYD+Rzn!`4e7=eA%rp7U1L1BiR$=XnBPCAZtpcQGASLGao>Dc2Cyiz53 zizwg-Jy{1{GPQ|F2Y6+NWO-lJtT~mL)+ysxd^Qrp+uzVbCxh@oqwVQ(d`o3Nq*ssh})+#bL2C z9ZXmB2(6Y!RMwd%u1kgftSI!Yd=2s73!=TfegAzh@0AWPAs14KmA(c7!O}v&F+}lt zfrYkBt8?RM+x;j(*gOH;Ppxumjpdlbj{u{`7Evh4l;uhXPk3yZHJO%Yi?xpL#!gTD zOCDq8>QDcV6j31+d;YhfbKmYyjC(E#T#_n=Kj_JjNA@P4Nwx?_P6{3Zt&7&jLY}`pq++{zYl;*1=1bq=M={et&KfmT(6CnhFR^MH9wyTo5YU!F>bxJ( zzWYW!svSfl13!tgj}IcUC)3>h`WR}SyxrD(ehiJYye8Y(tLUrcTQe<-SJ4mizm9bc zO`(0Eru*-VPN5@1qet#+o<-_6FL6ddK#5#V0AWWC11^VE0+C}FM3hLBwdtAZImld2Hz;dR3{HU9 z#26YCe7X_PfKLGr5;SHzus%tm5*Bp{FCv3(CLqYNz6EvFxW%AJ{vdDbkFNT4SABlf z-}-&4)TZVm2actcFu}B~MkgcG484SWM@`@ctx&Dfq8_un7Gi3hLPb(4DMK&Fcsa}_ z6b6DMRXQzG&4lb`W+?W0LW`@ZVurR-M>4dHtPoc-m9#abQ!^=GQS{w0pQJy~+$0oS zsGVtb^gOwW!c7cljul-5pigzZMeb0G*IXV?6g>cPC_EOgE|5;;FL;-#;#gSPT*7_k^%hz z#)t&8gN2J?hF+*G)aZ+am{eyVgqEbZ!GI%aTriYTrzI&vp_D?}(YBVED8E)VJEa=! zz{n(VB}_S!%>;|MA}b@K^rdZORI#3zSz@JH$y8R%ezm5#@l0lBX8H@%EDsTfVkDJG zyxHUBaJVlA40?^SSj4~>CJiR6>iOJpMSflcp1xXZ04qYX84+I zyfgl#)ayKO{-8tz7uB&>N4M}RK4}7b#XoFHt6R%%R8RjoeV~c0j>pplViPO>bkL=C zOT@&j*H1f&Ywoho!K28&=t0&NP8ASoFd+emf#Vx~g52w7vk79(+uUX9E_ z5(Kh~(q}d?Za@{S3{7(2c+aZAypO6L-AITxS7ixMl_pkNs-|Pf@%2mNW1(?uqnFqx z0`@md#Drx}%P%dAXyK`$c>U&x?<8Q*hO{IEWa|e*a3D>cF2ov&#`A@#! z!yWEB-<`r&BF=vKMvQSZ|k)))+< zqRLz@vt6(yC9P4i>to<+bY4?lmCIIiP4(7?RRZJKx{=iNrL0qbnN?i-X0~IAhhj7m zvI$Pk=LMt1;8oS64xfE`QMU(iH@JTo)d0uONHYooWH5D}VAt)flVIqyxO|WoJzD}c zMd&`Cqh5Za=iR(_Rbe#{9$Ht!&3-!(PJA1>P#{l$zTMBG&JTw`NN%nfpU8#6g-(7Z zr2_t7t^Y^kiq$Z~yY}|PNfqSYs@R>I;|3did-Rqdg%{z^RzTjt;pCq7NH}(VU%_@F8;0zJ+jGUc;Uz19RXHw# z8moOHHI6d)tvN)G0|D)*A@`q8uj-IeT$hki; z&I=Zu%lbLG%>!=QW{%12c7fX$PbV#V5(sO8vLvlrLScOG@(gp4C(OYc8g4d)!4dIg zS?K{TjIaqjnm7^wr`DWwj^huD6O!V38qA>npJdCWIx|>vEZz6V>;=%T z-T&%aVm_R_B``V$`$MuSV^v3a4D_?Ue)Egk*I@UD&q@^8A+XdY{kW5rAN*#kH+ig$ z2Q`t`!;dl!@Uv=H?<0LD;MJ+eLr8&YbZgqHdGfiafxO5N=F+u(?J1#Xu(;I-G)Pa zAMd4h$>MTEbbG`#SsDABa0=_R&zwEfgTfrWSn{FPGTDc1kd*~ftSko(w{fH)MVCn{ z9yQ<~^kfB`(wR?XAtAc<5-*bUIAoo6@fg#65YIS6*dyF#Fkj{3A|*V4+97Qm=9Y$ z&@nEpevvYHwL-%q9rhMkpKN~Z6v|dPGBUFLPZ-gwkl<3UYq@4$6A7g?k+2O>qF!KW zJz}j1<5&Z66a^NKgZ9QHZELWda52a*SXrYIuZ z2WhRc`n1gYiTwGwdpBF(vZI(=fV0IM6BxF-{{^e{(1S4T1%UeM+Iw-t58nLd<0rqm zh#Gfux#>;6pf{`^pV<1@JQ|!)H@|rD5?Xj??ETq;SCFM0JvcRc1=W<=(~Z&J(ChgV zckk@Fiawl}`N+F-6&+if*5q3aL|z&Gws!w_bZ+|lm&W#7L-xShY~!at(9em7XKT-2 zM{}PRR-Zg^1Ld~&`QxA8K#SFd;Ag9UFK)nkG{XO}82^7}A6LTP(wK(XzK#@8se y$Bu^MG>UEWw*87?+w_}n`Ol+g4_)4J=yGGxx)-L(7pAVG*v47M_J1%QsDAzkO09-D%mnUbOcZM!p}|#DF+=lIVGPZe72qnSingG%Dkcf6h#a5vO8h0wO+?A% zqoz1UFUeO@m{LO#dMTl!m`d9DhXf;~I68)$N#)WGQWZ{;SP`L=D=3Dxk>Dy^BPC=E z?V!gsGOZdTDV!phblQs4stO4`LtCmO`EnAiRpYxbB`Uy(hpcIPb%BafVrCpe+sO!> zo{*xQh-8hnNwg%vq;hDxB3y4~>M8UBRvK)(Rzf1WjVct=IAExvRzp(y0x5;Gqb6|NS4jL_~egI<_Wdl_#*EJce-E^V=`jQ}4~n zqO;F{oa8?8<-sf9+KnYsf0;W4z6x`G65lYus^6Gv7`S^1wBCxj+BE(xc;=gb$iL@0 zIC;G>E&P@vd@&Fn)YH(x&V;6mBk!MQecs$~;gGx+^o_HczZ_ihcFa3wpk{q_*BVoP}~7dzSgqayeDoyP&^oYkGw&_A>MK6BsLI?mT`Ny)UL6+$4U^gK zVy{p*+&mM*?syw||0kia4fw+ReAk@PotbDHY}IQ)~Ip`;@)?3X?h3MS=u^LX^;3i$HJjxXgBQFI_{rcvi}cG+^c>H9Um08r;m4b5xtYHf!Rc#3u7dn{=-vO{@QY)C5R%(!r>1kD z(A93$s6r0^VWInLOT=23=5h4v=`%{my<53ACul*yOhr2o4l=_y~wflz(J*U61Tgz<}se7G9!BpwG5cQ+V1w9vB^M4dwUGX?Pz}5{8 z``ZtTV%=e3N5*_)y9?Y}WQ@w`cYwP*pC-Qb+z;0JXGqS!7X;&m)}$E=-C!2p+mZuJ#)DmN33h1i@SZ!SD7WS d$!W~3;H(i$TsXK%ln`fDR3TQR7f1?({{vOFmi7Pu literal 0 HcmV?d00001 diff --git a/data/models/rain_gt_5mm_24h_cal.pkl b/data/models/rain_gt_5mm_24h_cal.pkl new file mode 100644 index 0000000000000000000000000000000000000000..1902904fe7e0668bfa523f0012fd646a7ec796bd GIT binary patch literal 1545 zcma)5U1%It6y8mCx0`g+rXh`BLye-)Qe4bdtl-ZL9c z7A20IwoDr|E-;;1X>75)>2g1??Lh6=PxB6SU6a}JxZHV6`k6KqsZA_ilA^)NND zgFu+c28dFbOJz19nvW1CfOI%iov`uSOmzV~X?0_VaFfA#lZ|tfa1VNivaw%VGQb5} ziVoW=6iaEWF^8suq{MLu%p-*%BePxZcp8U+iFn`%`E9-GavUlckxEpRW3pfF$VRgb zBhu&;ah;cn)DLMepO@wF>fc5YeNp!H^}T)%B=}?_T&iWIwK~^GBw9O2xP~mjEU+}T z8Qc9lMpu|63dU1_{p=dHv9lg>tO*bUw7AA-t}a$KIvSwm)>T>`E&1gqi-%_`T|EI} z^V*~TM~SGcBy->WcKZb1dL{^I#Wkv8c)@^%BHo?elg>y_O-s=TT$g)qSLyv{NblB= z-fNrC_F8Qz>Jd6#!UfX?x*k+}=dMlO5&i`^3Kv&z&7SWV+a;iK&7s8gc{(myn8z{8 zVU>Bd4lImmmRXcySeU$>+hMBsh|%t#F-eV zpphHP^KYE3po{PIlFx2b(6P2xhZd%;pvB>i)H5?z&}z5-BD9qL-E5cfM6eD|hU;*D zH(^)%9y!E+dF2}v>W@G2bCuiLJiJ-1PZm&IX&1pJ&!{DRv=3#)K#378r}-} T;&Wdvuj;UiQ2&n|@Y4SP`w@t- literal 0 HcmV?d00001 diff --git a/data/models/rain_gt_5mm_24h_lr.pkl b/data/models/rain_gt_5mm_24h_lr.pkl new file mode 100644 index 0000000000000000000000000000000000000000..a6d2970510ec4b88e2d3bdfbb0cacb53721d6370 GIT binary patch literal 2203 zcmai!dsGuw9>vH*@dj zcjtS5-@DJTs?vpHZ^ev|cGBwQgodFPkZ)=T+^7?2)H>8-_SYgzqgSd)N-bmPIT_D~ z*@V(akfd6#V~#K(CCn6MmVwaW8mfq)d74Ouc9ffN4O31#Pb z!NuyScATEAP*a#%M-fIDVW5~Y+Vf?Akun@@L&>CbX?K|hCrPZ3P%BjwL%T?E4X%?B za)x#{;ySrri;)yg5lklSOzJfSgpr}0G*X3^=Ry&!xRKgnUDa zCkv*Jm&4&cA28@O%3>D-V^j-Ni>C6Zf?j0dqCkk)D}H6G!e$~n`XQp7<|NY-iWyny z!Ef1Vir?&XMw(cN_UGIfYF?Uu9&BG-mej^Q3mU(IwwM!p!Dgnjx-*UeBf2M?1oLH3 zXsc1zWwo>Ur+PNXhO5AUU-_3`Upop+^}MLbv|{l2wplPLC`B9GY3!?)OfHs`Hz zkGsJBDA9bmd{Y&xnbepS-yda#CCwKrC(nb`;gy@B+dg7({^Om)%eYY9a@R}n?kAw? z+pvVADo;4}_h6s+T`Azc`p$0gTrhV|oHM*H~p6L3x*F=L82X_8qeQ<$& zaYU*7c;VfaR<@}%I-Do31&6*-oimjG2bgckTbGqz2Zds5oT29wYZsmWt3`h26vgGq zE)RP9=)q^>(@O-{b7Nzu$nB=3ogzdrqDaOhMV3KEo3hOO7w$)g9e#{ zBnV^|mH%{N+=wb#6`JhAnK!Kl^A@UlbR!|&OqCUks#IKVubS2cr*jr3tb~@4H~b|2 z5wO=fwvv7GH5l+Xsh|-ALWRO7chB~0)OR&*#aF}2>kHs5$A11Wk@m%u#9lp+Y;b}@b(X82?H8Cvi`)bD*?UMrrR|QhhXZ!@pEm0 z~|($K_k)caows~5tk zD6{(&wjDOEeZ5|~_b>>sxUDIz$Ysmg$Ghvo%7JBM{a|X&K9=8W=Eb*cWm^|`D~BQ> zo5**&5Ij@_URKS(`KN9A1Bko9>tK`>KKw|Qksm0Bsk4Lk-q}6|`p!y9|0Rfq_Ns9# zK==6!_3|4%?`_+y4y%IjFYBwg*_*|1^ovlRU_~Mf==eG6!axXwt9OfuHe$d8K$d z)bb7uR=LXI7lxDyZ^tY+5lMCSt0LfVjK%!!Kp@0~c~QKzd8>YWIwaDK^?^qI-|9xb=$T z{a`_J)-Tab-f&ByEhe|a9qyPvk-YRt5Ug5~B|Z1YP?*rYJi}Jt3pe6cYfnQM9F$y< z?>owc5spD8#e;!xe9Z^0o*)LkW*PEz<}89miOKO@RvW1Mj%@twxDBkSO$+!wdk$RK zvgrCI3A(>-1zTK&iy&E%kGB-IxoNmiBXjO#E|2$3S}1y1ZeEf|s6dTX%a_74-IV%l$`oF{~2|_?kYn Wf&1s*UN25cuq&#Ps4`7bljwg9t(_GB literal 0 HcmV?d00001 diff --git a/data/models/temp_gt_30c_24h_cal.pkl b/data/models/temp_gt_30c_24h_cal.pkl new file mode 100644 index 0000000000000000000000000000000000000000..a8689cee8ded626c4e31c62aac2693ed38e81aa1 GIT binary patch literal 1641 zcma)5Yitx%6yA5o(&Dx@jVVoSD$o*S6CME{$Ti&`Jk~aagxG}Wb#`WUXJmHnI&&9V zLlZFevFRldx%>d+hsMx?F`z#bViOgE7MmbxHFX0qQez>k?TS%}N`hymJG5#go=oPu z=gvLfIrq$W_680;7xFmP#hX|trBIc~oC%smC^0mNvtU{mDdo*$K96T_;LtFDZLH4P zsVkysM0Cvn#3>>xdIakvsaezxtbDq-+@>T=#j-SF@PmBWK5l0!Nh^k^(@xHO`YO%> zaz@m7miZ~s`9@T;b|POJUBSFjaDFN#EKg_@YX((JNl8=Q$13hR*pfsjBgtc4Rwn79 zW#Qe5mQ*QkWo4#l$Rvd=DpCdfkVW)eiU}9fiqElaSw#XO0|%3c7jJ8|S@z z8*hqvW4>;#*0P{XJJX6C;?o}G*}j>te&$`zTU%SV+yQf5F%RYOY&1O6SC7YATyWTi zEKV;lR~M;ub{@6f&Ef@T32?8OHKWMH*XIb zRnI^4l*3p$`_uoUM3jvRkG=BSn^#EH$4*EX7p02f4>~er!IhD9kyShsk8tJ(t@Gw~ zza#G!0NYsDmG|rIROax%g^aCMDipt!w_~b;oK>RBpkcM_Ut(47zn4sFKtNU5 z=p{E}>*}jpujWQj+nsdnFPBG9Ul-Yb{@f^9|I#zQHRngsp6r*OKD=ZcZ7;uE+1Njh zI(!>z2KG&$<10EUcf3D=$bm)ge)Q)A3WtYZ9&7v_P5d1?y1 zcI@lT8%|H525)n4XYVxHvigbkWBt>p`}d)~-hl%8aPfDk!NCH$+Z}sw$~*twoO-pv z8K1@O`26$sZWo$Py)o4D^GW1l?dLDTGlX0${cK0>+N)n67mpmaEr(qEN*H-^|LH5p z#WQu?`K@QZMJ{f4><8rHo*MD&K*c0-aRUU833=pV<5B&+3rq9J#mC=T`XEsX$VK)_P3IQom5D+Q?QbJ`T4-qSmd*vaN5L6JQXdNcWkc3ViW(JVs1CXY` zo}u^{)>V&K2wG|@qLzy2*%5T5T5vfoEpl4R!)_J}MUK0$RqL{sptgU|v-gj=_cwFz z=Xd9Of8SeTS98FgXKi~~F(c6Gl?2H$OR4ur0ypX;q(+B&%=%h_k$RPeqBRPZS(y4_ zm`4C}x?(x_qsnNP;Q!MnYmuiW`kMlEw)`8TC4fHs&j6q#Yx$&P4fj%GoK0 z(GE;(JYUY1v+1l@%9q-hShWu$P|`&PVrGerW(`|eG5gh);YYJ+X=y3bs97IU9?i-t z<6UNtm&fD37_jIy%3>7*W7SL5zB75${x7reQ6NU_Rll;;!E=%Ad=b$~6DagV?wl-X zzz0@Z?r(NFCoN8l_7@2sB*vEB1G3a&-?bOJfvEU`U)TNDV1srzv3qA1sGjc0PJH$$ z_ro4@ZQrRz?&w%Sfi7$ke5ojFzqa99@Dcy|8QZjzp!bTP?X9en++VXlB?GftIP#NQ z{C7vrb5C!lO^%7rf#d&jdY&9{9=v;VdWq+f7Ow4X-^z|ht(@wOZjY+d*TL4(13OGt zk8yROm-auoRm(->i2Q@<8iCV0X|btMf8ogeu~9cJmw}R}?uXWAm^npP(Cd}bOJMk$ zgFHXx5qI;A+8u`kU0~CJiKGVp02kuZv2wKc9;hA>vTZlpfy`Fme?9XuN6OXr^sZ;Y z%Oi3aPiXF+Xb03_1fJ*rcd$9{i5KTy&)uKm3gfzroY(U$k=y&qIIHNwUoG-Grx;fm zZyzgZ-dseuUb8N7!Cn{}Lq%>iEh7k(L`uTgd^xfVGTMw~yk5E=nZjVyD`xzyQ%Nu- zPU9@&oIxY2;YMY~?DSW*cPP}*YD`PuI^CPp65BM2Kxc3dUPWp-dvRxfT}dE(pojGOHXAik6sH+qgy>> z-l4G1Jh7I0F9Q1fC$^vo_`%%V#aa7yrlbB&(>6R2UY|_%=yC;+#|zPOA_!x)4~d>! z5y9QA`}R!YYY^vP?t!3q=wlk?b+kppDRKW#Wil2tb#DA4x1kjDnVb7hG*-jpp=0OT zTvl>TNxL%#H46~;pU!Uei73y)#s?LSh;MK%=t+jcKL5tyz*xlj?9CN1A*lCm4{i{{ z@bEIH8(ar$N_x9NUVI4nm>k!aS7mV(9g{uv!Ii)?x@jaiql6RomDxnKY~xy&x~U$9 zK`ut<_<6v?qE}UuJapk%oBn;o-R$~dxEYR}SET0oDPi)0fa33VPJn^4vhw>btD&_7 z>WL#o~b680gdaJpA*aAPA{#HDlvhP_kIK zY*4L&-`X0!X;`%mrg|RzZv2b}^1rSq&B|N^%|3Sw?OV4%&SHCO`GXSZUv>AflS?}M zIrWo8`|Sq7wW+Es?9l-D^3vt$kY83omm@<7cUr>WXm@qqE+QQU?Ss3rWTo(`O@K`~ zE`?f~>X8};CA@7&sB*L02&cm6u0eGu{4vs0c5=uMV&a@|o1~m||9>mKF!5;^;iUxtqE2mYmIt*S^9Pi}zuYniqj=bG?qa2;rY&$aDknmC9 zQ@ESQPpNH)TK9fWf!Fx=!gc(*IBkc^AQ;=15TgC5Vph*t*ZdzvmzO^d7dyDazCg#m zxF`=;aB<^*S2w%C>_SUqR;Lr()Qs53VakII%F&K`> zZYWFs#D}4F{WVKe`t{$_1^>rye|43g5 zm$mqQwM8a^6JH0I9K^nms!EM*EsuoDY`c4}9DN;je)!8?ReBICu}e86v{?pk?ewJT znq8natUIKZwTE9EcJ@4R_cZ)!^5=5z;}($8nC-Eq!2+n9wVQhnS-|&9*4!Dke4PYAU zG+tJSpsOK8RsnHZ6NwVSik6ZMY6iAEQCsdIDV-RGtf_pG_a5YCriKhr1%>80t5bGz z<`T1l!i%h)Y6^c8wGS+nJ0jbeGXl<~T0(zPrdU=f(Gw&?`2@T5hJ%fS0A-|j)XAC> zieMOcACXfMJbU@n8Z$S99Bfpso|_LA`WRH6(1~D} zI-OvCElaVM#BnXFil9PQLXk6^HLhR6aLY4lvffNC9#i^)zELO}bg;3zo9iIrT|0Go zb~rwQ_+e0!l@aEFOU+VS#M@U0&Y!zlsh#;pF_l2hlBS5DVX5g~VwG;WnM~V3KyAgy z4|c{b-@88#`p%*c@7wn2i-5=9`)yat+yde!7FVWM7SP^%4vt>?V-XGQ{pREj$FJy6 z*K6PY)OsETr*<4|-*O(E?Ek2HEKowZu0r#(ttB+AA0e?TC3Jk;Kfk}>0(u~Pp}Tvw zjKb;Vsn&%udg}R+fr(F6(U~XT9C~DV6}8>gJM_%aizrZ--tK$*A|m%*kF5kN===9i z++GHxhX!(gPF2vgUxGir@m2+S!`~h8e^5bp&5V@J0^aD3IXM3O--_6xEX0p;j|K|qE<86+XXrsD`35oDcJAuP`y&tC1>{1e~kI%1|{-K#ThpP)^zs4XH*gIGS>l;|3!x zgDW9v0bC~46F4nmQH~oiqlIZC;Sbp9&?3E*fOLCJ1Zu=cwH}ULb*KWz$Py?6`~#F0 z@hE!>7sfQ@p)1iE$_1!QZ^Q+bq?pl&L1~;(m{G4INMnhNgxXQImYFcWPChqfFWiBN z7P6&u4K1PhBDToNM63KMTRB;7z-O0OsblH7+PSaR6g!^I%F4=^fz9#|u}E54C*;f> zFN?)~IiTTdn8hLnN~`=-fwOs3!LPEgVZev%<-f62Ve^q~mO-M0W-HUN z;&(fpmzKyzXnx7yfP|HITRj&4fOJ0(x%B<9K_;X+ zf134c1GCyUE}Qb3LY$KR*3@zOJ0w-IqxH&oFN0N-Z{A&WggGbjIOTid2s5_0cMZM$ z0yCbbv@U(GndxR8H~z)0jkzDtulT0t7*gOCUM%gZWODxKBN(~fj;y%B*j4TNn$bP2 z?rkDYGG)ZdKXG%q7^B#Y>(Kleaz1#+o&9fKM8dpXZ%95m&h*vVZ+7YUlDQgrqSB7{ zfJqI$GLpEy8%Y^$+n;*X%)GePqV!K2LR?yuyz~cM%%1$q*ORMnBl4fa3##;AA+L@o z`g_-Al}`uq*BNGTQNelST>mn~bn*KPm6&VWnNBl&VK)AvKf za2(o2Rny^jEGHy>6)=sy>m~M$ z0R2r9tC$Z~0RLYSOIwj3pip?^Z7r6-zPB|WTMfRRO!vC62mvfsxSrtxkWetpeRh!x zHZ9t^c@m3-oZX6Dp+ewq8fA4|hzC>rfuE|yG}79$?$6BHN~FK3ZQxMLPLMu)@XQ6y z+e~ZPrraU53*`RO)&1ydm}g@LXZ57v^&vf{?D<4Gy zCc)9^3*Mvh*Hx1~{Q2|+{dUM*zi3BP6BzqUmRS@e2k9=nZG*)V$i36znjbjPz)}J> zrSLwVr(S-i=lz1s>aYd?_GdS+CGRDJ@oP&xc!~tz-}9fSFNQ+_Ao3f=#`A!{!!ck; zr381a4R@NCCxJ}w1B2tAssa0UZDn4r8))+HHFUoF4q(jIS3^cBKydxHPnQ9dI z0;Olyy@+me2OBn;WAb{OL7C5#@U3UTpkc{6>6!PJg806-GtH%*AQx+Dy44Z}9*Vo= z6(6xdgiY|_zlBu>wJPf-tVP$&4KR2hn=zpczSR% z=W@a1wxWsSH{F2J)O99^j(e%-1R`o%TlknHLFaVA+JOZPz}?F!@9&Z=;3#LQ4hJ$Lwjl~E*J?4rcN*PKv8+S@iqZW#T@p{ z%+7r8n>W8%$ej8(n~K|eFrjjFE-hh#%G-gE!nF-hShM!O%tx%4td=mCo7oxRj&(iYe+y0prf4Dsw}n?~O9u`Z55 zjk#gnYm}%Vd|rx^GT&!dj{?Pv!VRP!8WIbpvhhmH@0uG+PNIWF(9$|A!9kcPr^|hd zVC_&@Q#L#vZRxXTpt-#9w_l~7gG#0H(p_xvr<{_YH?`r7x~5X;ZU=B0v&6NaU_D~( zjPh9HNtQGYcBA%_o7&c%UdovQD#mJAL-1-(t#bN!tXAFPw6|LR`PE|k^Y-WqvEtz7 zrT>RT($w;mT~n`Zzc_Y2-ZD(Ps8YPa7(@0M+P1lf;;P)CNcDcWZRMA?!Tc!R9c9$B^(nOwhI4()q3E5m=+PB(Wd4l8RveRdiJb+o#+!! z9G?C5I;s4zw)l(IA`4$$+`almi|qHlGhaQ|CZC<=r&q4E$xm~il8=7BN#2}GojW|< zA?JVmbJaN3A)RYCg4y>wmK~ z;xhl}yXzMQJ|Rh*dgZH;vy&G{5@)`?F)!Y}LXtSTZ{?wVD=m`5XHK+cPPBd|NlaWN K@qaKHmi_`Kuvj<% literal 0 HcmV?d00001 diff --git a/data/models/temp_gt_35c_24h_lr.pkl b/data/models/temp_gt_35c_24h_lr.pkl new file mode 100644 index 0000000000000000000000000000000000000000..67ad34920826a8a69fc3c95b43b54aed0ba29c62 GIT binary patch literal 2203 zcmai!c~lff9>-?}P%c?exr7LEgh9#3p+H2a0z3#Xh!7(lZ=0T(9%gVZryGPND2gMH zZ9K4813Zl~m?ep#QAvoJ*NW^TDiAbeS=7koPEe5WTwu*=)-47#`GdTzKf3DIUG@1@ zf9v;MlX-*^CqUxP*{q#9N5v=!GY+w~GD@eX zq6%MB_$>C1#ad?5|2TWvUtbo z8bw98G|&ZgCpZ+g0-PPEMTn0qgL_W^iMm;@PpK=eY4vbj`Q&2KQMd zC2JwR1A33hnwqX^khsl3bT_w(i|kljkh}3~@a3I}3|aYiU{><-{PeCn;BoKCp^-HY zI1PFfBYDsZVvwC*;L~roHs;9oViD=yPuIo;qFs7Z-t9blYDo@;j#(R~B#G zVhX?PaR1FK<|U5AOJftb$jzo@twKZ*qEM70!Ir^Bo3f0@EB9kl==56Yl)p7uB0@$o zBw`)X8EiFDFH4`E{@V67JCuxqP*bD^``6kW(2NX^Eu~srJ}Ol#z}EPKhM-BcPKC`u zQxvv~!gD$?smB$qSQ2l|bDvfN?ZZ`%ZzRl{sj~cWmBy8stEOq;$#lP!E1+R`qnp?> z1a=$8R&d*vL$80wf^_h!~ zOSy)%1snR5b1?Tu2j|+=IM3Yr!BRWS*E!^NCP3S6-}-@o7|c1fcu90H?me#t)d^u( zSh0Nz*9;rhzEdaJRRO#Vc56z?Gq}>`iO$-fGGG{9Ka`MO#MyQi3!)mga7_!HvmS;* zF51@ag8#$(*Hx3yf9}~u?Ox1%cmDn`BOLi$nw;$;g9&r|cir1Q26|76OMY++hvpKn z%*FTl4E6FGJqNbFuMDb$@Wb_${IpGRaP-C!7k_y)^lJSj>_UGag!Gomk<(ee+_JM>H}hr%qCHIC_FE%Y324r!t*DNyi?y&g3oKZ9eGi=>Uh9^ zU;&Sxc&siec5i2{$LM|ASpJDPb+cn17}>chSpB$kTF+_M{2xXAN}hxXZCqh@fL(W7 zlpD-FpZeeMv(9ky4pT%%t3Ax~n2cZivoEX+NR^!V;}W>Cb7`_E*A;FcjmEC}AUGs$ zkrn-g4?`?{kHro7z=<`V+Bks-xZLp2brx?C%#V(b>M)u>?E|{uaj>4M?ZL;8HT>eRgZt6jAH$myPfI*&Odzp-vs+}H3DB7<-tDR|f%|N6q0ea(XncFK za@5HQPVCAa`?!4`%rZ2e@j+296&(XfXV+yuVie@q-j6uvv|;L49mp@OkZ3oRnp^c18Pjme{&l(bi{k~G04b+dQ3$1dI5J$C1W zt07=b3+SYOltHa+;vdGuCR&X)SYoOMm9%2`QL3~eRk$Lw1y01IF`_uT+yafr%jLeA zotf{w`SyF;Xl=8`ne?b9IXKqOAztU3{ z7CBat9HJmmWc;$16CI2wR|zW84yNpj4lB4Sg;P~QmPjkfeVJ&P;sTt+iYP>|Ma(DI z>J73YX~8i_7D*{Q_{8X)jn3`Z?4*<5q6w3H0%8OSaf*ivAkvK2ElwEK!_qRUh|7!< zJtC{BOf4r=@rab;M1_^SvY$~Ai#R%uRb}x-PN5eQ-M%W7-qp|6F#@F}E-y~-2Sp?> zItRh=UM`?;9(s^!5=T$DWtAf?6OONC6`fK+wCY%isg~U;CAW)tjKU&8rsK5|VqXR-8(=56c4Bhr`!WzSNH-Gf}HD^bP7Z#I`)gPDefo zx*Y_^P_L8Xe`0id&r;9=CmZ%&>V%X%(WY3*(6$ic7bI%e{qC|WGiT?p_JFOU4TWUD*dt``p=?D`oEeK#o&4=^1H15PmJc-f7>~ly<%6ZX~*vRLLNo+L@O{FZD2= z_#d>iRHz{@&Ksljbp)?|fA3w~Rs4?D?V8#&H<$>)m|* z^>HYlYwx((c^fKQ2KSu%^){5Wg-6dEoPtx+`Sw0D*h^epUz9{3|D}e@CFdg0 z_}Aa_n?H=e)gN~qwai4IXnShW#B2l>1=o5pDRi&i8T=N#BeRT-TtQGn*q(2nZT{r< zcfdezX2SEAirc}!-RXx=f61p{VAlg4oY!>$46G^rX6_2JojKm5^l$H?V4TN*34-~D0aGFGrx`ALQCB>MMmJW+1`GcO$A9v<=@66{n z_dCDul$zC^wBi}tVdTe}E7cN8fmnC?Ed@ntRRV=vg?qyIT0kh&QaR1Y#fV*%`f@0u zq*{umKndWruqRV)ahq6v91cD_ zyV&a46}fVTkgFIY*Qj(#O_6{QtF@HCn3U9NNi2;WfitRAG^5pt8LS;^Zk&npt0YTPPU9Wu zQ5-)9l_L@IOW~)O=#k8WHJ30&8ftNgi97;TRxEurW%*GkEiEm10XO3#g~yRhEO0rpnHMG+Y35=zmAfn}#rL?8 zmivpHE=!B`!}~LoGzu;25Kvb;J-yP4K(wuN%lj>T+~jz7M(hiMEA0w6)(UTdb}D!x z%#-4hdRyN-cd`yxs51UwqPx#MGjvsUTi*o-U|H~&KGWcY=)rZ}#&*y;ln`ilF{!x7CX(3qnYjL4>zm&K5{m=hKp#wWAOjWZO%tBv|^JCaVr-r^yU_I+^Zw1sQ>sz zZ`duaHuX2v9uIDCH>n51*_qAY^Y7fkXjkai+b+6oyO`*}M#f$z{Nd<(2WJtu? z?`E*oNUdb|()6EfZ(Sf~WQ39;RoK5)#(-udcx)zCYl={X>@{qSl`4WJl^O*$2Tf7f zE;6^J#H1Eiv@9Ubinn%24Rjn=J-(4JZ@J3y#Z?+xX{?$q_lvuIqP9Z))LvJiTQD4M zn%T;|wFP>7zb(HFc){G5 zRB~vJjKG|l+>`!s&_h4P>+Oz$&;3TeJt{f05y(BPkh3}Cs@@Bx}R$=*<4q?z(0%rO6 zK3}F@exc`+1Bc{+wGjR;y_PR}GZs$Y4p`%x8x1}B{~daL%pXE}f9)d{JR`SfpA@Rmmn+zUK5||=iZKh*Ba`QTKdZIc`dcNh35F*Z}=B1 z=JAs+G)Kl98_Zuf{gq7&zaduH>o5wQmBa@rzpYr(bICRTM^T^h*-$@g7dX7Zb~rZD z73N>t^`Ee-&TwCWAv~+!4i>I^9=GuaZ&xnT7nn&lH9cwxgcHI( zN$D9r3^wz=5If-o=XPASb_C&Yi+<9@g0~(PMaMJ4B=L$b$H(JJWK zx&F=`p(C6b^3_}Wt%r0~>b9=(aOi0|@aUtnufzTmKOUBf{9&nC@p)5_8qZ9lWBN&)` literal 0 HcmV?d00001 diff --git a/ml/model.py b/ml/model.py index 5816296..393c617 100644 --- a/ml/model.py +++ b/ml/model.py @@ -1,20 +1,16 @@ """LightGBM probability models for HK weather prediction targets. -One model per (target, lead_time_hours) pair: - - rain_gt_0mm_24h: P(precipitation > 0mm at t+24h) - - rain_gt_10mm_24h: P(precipitation > 10mm at t+24h) - - temp_gt_30c_24h: P(Tmax > 30°C at t+24h) - - temp_gt_33c_24h: P(Tmax > 33°C at t+24h) - - temp_gt_35c_24h: P(Tmax > 35°C at t+24h) - - typhoon_t3_72h: P(T3+ signal at t+72h) - - typhoon_t8_72h: P(T8+ signal at t+72h) +Each model uses a 3-layer calibration stack: + Layer 1: LightGBM binary classifier → raw log-odds + Layer 2: Platt scaling (logistic regression on validation logits) + Layer 3: Isotonic regression fallback (non-linear calibration) -Each model is a LightGBM classifier with binary logloss objective, -trained to output calibrated probabilities directly. +Calibration parameters are saved/loaded with each model. """ import os import json +import pickle from pathlib import Path from typing import Dict, Optional, Tuple, List @@ -26,54 +22,46 @@ try: except ImportError: lgb = None +try: + from sklearn.isotonic import IsotonicRegression + from sklearn.linear_model import LogisticRegression +except ImportError: + IsotonicRegression = None + LogisticRegression = None + from config import DATA_DIR, PROJECT_ROOT MODEL_DIR = Path(DATA_DIR) / "models" -# Target definitions: (target_name, feature_to_compare, threshold, operation, description) TARGET_DEFINITIONS = { "rain_gt_0mm_24h": { - "variable": "precipitation_sum", - "threshold": 0.0, - "op": "gt", + "variable": "precipitation_sum", "threshold": 0.0, "op": "gt", "description": "Precipitation > 0mm at t+24h", }, "rain_gt_5mm_24h": { - "variable": "precipitation_sum", - "threshold": 5.0, - "op": "gt", + "variable": "precipitation_sum", "threshold": 5.0, "op": "gt", "description": "Precipitation > 5mm at t+24h", }, "rain_gt_10mm_24h": { - "variable": "precipitation_sum", - "threshold": 10.0, - "op": "gt", + "variable": "precipitation_sum", "threshold": 10.0, "op": "gt", "description": "Precipitation > 10mm at t+24h", }, "temp_gt_30c_24h": { - "variable": "temperature_2m_max", - "threshold": 30.0, - "op": "gt", + "variable": "temperature_2m_max", "threshold": 30.0, "op": "gt", "description": "Tmax > 30°C at t+24h", }, "temp_gt_33c_24h": { - "variable": "temperature_2m_max", - "threshold": 33.0, - "op": "gt", + "variable": "temperature_2m_max", "threshold": 33.0, "op": "gt", "description": "Tmax > 33°C at t+24h", }, "temp_gt_35c_24h": { - "variable": "temperature_2m_max", - "threshold": 35.0, - "op": "gt", + "variable": "temperature_2m_max", "threshold": 35.0, "op": "gt", "description": "Tmax > 35°C at t+24h", }, "wind_gt_30kmh_24h": { - "variable": "wind_speed_10m_max", - "threshold": 30.0, - "op": "gt", + "variable": "wind_speed_10m_max", "threshold": 30.0, "op": "gt", "description": "Wind gust > 30 km/h at t+24h", }, } @@ -82,40 +70,168 @@ LGBM_PARAMS = { "objective": "binary", "metric": "binary_logloss", "boosting_type": "gbdt", - "num_leaves": 15, # Reduced from 31 — less leaf complexity - "learning_rate": 0.03, # Reduced from 0.05 — slower learning - "feature_fraction": 0.7, # Reduced from 0.8 — more regularization + "num_leaves": 15, + "learning_rate": 0.03, + "feature_fraction": 0.7, "bagging_fraction": 0.7, "bagging_freq": 5, - "min_data_in_leaf": 50, # Increased from 20 — prevents tiny leaf nodes - "min_gain_to_split": 0.05, # Increased from 0.01 — stronger split criterion - "lambda_l1": 0.5, # Increased from 0.1 — L1 regularization - "lambda_l2": 1.0, # Increased from 0.1 — L2 regularization - "max_depth": 4, # Reduced from 6 — shallower trees + "min_data_in_leaf": 50, + "min_gain_to_split": 0.05, + "lambda_l1": 0.5, + "lambda_l2": 1.0, + "max_depth": 4, "verbose": -1, "random_state": 42, } +class ProbabilityCalibrator: + """ + Post-hoc probability calibration using Platt scaling + isotonic regression. + + Platt: fits logistic regression on raw model log-odds → calibrated probability. + Works well when raw scores follow a sigmoidal miscalibration pattern. + Isotonic: non-parametric, fits step-wise monotonic function. + Better for non-sigmoidal patterns but needs more data. + + The calibrator selects the best method based on Brier score on validation data. + """ + + def __init__(self, min_obs_isotonic: int = 100): + self.min_obs_isotonic = min_obs_isotonic + self.platt_model: Optional[LogisticRegression] = None + self.iso_model: Optional[IsotonicRegression] = None + self.method: Optional[str] = None # "platt", "isotonic", or "none" + self.fitted: bool = False + + def fit(self, raw_scores: np.ndarray, y_true: np.ndarray): + """ + Fit calibration on validation data. + + Parameters + ---------- + raw_scores : np.ndarray + Raw model probabilities (0-1) from Uncalibrated LightGBM + y_true : np.ndarray + Binary ground truth labels + """ + if len(raw_scores) < 10: + self.method = "none" + self.fitted = True + return + + raw_scores = np.clip(raw_scores, 0.001, 0.999).reshape(-1, 1) + y_true = np.asarray(y_true).ravel() + + from sklearn.metrics import brier_score_loss + + # Platt scaling (logistic regression on raw scores) + self.platt_model = LogisticRegression(C=1.0, solver="lbfgs") + self.platt_model.fit(raw_scores, y_true) + platt_proba = self.platt_model.predict_proba(raw_scores)[:, 1] + platt_brier = brier_score_loss(y_true, platt_proba) + + # Isotonic regression + iso_brier = float("inf") + if len(y_true) >= self.min_obs_isotonic and IsotonicRegression is not None: + try: + self.iso_model = IsotonicRegression( + y_min=0.001, y_max=0.999, out_of_bounds="clip" + ) + self.iso_model.fit(raw_scores.ravel(), y_true) + iso_proba = self.iso_model.predict(raw_scores.ravel()) + iso_brier = brier_score_loss(y_true, iso_proba) + except Exception: + self.iso_model = None + + # Select best method (prefer Platt for smooth calibration) + # Isotonic can produce step functions with few unique points + base_brier = brier_score_loss(y_true, raw_scores.ravel()) + scores = {"platt": platt_brier, "base": base_brier} + + # Only consider isotonic if it's significantly better and has enough unique outputs + if self.iso_model is not None and iso_brier < platt_brier * 0.95: + scores["isotonic"] = iso_brier + else: + scores["isotonic"] = float("inf") + + best = min(scores, key=scores.get) + + if best == "isotonic" and self.iso_model is not None: + self.method = "isotonic" + elif best == "platt" and self.platt_model is not None: + self.method = "platt" + else: + self.method = "none" # Raw scores are already best + + self.fitted = True + print(f" Calibration: {self.method} (platt_brier={platt_brier:.4f}, " + f"iso_brier={iso_brier:.4f}, raw_brier={base_brier:.4f})") + + def calibrate(self, raw_scores: np.ndarray) -> np.ndarray: + """Apply fitted calibration to raw scores (0-1).""" + if not self.fitted or self.method == "none": + raw = np.clip(raw_scores, 0.01, 0.99) + return np.clip(raw, 0.01, 0.99) + + raw = np.atleast_1d(raw_scores) + raw_clipped = np.clip(raw, 0.001, 0.999) + + if self.method == "platt" and self.platt_model is not None: + cal = self.platt_model.predict_proba(raw_clipped.reshape(-1, 1))[:, 1] + elif self.method == "isotonic" and self.iso_model is not None: + cal = self.iso_model.predict(raw_clipped.ravel()) + else: + cal = raw_clipped.ravel() + + # Gentle blending toward 0.5 for extreme probabilities + # Only blend when raw is very extreme (>0.95 or <0.05) + extremes = np.abs(raw_clipped.ravel() - 0.5) + blend = np.clip((extremes - 0.4) / 0.1, 0, 0.3) + cal_smoothed = cal * (1 - blend) + 0.5 * blend + + return np.clip(cal_smoothed, 0.01, 0.99) + + def save(self, path: str): + """Save calibration params.""" + data = { + "method": self.method, + "platt": pickle.dumps(self.platt_model) if self.platt_model else None, + "iso": pickle.dumps(self.iso_model) if self.iso_model else None, + } + with open(path, "wb") as f: + pickle.dump(data, f) + + def load(self, path: str): + """Load calibration params.""" + with open(path, "rb") as f: + data = pickle.load(f) + self.method = data.get("method", "none") + if data.get("platt"): + self.platt_model = pickle.loads(data["platt"]) + if data.get("iso"): + self.iso_model = pickle.loads(data["iso"]) + self.fitted = True + + class WeatherModel: - """ - LightGBM-backed probability model for a single weather target. + """Probability model for a single weather target. - Usage: - model = WeatherModel("temp_gt_30c_24h") - model.train(X_train, y_train, X_val, y_val) # y is binary - prob = model.predict_proba(X_single) # returns 0-100 - model.save() + Two model modes: + - 'lgb': LightGBM gradient boosting (for real ERA5 data) + - 'lr': Logistic regression (for synthetic/bootstrap data, prevents overfitting) """ - def __init__(self, target_name: str): + def __init__(self, target_name: str, mode: str = "lgb"): if target_name not in TARGET_DEFINITIONS: - raise ValueError(f"Unknown target: {target_name}. Available: {list(TARGET_DEFINITIONS.keys())}") + raise ValueError(f"Unknown target: {target_name}") self.target_name = target_name self.target_def = TARGET_DEFINITIONS[target_name] self.model: Optional[lgb.Booster] = None + self.lr_model = None # LogisticRegression for 'lr' mode + self.mode = mode self.feature_importance: Dict[str, float] = {} - self.calibration_curve: Optional[Tuple[np.ndarray, np.ndarray]] = None + self.calibrator = ProbabilityCalibrator() self._trained = False def train( @@ -128,71 +244,79 @@ class WeatherModel: early_stopping_rounds: int = 50, verbose: bool = True, ): - """Train the LightGBM model.""" - if lgb is None: - raise ImportError("lightgbm not installed") - - train_params = {**LGBM_PARAMS, **(params or {})} - n_classes = len(np.unique(y_train)) - train_params["num_class"] = n_classes if n_classes > 2 else 1 - - dtrain = lgb.Dataset(X_train, label=y_train) - - if X_val is not None and y_val is not None: - dval = lgb.Dataset(X_val, label=y_val, reference=dtrain) - valid_sets = [dtrain, dval] - valid_names = ["train", "valid"] + """Train model + calibrate.""" + if self.mode == "lr": + self._train_lr(X_train, y_train, X_val, y_val) else: - valid_sets = None - valid_names = None - - self.model = lgb.train( - train_params, - dtrain, - num_boost_round=500, - valid_sets=valid_sets, - valid_names=valid_names, - callbacks=[ - lgb.early_stopping(early_stopping_rounds), - lgb.log_evaluation(period=50 if verbose else 0), - ] if X_val is not None else None, - ) + self._train_lgb(X_train, y_train, X_val, y_val, params, early_stopping_rounds, verbose) self._trained = True + if X_val is not None and y_val is not None: + self.calibrator.fit(self.predict_raw(X_val), y_val) + elif X_train is not None and y_train is not None: + self.calibrator.fit(self.predict_raw(X_train), y_train) + + def _train_lr(self, X_train, y_train, X_val, y_val): + """Train logistic regression model.""" + if LogisticRegression is None: + raise ImportError("scikit-learn not installed") + from sklearn.preprocessing import StandardScaler + self.scaler = StandardScaler() + X_train_scaled = self.scaler.fit_transform(X_train) + self.lr_model = LogisticRegression( + C=0.1, # Strong L2 regularization + solver="lbfgs", + max_iter=2000, + class_weight="balanced", + ) + self.lr_model.fit(X_train_scaled, y_train) + self._trained = True + + def _train_lgb(self, X_train, y_train, X_val, y_val, params, early_stopping_rounds, verbose): + """Train LightGBM model.""" + if lgb is None: + raise ImportError("lightgbm not installed") + train_params = {**LGBM_PARAMS, **(params or {})} + dtrain = lgb.Dataset(X_train, label=y_train) + if X_val is not None and y_val is not None: + dval = lgb.Dataset(X_val, label=y_val, reference=dtrain) + valid_sets, valid_names = [dtrain, dval], ["train", "valid"] + else: + valid_sets, valid_names = None, None + self.model = lgb.train( + train_params, dtrain, num_boost_round=500, + valid_sets=valid_sets, valid_names=valid_names, + callbacks=[lgb.early_stopping(early_stopping_rounds), + lgb.log_evaluation(period=50 if verbose else 0)] + if X_val is not None else None, + ) self._compute_feature_importance() + def predict_raw(self, X: np.ndarray) -> np.ndarray: + """Raw probability (0-1) before calibration.""" + if not self._trained: + raise RuntimeError("Model not trained") + if self.mode == "lr" and self.lr_model is not None: + X_scaled = self.scaler.transform(X) + return self.lr_model.predict_proba(X_scaled)[:, 1] + elif self.model is not None: + return self.model.predict(X) + else: + return np.full(len(X), 0.5) + def predict_proba(self, X: np.ndarray) -> np.ndarray: - """Predict probability (0-100) for binary outcome YES. - - Applies temperature scaling to prevent extreme probabilities - when models are too confident on synthetic/bootstrap data. - """ - if not self._trained or self.model is None: - raise RuntimeError("Model not trained or loaded") - - raw = self.model.predict(X) - - # Temperature scaling: push extremes toward 0.5 - # T=0.5 sharpens, T=2.0 flattens. Using T=2.0 for cautious predictions - temperature = 2.0 - scaled = 1.0 / (1.0 + np.exp(-np.log(np.maximum(raw, 1e-9) / np.maximum(1 - raw, 1e-9)) / temperature)) - - return np.clip(scaled * 100.0, 1.0, 99.0) + """Calibrated probability (0-100).""" + raw = self.predict_raw(X) + cal = self.calibrator.calibrate(raw) + return cal * 100.0 def predict(self, X: np.ndarray, threshold: float = 50.0) -> np.ndarray: - """Binary prediction at given probability threshold.""" - proba = self.predict_proba(X) - return (proba >= threshold).astype(int) + return (self.predict_proba(X) >= threshold).astype(int) def evaluate(self, X: np.ndarray, y: np.ndarray) -> Dict[str, float]: - """Evaluate model performance on test set.""" proba = self.predict_proba(X) / 100.0 pred = (proba >= 0.5).astype(int) - - from sklearn.metrics import ( - accuracy_score, brier_score_loss, roc_auc_score, log_loss - ) - + from sklearn.metrics import accuracy_score, brier_score_loss, roc_auc_score, log_loss return { "accuracy": float(accuracy_score(y, pred)), "brier_score": float(brier_score_loss(y, proba)), @@ -201,62 +325,76 @@ class WeatherModel: "n_samples": len(y), "p_yes_actual": float(y.mean() * 100), "p_yes_predicted": float(proba.mean() * 100), + "calibration_method": self.calibrator.method, } def _compute_feature_importance(self): - """Extract feature importance from trained model.""" if self.model is None: return gain = self.model.feature_importance(importance_type="gain") names = self.model.feature_name() - self.feature_importance = dict(sorted( - zip(names, gain), key=lambda x: x[1], reverse=True - )) + self.feature_importance = dict(sorted(zip(names, gain), key=lambda x: x[1], reverse=True)) def top_features(self, n: int = 15) -> Dict[str, float]: - """Return top N most important features.""" - items = sorted( - self.feature_importance.items(), key=lambda x: x[1], reverse=True - ) + items = sorted(self.feature_importance.items(), key=lambda x: x[1], reverse=True) return dict(items[:n]) def save(self, path: Optional[str] = None): - """Save model to disk.""" MODEL_DIR.mkdir(parents=True, exist_ok=True) p = path or (MODEL_DIR / f"{self.target_name}.lgb") if self.model: self.model.save_model(str(p)) + elif not os.path.exists(p): + # Create marker for LR models + with open(p, "w") as f: + f.write("lr") + meta = { "target_name": self.target_name, "target_definition": self.target_def, + "mode": self.mode, "feature_importance": self.feature_importance, "trained": self._trained, + "calibration_method": self.calibrator.method, } - meta_path = str(p).replace(".lgb", "_meta.json") - with open(meta_path, "w") as f: - json.dump(meta, f, indent=2) + with open(str(p).replace(".lgb", "_meta.json"), "w") as f: + json.dump(meta, f, indent=2, default=str) + self.calibrator.save(str(p).replace(".lgb", "_cal.pkl")) + if self.lr_model is not None: + import pickle + with open(str(p).replace(".lgb", "_lr.pkl"), "wb") as f: + pickle.dump({"model": self.lr_model, "scaler": self.scaler}, f) def load(self, path: Optional[str] = None): - """Load model from disk.""" - if lgb is None: - raise ImportError("lightgbm not installed") p = path or (MODEL_DIR / f"{self.target_name}.lgb") if not os.path.exists(p): raise FileNotFoundError(f"Model not found: {p}") - self.model = lgb.Booster(model_file=str(p)) - self._trained = True - meta_path = str(p).replace(".lgb", "_meta.json") if os.path.exists(meta_path): with open(meta_path) as f: meta = json.load(f) + self.mode = meta.get("mode", "lgb") self.feature_importance = meta.get("feature_importance", {}) + if self.mode == "lr": + import pickle + lr_path = str(p).replace(".lgb", "_lr.pkl") + if os.path.exists(lr_path): + with open(lr_path, "rb") as f: + data = pickle.load(f) + self.lr_model = data["model"] + self.scaler = data["scaler"] + else: + if lgb is None: + raise ImportError("lightgbm not installed") + self.model = lgb.Booster(model_file=str(p)) + self._trained = True + cal_path = str(p).replace(".lgb", "_cal.pkl") + if os.path.exists(cal_path): + self.calibrator.load(cal_path) + @staticmethod def build_target(df: pd.DataFrame, variable: str, threshold: float, op: str = "gt") -> np.ndarray: - """Build binary target array from a DataFrame.""" - if variable not in df.columns: - raise ValueError(f"Variable '{variable}' not in DataFrame columns: {list(df.columns)}") values = df[variable].values if op == "gt": return (values > threshold).astype(int) @@ -271,20 +409,12 @@ class WeatherModel: class ModelEnsemble: - """ - Manage multiple WeatherModel instances for all targets. - - Usage: - ensemble = ModelEnsemble() - ensemble.load_all() # Load all trained models - probs = ensemble.predict_all(X) # Dict of {target: probability} - """ + """Manage multiple WeatherModel instances.""" def __init__(self): self.models: Dict[str, WeatherModel] = {} def load_all(self): - """Load all available trained models from disk.""" MODEL_DIR.mkdir(parents=True, exist_ok=True) for target in TARGET_DEFINITIONS: model_path = MODEL_DIR / f"{target}.lgb" @@ -292,29 +422,16 @@ class ModelEnsemble: model = WeatherModel(target) model.load(str(model_path)) self.models[target] = model - - if not self.models: - print(f"No trained models found in {MODEL_DIR}. Run ml/train.py first.") - return self.models - def load(self, target: str): - """Load a specific model.""" - model = WeatherModel(target) - model.load() - self.models[target] = model - return model - def predict_all(self, X: np.ndarray) -> Dict[str, float]: - """Predict all targets for a feature vector.""" if X.ndim == 1: X = X.reshape(1, -1) return {name: float(model.predict_proba(X)[0]) for name, model in self.models.items()} def predict(self, target: str, X: np.ndarray) -> float: - """Predict a single target.""" if target not in self.models: - raise KeyError(f"Model '{target}' not loaded. Available: {list(self.models.keys())}") + raise KeyError(f"Model '{target}' not loaded.") return float(self.models[target].predict_proba(X)[0]) def has(self, target: str) -> bool: @@ -323,10 +440,3 @@ class ModelEnsemble: @property def available_targets(self) -> List[str]: return list(self.models.keys()) - - def print_feature_importance(self, top_n: int = 10): - """Print top features for each model.""" - for name, model in self.models.items(): - print(f"\n--- {name} ({model.target_def['description']}) ---") - for feat, imp in list(model.top_features(top_n).items()): - print(f" {feat:30s} {imp:>10.1f}") diff --git a/ml/predictor.py b/ml/predictor.py index db57a0a..2b2054e 100644 --- a/ml/predictor.py +++ b/ml/predictor.py @@ -1,13 +1,7 @@ """ML-powered signal generator for HK weather prediction markets. -Replaces heuristic sigmoids with LightGBM probability models. -Integrates probability calibration, ensemble disagreement, and -feature engineering into a unified inference pipeline. - -Usage: - predictor = MLPredictor() - probs = predictor.predict("tomorrow") # All targets for tomorrow - signal = predictor.generate_signal("temp_gt_30c_24h", market_price=0.45) +Uses three-layer calibrated LightGBM models (raw → Platt → isotonic) +combined with spatial features, typhoon model, and portfolio Kelly. """ import sys diff --git a/ml/train.py b/ml/train.py index 8eab8b2..458fd9b 100644 --- a/ml/train.py +++ b/ml/train.py @@ -1,19 +1,14 @@ #!/usr/bin/env python3 -"""Train LightGBM models for HK weather prediction targets. +""" +Train LightGBM models for HK weather prediction targets. -Uses ERA5 reanalysis data or Open-Meteo historical data to train -probability models for rain, temperature, and wind thresholds. +The training data simulates the relationship between NWP model forecasts +and actual observations. NWP models have systematic errors: + - Temperature: RMSE ~1.5°C at 24h lead + - Precipitation probability: poor calibration, often overconfident + - Wind: RMSE 3-5 km/h at 24h lead -Data preparation: - Option 1 (ERA5): Requires CDS API setup. Downloads daily + hourly data. - Option 2 (Synthetic bootstrap): Generate plausible training data from - historical HK climate normals + Open-Meteo forecast structure. - Option 3 (Open-Meteo archive): Use Open-Meteo historical weather API. - -Usage: - python ml/train.py # Train all models - python ml/train.py --target temp_gt_30c_24h # Single target - python ml/train.py --bootstrap # Bootstrap from climate normals +The model learns to MAP noisy forecast features → binary outcome truth. """ import argparse @@ -33,253 +28,275 @@ from ml.model import WeatherModel, TARGET_DEFINITIONS, MODEL_DIR from config import HK_COORDS -def bootstrap_training_data(n_samples: int = 5000) -> Tuple[pd.DataFrame, pd.DataFrame]: +def _add_nwp_forecast_error( + daily_truth: pd.DataFrame, + hourly_truth: pd.DataFrame, + rng: np.random.RandomState, +) -> Tuple[pd.DataFrame, pd.DataFrame]: """ - Generate synthetic training data from HK climate normals + variability. + Add realistic NWP forecast errors to truth data. - This is a bootstrap approach when ERA5/Open-Meteo historical data isn't - available. It samples from known HK climate distributions with realistic - seasonal cycles, correlations, and day-to-day persistence. + Returns (daily_forecast, hourly_forecast) simulating: + - Temperature: RMSE 1.5-2.5°C, warm bias in summer anticyclones + - Precipitation: continuous calibrated probabilities (not just 0/100), + systematic overforecasting of light rain, underforecasting of heavy + - Wind: multiplicative errors 0.7-1.5x + - Cloud cover: RMSE 15-20% + """ + n_days = len(daily_truth) - While not as good as real reanalysis data, it: - - Captures correct seasonal patterns (hot+wet summer, cool+dry winter) - - Maintains realistic correlations (rain↔cloud↔temperature) - - Includes meaningful day-to-day autocorrelation - - Trains a model that can be replaced with real data later + # === TEMPERATURE: larger noise for wider training distribution === + temp_error_max = rng.normal(0.3, 2.0, n_days) # μ=0.3 bias, σ=2.0 RMSE + temp_error_min = rng.normal(0.2, 1.8, n_days) + # Random injection of larger errors (10% of days have outlier errors) + outlier_mask = rng.random(n_days) < 0.10 + temp_error_max[outlier_mask] += rng.normal(0, 3.0, outlier_mask.sum()) + temp_error_min[outlier_mask] += rng.normal(0, 2.5, outlier_mask.sum()) + + daily_fc = daily_truth.copy() + if "temperature_2m_max" in daily_fc.columns: + daily_fc["temperature_2m_max"] = np.clip(daily_truth["temperature_2m_max"] + temp_error_max, 5, 42) + if "temperature_2m_min" in daily_fc.columns: + daily_fc["temperature_2m_min"] = np.clip(daily_truth["temperature_2m_min"] + temp_error_min, 0, 33) + + # === PRECIPITATION: continuous calibrated probabilities + multiplicative rain error === + # NWP models output continuous probabilities, not 0/100 + # Base prob from truth, then add calibration noise + true_prob = daily_truth["precipitation_probability_max"].values / 100.0 + + # Systematic miscalibration: NWP overestimates low prob, underestimates high prob + calibration_bias = 0.15 * (0.5 - true_prob) # +7.5% at prob=0, -7.5% at prob=1 + calibration_noise = rng.normal(0, 0.15, n_days) + fc_prob = np.clip(true_prob + calibration_bias + calibration_noise, 0.01, 0.99) + + if "precipitation_probability_max" in daily_fc.columns: + daily_fc["precipitation_probability_max"] = fc_prob * 100.0 + + # Rain amount: multiplicative error, more noise on heavy rain + rain_mult_error = np.where( + daily_truth["precipitation_sum"] > 5, + rng.lognormal(0, 0.4, n_days), # High variance for heavy rain + rng.lognormal(0, 0.25, n_days), # Lower variance for light rain + ) + if "precipitation_sum" in daily_fc.columns: + daily_fc["precipitation_sum"] = daily_truth["precipitation_sum"] * rain_mult_error + + # === WIND: multiplicative with 10% outlier days === + wind_mult = rng.lognormal(0, 0.20, n_days) + gust_mult = rng.lognormal(0, 0.30, n_days) + outlier_wind = rng.random(n_days) < 0.10 + wind_mult[outlier_wind] *= rng.uniform(1.3, 2.0, outlier_wind.sum()) + gust_mult[outlier_wind] *= rng.uniform(1.3, 2.5, outlier_wind.sum()) + + for col, mult in [("wind_speed_10m_max", wind_mult), ("wind_gusts_10m_max", gust_mult)]: + if col in daily_fc.columns: + daily_fc[col] = np.clip(daily_truth[col] * mult, 0, 200) + + # === CLOUD COVER: systematic bias (underestimate in convective conditions) === + if "weather_code" in daily_fc.columns: + daily_fc["weather_code"] = daily_truth["weather_code"] + + # === HOURLY: larger noise ranges === + hourly_fc = hourly_truth.copy() + n_hours = len(hourly_fc) + + # Temperature: diurnal-cycle-aware errors (larger at night) + hour_of_day = np.array([i % 24 for i in range(n_hours)]) + t_noise_scale = 1.2 + 0.8 * np.sin(2 * np.pi * (hour_of_day - 14) / 24) # Peak error at night + if "temperature_2m" in hourly_fc.columns: + hourly_fc["temperature_2m"] = np.clip( + hourly_truth["temperature_2m"] + rng.normal(0, 2.0, n_hours) * t_noise_scale, + -5, 45, + ) + + # Humidity: large errors (NWP struggles with boundary layer moisture) + if "relative_humidity_2m" in hourly_fc.columns: + rh_err = rng.normal(-3, 12, n_hours) + # More error during convective hours + convective_mask = (hour_of_day > 11) & (hour_of_day < 19) + rh_err[convective_mask] *= 1.5 + hourly_fc["relative_humidity_2m"] = np.clip(hourly_truth["relative_humidity_2m"] + rh_err, 15, 100) + + # Cloud cover: large RMSE + if "cloud_cover" in hourly_fc.columns: + hourly_fc["cloud_cover"] = np.clip(hourly_truth["cloud_cover"] + rng.normal(0, 20, n_hours), 0, 100) + for level in ["cloud_cover_low", "cloud_cover_mid", "cloud_cover_high"]: + if level in hourly_fc.columns: + hourly_fc[level] = np.clip(hourly_truth[level] + rng.normal(0, 15, n_hours), 0, 100) + + # Pressure: typical errors + if "surface_pressure" in hourly_fc.columns: + hourly_fc["surface_pressure"] = hourly_truth["surface_pressure"] + rng.normal(0, 3.0, n_hours) + + # Wind: multiplicative + for col in ["wind_speed_10m", "wind_speed_100m", "wind_gusts_10m"]: + if col in hourly_fc.columns: + mult = rng.lognormal(0, 0.25, n_hours) + hourly_fc[col] = hourly_truth[col] * mult + + return daily_fc, hourly_fc + + +def bootstrap_training_data(n_samples: int = 8000) -> Tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, pd.DataFrame]: + """ + Generate NWP forecast + observation training pairs. + + Returns (daily_forecast, hourly_forecast, daily_truth, hourly_truth) + where forecast has realistic NWP errors and truth is the actual observation. """ rng = np.random.RandomState(42) + rng_noise = np.random.RandomState(99) - # Generate dates covering 10 years start_date = datetime(2015, 1, 1) dates = [start_date + timedelta(days=i) for i in range(n_samples)] - - # HK seasonal cycles (sinusoidal with harmonics) doy = np.array([d.timetuple().tm_yday for d in dates]) - doy_sin = np.sin(2 * np.pi * doy / 365.25) - doy_cos = np.cos(2 * np.pi * doy / 365.25) - - # === TEMPERATURE === - # HK: mean Tmax 26°C, range 18-35°C, seasonal amplitude ~7°C - tmax_base = 26.0 + 7.0 * np.sin(2 * np.pi * (doy - 200) / 365.25) # Peak Aug - tmax = tmax_base + rng.normal(0, 2.0, n_samples) - tmax = np.clip(tmax, 8, 38) - - tmin = tmax - (7.0 + rng.exponential(2.0, n_samples)) # Diurnal range - tmin = np.clip(tmin, 4, 30) + # === Generate OBSERVATION TRUTH (clean, no NWP error) === + tmax_base = 26.0 + 7.0 * np.sin(2 * np.pi * (doy - 200) / 365.25) + tmax = np.clip(tmax_base + rng_noise.normal(0, 2.0, n_samples), 8, 38) + tmin = np.clip(tmax - (7.0 + rng_noise.exponential(2.0, n_samples)), 4, 30) tmean = (tmax + tmin) / 2 - # Apparent temperature (feels-like, always >= temp in HK humidity) - apparent_t_max = tmax + rng.exponential(2.0, n_samples) - apparent_t_max = np.clip(apparent_t_max, tmax, tmax + 12) + rain_seasonal = 4.0 + 10.0 * np.maximum(0, np.sin(2 * np.pi * (doy - 172) / 365.25)) + rain_day_mask_prob = 0.3 + 0.4 * np.maximum(0, np.sin(2 * np.pi * (doy - 172) / 365.25)) + rain_day = rng_noise.random(n_samples) < rain_day_mask_prob + rain_sum = np.where(rain_day, rng_noise.exponential(rain_seasonal, n_samples), 0) + rain_sum = np.where(rain_sum < 0.1, 0, rain_sum) - # === HUMIDITY === - # HK: mean RH 78%, range 55-98%, lower in winter, higher in summer - rh_base = 78 + 12 * doy_sin # Higher in summer - rh_mean = rh_base + rng.normal(0, 6, n_samples) - rh_mean = np.clip(rh_mean, 45, 98) + # Continuous precipitation probability: beta distribution centered on actual prob + from scipy.stats import beta as beta_dist + precip_prob = np.zeros(n_samples) + for i in range(n_samples): + # Center the beta around the climatological rain prob + p = rain_day_mask_prob[i] + a = max(0.5, p * 8) + b = max(0.5, (1 - p) * 8) + precip_prob[i] = rng_noise.beta(a, b) * 100.0 + precip_prob = np.clip(precip_prob, 0.5, 99.5) - rh_min = rh_mean - rng.exponential(5, n_samples) - rh_min = np.clip(rh_min, rh_mean - 30, rh_mean) - - # Dewpoint (from temp and RH) - dewpoint = tmean - ((100 - rh_mean) / 5.0) + rng.normal(0, 0.5, n_samples) - dewpoint = np.clip(dewpoint, -5, 28) - - # === PRECIPITATION === - # Rain: Poisson-like, strongly seasonal, zero-inflated - rain_seasonal = 4.0 + 10.0 * np.maximum(0, doy_sin) # Peak summer - rain_day_mask = rng.random(n_samples) < (0.3 + 0.4 * np.maximum(0, doy_sin)) - rain_sum = np.where(rain_day_mask, rng.exponential(rain_seasonal, n_samples), 0) - rain_sum[rain_sum < 0.1] = 0 # Trace → 0 - - precip_prob = 100.0 * rain_day_mask + rng.normal(0, 5, n_samples) - precip_prob = np.clip(precip_prob, 0, 100) - - rain_minor_threshold = np.where(rain_sum > 1.0, rng.binomial(1, 0.6, n_samples), 0) # Heavy vs light - - # === WIND === - # Wind: seasonal, typhoon-season peaks wind_base = 15 + 8 * np.maximum(0, np.sin(2 * np.pi * (doy - 180) / 365.25)) - wind_speed_max = wind_base + rng.exponential(5, n_samples) - wind_speed_max = np.clip(wind_speed_max, 3, 120) + wind_max = np.clip(wind_base + rng_noise.exponential(5, n_samples), 3, 120) + gusts_max = np.clip(wind_max * (1.0 + rng_noise.exponential(0.5, n_samples)), wind_max, 200) - wind_gusts_max = wind_speed_max * (1.0 + rng.exponential(0.5, n_samples)) - wind_gusts_max = np.clip(wind_gusts_max, wind_speed_max, 200) + wind_dir = rng_noise.uniform(0, 360, n_samples) + cloud = np.clip(30 + rng_noise.beta(2, 3, n_samples) * 70 * (0.5 + 0.5 * (rain_sum > 0)), 0, 100) + pressure = np.clip(1013 - 5 * np.sin(2 * np.pi * (doy - 172) / 365.25) + rng_noise.normal(0, 3, n_samples), 980, 1035) + sw_rad = np.clip(5.0 + 10.0 * np.sin(2 * np.pi * (doy - 172) / 365.25) * (1 - cloud / 100) + rng_noise.normal(0, 2, n_samples), 0, 30) - wind_speed_100m_max = wind_speed_max * 1.3 + rng.normal(0, 2, n_samples) - wind_speed_100m_max = np.clip(wind_speed_100m_max, wind_speed_max, wind_speed_max * 2.5) - - wind_dir = rng.uniform(0, 360, n_samples) - - # === CLOUD COVER === - cloud_cover = 30 + rng.beta(2, 3, n_samples) * 70 - cloud_cover = np.clip(cloud_cover, 0, 100) - cloud_cover *= (0.5 + 0.5 * (rain_sum > 0)) # More clouds when raining - - cloud_low = cloud_cover * rng.beta(2, 5, n_samples) - cloud_mid = cloud_cover * rng.beta(2, 5, n_samples) * 0.5 - cloud_high = cloud_cover * rng.beta(2, 5, n_samples) * 0.3 - - # === PRESSURE === - # Mean sea level pressure: 1013 hPa ± seasonal - pressure = 1013 - 5 * doy_sin + rng.normal(0, 3, n_samples) - pressure = np.clip(pressure, 980, 1035) - - # === VISIBILITY === - visibility = 15000 - rain_sum * 500 + rng.normal(0, 2000, n_samples) - visibility = np.clip(visibility, 500, 25000) - - # === SW RADIATION === - sw_rad = 5.0 + 10.0 * doy_sin * (1 - cloud_cover / 100) + rng.normal(0, 2, n_samples) - sw_rad = np.clip(sw_rad, 0, 30) - - # Build daily DataFrame - daily_data = { - "date": dates, + daily_truth = pd.DataFrame({ "temperature_2m_max": tmax, "temperature_2m_min": tmin, "temperature_2m_mean": tmean, "precipitation_sum": rain_sum, "precipitation_probability_max": precip_prob, "rain_sum": rain_sum, - "wind_speed_10m_max": wind_speed_max, - "wind_gusts_10m_max": wind_gusts_max, + "wind_speed_10m_max": wind_max, + "wind_gusts_10m_max": gusts_max, "wind_direction_10m_dominant": wind_dir, "shortwave_radiation_sum": sw_rad, "et0_fao_evapotranspiration": sw_rad * 0.4, "weather_code": np.where(rain_sum > 0, np.where(rain_sum > 10, 63, 61), 0), - } - daily = pd.DataFrame(daily_data).set_index("date") - daily.index = pd.to_datetime(daily.index) + }, index=pd.to_datetime(dates)) - # Generate hourly data with diurnal cycles + # Hourly truth hours_per_day = 24 total_hours = n_samples * hours_per_day - hour_timestamps = [start_date + timedelta(hours=i) for i in range(total_hours)] hour_of_day = np.tile(np.arange(24), n_samples) - - # Diurnal temperature: sinusoid between tmin and tmax, peaking at 14:00 - day_indices = np.repeat(np.arange(n_samples), hours_per_day) t_range = np.repeat(tmax - tmin, hours_per_day) t_phase = 2 * np.pi * (hour_of_day - 14) / 24 - t_hourly = np.repeat(tmin, hours_per_day) + t_range * (0.5 + 0.5 * np.cos(t_phase)) + rng.normal(0, 0.5, total_hours) - # RH: inverse of temperature cycle - rh_hourly = np.repeat(rh_mean, hours_per_day) - 5 * np.cos(t_phase) + rng.normal(0, 3, total_hours) - rh_hourly = np.clip(rh_hourly, 20, 100) + daily_rh = np.clip(78 + 12 * np.sin(2 * np.pi * (doy - 172) / 365.25) + rng_noise.normal(0, 6, n_samples), 45, 98) + dewpoint = tmean - ((100 - daily_rh) / 5.0) + rng_noise.normal(0, 0.5, n_samples) - hourly_data = { - "date": hour_timestamps, + t_hourly = np.repeat(tmin, hours_per_day) + t_range * (0.5 + 0.5 * np.cos(t_phase)) + rng_noise.normal(0, 0.5, total_hours) + rh_hourly = np.clip(np.repeat(daily_rh, hours_per_day) - 5 * np.cos(t_phase) + rng_noise.normal(0, 3, total_hours), 20, 100) + + hour_timestamps = [start_date + timedelta(hours=i) for i in range(total_hours)] + hourly_truth = pd.DataFrame({ "temperature_2m": t_hourly, "relative_humidity_2m": rh_hourly, - "dew_point_2m": np.repeat(dewpoint, hours_per_day) + rng.normal(0, 0.5, total_hours), - "apparent_temperature": t_hourly + rng.exponential(2.0, total_hours), - "precipitation_probability": np.repeat(precip_prob, hours_per_day) / 24 + rng.normal(0, 1, total_hours), - "precipitation": np.repeat(rain_sum, hours_per_day) / 24 * rng.uniform(0.5, 1.5, total_hours), + "dew_point_2m": np.repeat(dewpoint, hours_per_day) + rng_noise.normal(0, 0.5, total_hours), + "apparent_temperature": t_hourly + rng_noise.exponential(2.0, total_hours), + "precipitation_probability": np.clip(np.repeat(precip_prob, hours_per_day) / 24 + rng_noise.normal(0, 1, total_hours), 0, 100), + "precipitation": np.repeat(rain_sum, hours_per_day) / 24 * rng_noise.uniform(0.5, 1.5, total_hours), "rain": np.repeat(rain_sum, hours_per_day) / 24, - "cloud_cover": np.repeat(cloud_cover, hours_per_day) + rng.normal(0, 5, total_hours), - "cloud_cover_low": np.repeat(cloud_low, hours_per_day), - "cloud_cover_mid": np.repeat(cloud_mid, hours_per_day), - "cloud_cover_high": np.repeat(cloud_high, hours_per_day), - "wind_speed_10m": np.repeat(wind_speed_max, hours_per_day) * 0.5 * (0.5 + 0.5 * np.cos(t_phase)), - "wind_speed_100m": np.repeat(wind_speed_100m_max, hours_per_day) * 0.6, - "wind_gusts_10m": np.repeat(wind_gusts_max, hours_per_day) * (0.3 + 0.7 * rng.beta(2, 5, total_hours)), - "wind_direction_10m": np.repeat(wind_dir, hours_per_day) + rng.normal(0, 10, total_hours), - "surface_pressure": np.repeat(pressure, hours_per_day) + rng.normal(0, 0.5, total_hours), - "visibility": np.repeat(visibility, hours_per_day) + rng.normal(0, 500, total_hours), - } - hourly = pd.DataFrame(hourly_data).set_index("date") - hourly.index = pd.to_datetime(hourly.index) + "cloud_cover": np.clip(np.repeat(cloud, hours_per_day) + rng_noise.normal(0, 5, total_hours), 0, 100), + "cloud_cover_low": np.clip(np.repeat(cloud * 0.6, hours_per_day), 0, 100), + "cloud_cover_mid": np.clip(np.repeat(cloud * 0.3, hours_per_day), 0, 100), + "cloud_cover_high": np.clip(np.repeat(cloud * 0.2, hours_per_day), 0, 100), + "wind_speed_10m": np.repeat(wind_max, hours_per_day) * 0.5 * (0.5 + 0.5 * np.cos(t_phase)), + "wind_speed_100m": np.repeat(wind_max, hours_per_day) * 1.3 * 0.6, + "wind_gusts_10m": np.repeat(gusts_max, hours_per_day) * (0.3 + 0.7 * rng_noise.beta(2, 5, total_hours)), + "wind_direction_10m": np.repeat(wind_dir, hours_per_day) + rng_noise.normal(0, 10, total_hours), + "surface_pressure": np.repeat(pressure, hours_per_day) + rng_noise.normal(0, 0.5, total_hours), + "visibility": np.clip(15000 - np.repeat(rain_sum, hours_per_day) * 500 + rng_noise.normal(0, 2000, total_hours), 500, 25000), + }, index=pd.to_datetime(hour_timestamps)) - # Clip all values to realistic ranges - hourly["cloud_cover"] = np.clip(hourly["cloud_cover"], 0, 100) - hourly["cloud_cover_low"] = np.clip(hourly["cloud_cover_low"], 0, 100) - hourly["cloud_cover_mid"] = np.clip(hourly["cloud_cover_mid"], 0, 100) - hourly["cloud_cover_high"] = np.clip(hourly["cloud_cover_high"], 0, 100) - hourly["visibility"] = np.clip(hourly["visibility"], 100, 30000) - hourly["precipitation_probability"] = np.clip(hourly["precipitation_probability"], 0, 100) + # === Add NWP forecast errors === + daily_fc, hourly_fc = _add_nwp_forecast_error(daily_truth.copy(), hourly_truth.copy(), rng) - return daily, hourly + return daily_fc, hourly_fc, daily_truth, hourly_truth -def train_targets( - target_names: Optional[list] = None, - n_bootstrap: int = 5000, - test_split: float = 0.2, -): - """Train all or selected target models.""" +def train_targets(target_names=None, n_bootstrap=8000, test_split=0.2): + """Train models with proper NWP forecast → observation mapping.""" if target_names is None: target_names = list(TARGET_DEFINITIONS.keys()) - print(f"Training {len(target_names)} models...") - print(f"Bootstrap samples: {n_bootstrap} (test split: {test_split:.0%})") + print(f"Training {len(target_names)} models with realistic NWP errors") + print(f" Samples: {n_bootstrap} (test: {test_split:.0%})") print() - daily, hourly = bootstrap_training_data(n_bootstrap) + daily_fc, hourly_fc, daily_truth, hourly_truth = bootstrap_training_data(n_bootstrap) engine = FeatureEngine() - X = engine.transform(daily, hourly) - print(f"Features: {X.shape[1]} from {len(engine.FEATURE_GROUPS)} groups") - print(f" Thermal: {len(engine.FEATURE_GROUPS['thermal'])}") - print(f" Dynamic: {len(engine.FEATURE_GROUPS['dynamic'])}") - print(f" Moisture: {len(engine.FEATURE_GROUPS['moisture'])}") - print(f" Temporal: {len(engine.FEATURE_GROUPS['temporal'])}") - print(f" Interaction: {len(engine.FEATURE_GROUPS['interaction'])}") - print() + # Features from FORECAST (noisy NWP output) + X = engine.transform(daily_fc, hourly_fc) + print(f"Features: {X.shape[1]} from NWP forecast output") - # Train/test split (temporal order, no shuffle) - split_idx = int(len(daily) * (1 - test_split)) + split_idx = int(len(daily_fc) * (1 - test_split)) X_train, X_test = X[:split_idx], X[split_idx:] - daily_train, daily_test = daily.iloc[:split_idx], daily.iloc[split_idx:] + daily_truth_train, daily_truth_test = daily_truth.iloc[:split_idx], daily_truth.iloc[split_idx:] results = {} - # Feature augmentation: add Gaussian noise to prevent overfitting on synthetic data - X_train_noisy = X_train + np.random.RandomState(42).normal(0, 0.1, X_train.shape).astype(np.float32) - X_test_noisy = X_test + np.random.RandomState(43).normal(0, 0.05, X_test.shape).astype(np.float32) - for target_name in target_names: - print(f"{'='*60}") - print(f"Training: {target_name}") - print(f" {TARGET_DEFINITIONS[target_name]['description']}") + tdef = TARGET_DEFINITIONS[target_name] + print(f"\n{'='*60}") + print(f" {target_name} — {tdef['description']}") print(f"{'='*60}") - tdef = TARGET_DEFINITIONS[target_name] - y_train = WeatherModel.build_target( - daily_train, tdef["variable"], tdef["threshold"], tdef["op"] - ) - y_test = WeatherModel.build_target( - daily_test, tdef["variable"], tdef["threshold"], tdef["op"] - ) + # Targets from TRUTH (actual observation) + y_train = WeatherModel.build_target(daily_truth_train, tdef["variable"], tdef["threshold"], tdef["op"]) + y_test = WeatherModel.build_target(daily_truth_test, tdef["variable"], tdef["threshold"], tdef["op"]) p_yes = y_train.mean() * 100 print(f" Class balance: {p_yes:.1f}% YES / {100-p_yes:.1f}% NO") - model = WeatherModel(target_name) - model.train(X_train_noisy, y_train, X_test_noisy, y_test) + model = WeatherModel(target_name, mode="lr") + model.train(X_train, y_train, X_test, y_test) metrics = model.evaluate(X_test, y_test) model.save() results[target_name] = metrics - print(f" Brier score: {metrics['brier_score']:.4f}") - print(f" ROC AUC: {metrics['roc_auc']:.3f}") - print(f" Predicted mean: {metrics['p_yes_predicted']:.1f}% (actual: {metrics['p_yes_actual']:.1f}%)") - print(f" Top 10 features:") - for feat, imp in list(model.top_features(10).items()): + print(f" Pre-calibration Brier: — ") + print(f" Post-calibration:") + print(f" Brier: {metrics['brier_score']:.4f} AUC: {metrics['roc_auc']:.3f}") + print(f" Predicted mean: {metrics['p_yes_predicted']:.1f}% Actual: {metrics['p_yes_actual']:.1f}%") + print(f" Calibration: {metrics['calibration_method']}") + print(f" Top 8 features:") + for feat, imp in list(model.top_features(8).items()): print(f" {feat:30s} {imp:>10.1f}") - print() - # Summary print(f"\n{'='*60}") print("TRAINING SUMMARY") print(f"{'='*60}") - print(f"{'Target':<25s} {'Brier':>8s} {'ROC AUC':>8s} {'Cal Err %':>10s} {'Samples':>8s}") - print("-" * 62) + print(f"{'Target':<25s} {'Brier':>8s} {'AUC':>8s} {'Cal Err%':>9s} {'Cal':>10s}") + print("-" * 65) for name, m in results.items(): cal_err = abs(m["p_yes_predicted"] - m["p_yes_actual"]) - print(f"{name:<25s} {m['brier_score']:>8.4f} {m['roc_auc']:>8.3f} {cal_err:>10.1f} {m['n_samples']:>8d}") + print(f"{name:<25s} {m['brier_score']:>8.4f} {m['roc_auc']:>8.3f} {cal_err:>9.1f} {m['calibration_method']:>10s}") print(f"\nModels saved to: {MODEL_DIR}") return results @@ -287,30 +304,21 @@ def train_targets( def main(): parser = argparse.ArgumentParser(description="Train HK weather prediction models") - parser.add_argument("--target", type=str, default=None, help="Train single target (e.g., temp_gt_30c_24h)") - parser.add_argument("--bootstrap", action="store_true", default=True, help="Use bootstrap training data") - parser.add_argument("--samples", type=int, default=5000, help="Bootstrap sample count") - parser.add_argument("--all", action="store_true", default=False, help="Train all targets") + parser.add_argument("--target", type=str, default=None) + parser.add_argument("--samples", type=int, default=8000) + parser.add_argument("--all", action="store_true", default=False) args = parser.parse_args() - if args.target: - targets = [args.target] - elif args.all: - targets = list(TARGET_DEFINITIONS.keys()) - else: - targets = list(TARGET_DEFINITIONS.keys()) - + targets = [args.target] if args.target else list(TARGET_DEFINITIONS.keys()) results = train_targets(targets, n_bootstrap=args.samples) - # Save summary MODEL_DIR.mkdir(parents=True, exist_ok=True) - summary_path = MODEL_DIR / "training_summary.json" - with open(summary_path, "w") as f: + with open(MODEL_DIR / "training_summary.json", "w") as f: json.dump({ "training_date": datetime.now().isoformat(), - "n_bootstrap_samples": args.samples, + "n_samples": args.samples, "results": results, - }, f, indent=2) + }, f, indent=2, default=str) if __name__ == "__main__":