From 7a040fd599eda42ee451305c6caa6412dbb4e6cb Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 12 Jul 2020 16:23:09 -1000 Subject: [PATCH] Confidence band on regression missed multipler of sigma (#3650) * Confidence band on regression missed multipler of sigma The confidence band is proportional to the rms misfit (sum of squared misfits w.r.t. the regession line) divided by n-2. then sqrt. This term was in fact missing. * Update PS origs and set 99% confidence interval --- doc/examples/ex43/ex43.ps | Bin 44533 -> 44397 bytes doc/examples/ex43/ex43.sh | 2 +- doc/rst/source/gallery/ex43.rst | 2 +- src/gmtregress.c | 4 ++-- test/gmtregress/draper.ps | Bin 34471 -> 34467 bytes 5 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/examples/ex43/ex43.ps b/doc/examples/ex43/ex43.ps index 88a24002ae0eb249ed540a14a637fa7b25b83bba..80455cc2a6f9f229ee296e467ab88bde8bf41897 100644 GIT binary patch delta 663 zcmYjOyGlbr5GC;u4O$5{f^1_EOoMvuDnnGjsQ6l>Pq9-qn_&NDtxY*h0;p6ljs^B7;!hkSY61F6So6B^{Icg;n8gdr&F93hBNA^IMHprjBuFC}W(2YUoLmNcLYNy~=MtWNLAH-)g} z?3AMd@*^ejv_Fuq&J&uRvqdoUjJRaykdX4(Ina#VT$owkdPmc+rGSl5HgjUkNx979 z7zHSl2A>>lXSvoqfkw|^&E?vWu3*oYRW$~}*0C<(z*u;;pb@!A;gQ{H;Ttqk&?T9! z8NGuBv8uUmF~c$t#R9ERCtgq|0zt%?l_xB$F7y!29o9K`H^7qT5~J z&r8aClRO)Qa|YmSBKQw;W5d3M9k0|+pV3C|)Q=uNm_C2Cw%-0|E?4v;58luDV~*#$ zblVxu?zXo-Y+NYt model.txt - gmt regress -Ey -Nw -i0:1+l $file -Fxmc -T-2/6/0.1 > rls_line.txt + gmt regress -Ey -Nw -i0:1+l $file -Fxmc -T-2/6/0.1 -C99 > rls_line.txt gmt regress -Ey -N2 -i0:1+l $file -Fxm -T-2/6/2+n > ls_line.txt grep -v '^>' model.txt > A.txt grep -v '^#' $file > B.txt diff --git a/doc/rst/source/gallery/ex43.rst b/doc/rst/source/gallery/ex43.rst index 99ddf2a2dcb..9ee6c4cba1c 100644 --- a/doc/rst/source/gallery/ex43.rst +++ b/doc/rst/source/gallery/ex43.rst @@ -8,7 +8,7 @@ robust regression line using *reweighted least squares* and from this fit we are able to identify outliers with respect to the main trend. The result shows dinosaurs were large and dumb, humans and some monkeys pretty smart, and the rest of the -mammals doing alright. +mammals doing alright. Band shows 99% confidence internal on regression. .. literalinclude:: /_verbatim/ex43.txt :language: bash diff --git a/src/gmtregress.c b/src/gmtregress.c index db812d5c5fb..4cdd4dd5a68 100644 --- a/src/gmtregress.c +++ b/src/gmtregress.c @@ -1335,8 +1335,8 @@ EXTERN_MSC int GMT_gmtregress (void *V_API, int mode, void *args) { case 'r': /* Residual */ out[col] = S->data[GMT_Y][row] - gmtregress_model (x[row], par); break; - case 'c': /* Model confidence limit (add x and y uncertainties in quadrature since uncorrelated) */ - out[col] = t_scale * hypot (par[GMTREGRESS_SIGIC], par[GMTREGRESS_SIGSL] * fabs (x[row] - par[GMTREGRESS_XMEAN])); + case 'c': /* Model confidence limit (add slope and intercept uncertainties in quadrature since uncorrelated) */ + out[col] = t_scale * sqrt (par[GMTREGRESS_MISFT]) * hypot (par[GMTREGRESS_SIGIC], par[GMTREGRESS_SIGSL] * fabs (x[row] - par[GMTREGRESS_XMEAN])); break; case 'z': /* Standardized residuals (z-scores) */ out[col] = z_score[row]; diff --git a/test/gmtregress/draper.ps b/test/gmtregress/draper.ps index 81e8b37fe51b4945121bfc7bff235a8015eb3e8a..b8d3917ad680907095c69b6c4ef4af288f24d828 100644 GIT binary patch delta 2975 zcmc&$OK%%h6jn`XsVh}2h^k3nm!?)#g>mNIvBzVVMs3ulNfgH+F$5au)UhYdG!u`p zCy7f`+(%m26bN@)!GcOGK;mJRpdwhn0wET#VME~$lnqoYpxkq=J)Wp(p(_%_J~{V1 zzwg|~dsjp6UJbo@yf3PXACgEiDUiTVI1 zQnpf1x9Y;S->sy`d}+4iE|ylVAKEcKUm`>E4xx~rTu#8R=(t)+;< zv=~V+64MRB41@2z(#gLW*{7(zf+CyH%@Mm z4RCk5gCCmg+Sn|Qo0A9lTlpt9hXIGYW0%QDa0KDRYW~P|b*>GTrjZH0l2?Moab;>J zzj=CatIrN8S3Rrd@t*Vh_=i)E25K?$t84iz|8=mITj~9LkGFw;Vi~^EQ~wul+4YYD zt#?t%GqVm)%sj?_nti|VYx#rFCLjifhre2_u7WnDnwEeq8Nt7DiDCux&dmG z8N^Zm1gb*?3KVIy*EL_6grEvEhJB$Sl?lLys=g(qP_;D0NVqN>g%Xnv$pjEf#|qIv z!&1ONoULFBRj~y0tsr4>42xO=`W_AUXZX*jhXR9;UiXbC1z4n&(Iy*Fky_*(N0aG( z9a=lNL`Wq`-AAg)A0WkWh9#>Bqz$6hspiWRD@KZ9K1InkgOUL7qT;VbCUS+!1S*j> z1V{2}U@>T5U@A#`GpizkB{2k2s3V(5thiM~3ONKo_Zr0ZyxhE&@EqN^^g!om@6yh? zE+)N2EpAcR6$(p2J`^?FR=9y!959ja<19l<3dLEKT5-3v2o(2WD*_H^5^Q65U|=x0 z3F|6}p(S80WCO7_#B|7wD~BDov&0fhkXnf_g2oUuNdbb5bSF`SE;Q)|HN*_uIk?Fv z+0aC_d8Gul%$;iRUV|(B({rou61dl*!8^#YidLyWlw)VIW8I51s>ki_Gh^9qSB+A; z2R^uaa`z+s-DT=ue-QZh^IhiJ2Q|KQWLR^`t(l8mv#@26Hp!<9Xz6pD}%~0@tajf6nHWU-W6qZEHh*d zA0Ii99Wb968O$EGOX6jsD8k7c7#z=L#-2YsopY)dT<70j7{#|w;bI%8*Dih@0Cnh+ zcCeLaq?L#LTf9BX?Rb65?}Vltesg&beElWAIxiorhl{Oe3fzwOcpd(Smt))e-Vet> an3hrZ{hyb8R?#k**-FlKtO{SeqW%L`fi~9w delta 3157 zcmcImUrd`-6b}{MsM9QNGsAzo88lnMzH;yV{+FRLXc-GFjDC#8Mv%e>l+{vOCIqrI zWlW!zxLhBXn8g?$=7Xlb8T3WtBfgrLCBAEn555}jIk$hXBQa?s3GL~*=l;&Qzu&p% z`#!xF`21erqv#o8yF3^Qd71UBA0eZ0kGyY&d6>=>nC`gF{SWDJQ?AF>3K5du@QX!1 zd;iCQ)85KWa%Db8h(nkY(RqZ~gi~f#Uh4R}=k?H&e`I)i>RNP!TpCJ7$Hx*;;zh4} zWO96HI2w;8JQABu40~e}iDa-}=Pwg~b76(VZ;?fRsq$!f@2T!XVyesK8_CWJxxO2y z);c%adt95LVQ%Ix6#T^tnpCjWvL{nL9ZYPc&Q*R+X@O@_Zxc!`iK`otdagHG=j8c0 z03WS;O_f2f^i-%q!S)-4_`{f|B4b)IT!Db-G*d?PC|D<{OU?4!!b_r3^-pNiRr zHsY79V@H<@zMpGDpx@SePBaT-$h+{)%`+v&ImiR*Fx!>l7!H_5!!%9l<%aY+Af9f) zI=Aaih{R=Amlor4*>Gw;CV{I?lQ*38^%7&S3->cthhn^rWr3xVx>h|@hoaESDS#=H zK?*zocK}OyWkt0NSQsVDumDQdS%M#sL4#PvAXK`giREfg7}A0k5}+bEWK~YcgccH@ zEL;UHOyCTCjABkeWOxn0yDDvFl$cd(ELAOv!Q>FG^_UH(RqbX2 z2P#0XH-QVG*conJOorX4A}!!RnRcT|YE4>_UFra=>IXT6YNiFaY93)S6<~u`N<$^5 zRpzQJ+fom1Hkw37l^`+G>xwWbtK5cO?6Oi-bliuc(}vKh&r~ZwUKPu*YY@UgebICg z7-KoM8mf+mQ2OeKi&7h~EjA$)jf&|F>tA!Eji+gPyQkDVKybouGeW~A<<-yOTO^iL zW8ZXPt{|J5A!A2^K%>S{R5+K$@MZL|5jO4+Rd1@Fg)KAF zKjrn)Fw5YF*xMpUV|{tn&U*nX_x3k;jcMaMfxMVr(3bO(x$?TXvc+}RCdq>Mx_qVD zv`2QdF8Q)YDTDdcw>#O&YvAJ@IDOZhy(;Mb9bJLKYj~@mAMXO_ce|g6`Ma-*viyH| p_q>?h(CXLiY^T@yPz^9cTuq&QxtTh7r*n%F=~aJWT@*i{{{Sqi2XO!Z