<div dir="ltr"><div>You sent your NCL script and test data offline. <br><br></div><div>The time series that it fails on has 29737 values... ~84% (24878) are 0.0<br></div><div>The is a really unusual distribution.<br></div><div><br></div><div>NCL's Mann-Kendall code is a fortran subroutine which does two independent tasks:<br></div><div>[1] Compute the Mann-Kendall (MK) trend statistic; [2] Compute the Theil-Sen (TS) robust estimate of linear trend. <br><br>[1] I am sure the MK statistic and probability estimate were computed. <br>[2] The TS estimate requires that the 'slopes' calulated within the do loop be sorted into ascending order. The code fails in the fcode used for sorting. Specifically, at the following location of DPSORTU<br><br>[SNIP]</div><div><br> 65 I = I + 1
<br>
IF (I.EQ.J) GO TO 55
<br>
1341 ----> T = A(I+1)
<br>
IF (A(I).LE.T) GO TO 65
<br>
K = I <br></div><div>[SNIP]<br><br></div><div>The subroutine used is a 'famous' one translated from ALGOL (original author: ?Wilkinson?) to fortran (Singleton).
It is part of the ACM codes and is available via netlib.
<br>
<br>
C Sorts A into ascending order, from A(1) to A(N).
<br>
C Reference: Richard C. Singleton, Algorithm 347, SORT. <=======
<br>
C Comm. ACM 3, 321 (March 1969).
<br>
C Algorithm is Copyright 1969 Association of Computing Machinery,
<br>
<br><br></div><div>I speculate that the preponderance of 0.0 values has 'somehow' causes this code to fail. Why it required > 25,000 values to fail .. I do not know. Again, I speculate, if there had been more 0.0 in the 1st 25,000 it may have failed with fewer points.<br><br></div><div>There is no chance that I am going to change this sorting code. It requires more knowledge about sorting than I have. NOTE: the code used by NCL's 'dim_pqsort' is based on this code so I doubt this would be an acceptable alternative. That code is also downloadable from netlib.<br><br></div>----<br><div>OK .... MK is 'distribution free' BUT I think some knowledge of how MK works should be used. Basically, it checks to see if successive values go up (+1), down (-1) or are tied (0). The *amount *of the change is not used for MK. With series like 0,0,0,#,0,0,0,0,0,#,#,0,0,0,0,0,0,#,0,0,#,0,0,... etc. It will often go from 0.0 to some positive number (+1) and from some positive number to 0.0 (-1). I would speculate that the +1 and -1 counts almost nullify each other. All cases where the series has 0,0 would be a tie. This will be *by far* the largest number. The algorithm uses a correction for ties but I am not sure if was developed for cases where the number of ties is >> greater that the sum of the +1 and -1. <br><br></div><div>Also, some knowledge of hoe Theil-Sen estimates are made. In the process of calculating the sums and ties 'slopea' are calculated. The median of the slopes is the TS estimate. Since the majority of slopes are 0, the TS median value will be 0.0.<br><br></div><div>Just using regline ...<br><br> rc = regline(year_fraction,precip_ave) ; slope <br> print(rc)<br><br></div><div><br>Variable: rc <==== ALL <br>Type: float<br>Total Size: 4 bytes<br> 1 values<br>Number of Dimensions: 1<br>Dimensions and sizes: [1]<br>Coordinates: <br>Number Of Attributes: 8<br> units : mmhr^-1/3hrs<br> _FillValue : -999.99<br> yintercept : 2.898926<br> yave : 0.112996<br> xave : 2006.162<br> nptxy : 29737<br> rstd : 0.000660348<br> tval : -2.102962<br>(0) -0.001388687<br><br><br><br> ii = ind(precip_ave.gt.0.0) ; observations of when it rains<br> RC = regline(year_fraction(ii),precip_ave(ii))<br> print(RC)<br><br>Variable: RC <=== Conditioned on rain only times.<br>Type: float<br>Total Size: 4 bytes<br> 1 values<br>Number of Dimensions: 1<br>Dimensions and sizes: [1]<br>Coordinates: <br>Number Of Attributes: 7<br> _FillValue : -999.99<br> yintercept : -6.783739<br> yave : 0.6915334<br> xave : 2005.787<br> nptxy : 4859<br> rstd : 0.003920668<br> tval : 0.9505657<br>(0) 0.003726852<br><br></div><div><br><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Apr 5, 2015 at 9:04 AM, Mary Haley <span dir="ltr"><<a href="mailto:haley@ucar.edu" target="_blank">haley@ucar.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-size:small">Dave,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Could you write the precip_ave variable to a NetCDF file and send us the file, so we can try trend_manken on it here?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">You can do something like:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">fout = addfile("<a href="http://test_trend_manken.nc" target="_blank">test_trend_manken.nc</a>","c")</div><div class="gmail_default" style="font-size:small">fout->pcp = precip_ave</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">and then if the file is too large (say over 25 megabytes), then you can use our ftp to upload it. See </div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default"><a href="http://www.ncl.ucar.edu/report_bug.shtml#HowToFTP" target="_blank">http://www.ncl.ucar.edu/report_bug.shtml#HowToFTP</a><br></div><div class="gmail_default"><br></div><div class="gmail_default" style="font-size:small">Thanks,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">--Mary</div><div class="gmail_default" style="font-size:small"><br></div></div><div class="gmail_extra"><br><div class="gmail_quote"><span class="">On Sat, Apr 4, 2015 at 12:56 PM, David Adams <span dir="ltr"><<a href="mailto:dave.k.adams@gmail.com" target="_blank">dave.k.adams@gmail.com</a>></span> wrote:<br></span><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><div dir="ltr"><div><div><div><div><div>Hi all, <br></div>I am using the trend_manken<br></div>function to look at the significance of precipitation intensity trends (using TRMM data) for 15 years (every 3 hours data). When I use 10,000 data points it takes about 3 seconds to run, 20,000 data points about 10 seconds, 25,000 about 15 seconds. When I use the entire time series 29,737 data points, it never finishes. I have waited more than 20 minutes and it keeps spinning. I switched from double to floating point, but that makes no difference.<br><br>;---------------------------------<br>; Mann Kendall test<br>;---------------------------------<br>pa = trend_manken(precip_ave, False, 0)<br>print(pa)<br><br><br><br></div>Any ideas?<br><br></div>saludos,<br></div>Dave<br></div>
<br></div></div><span class="">_______________________________________________<br>
ncl-talk mailing list<br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a><br>
<br></span></blockquote></div><br></div>
<br>_______________________________________________<br>
ncl-talk mailing list<br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>