Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: scipy.stats.theilslopes returns invalid data when input data is np.uint8 #19678

Closed
korypostma opened this issue Dec 11, 2023 · 0 comments · Fixed by #19679
Closed

BUG: scipy.stats.theilslopes returns invalid data when input data is np.uint8 #19678

korypostma opened this issue Dec 11, 2023 · 0 comments · Fixed by #19679
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.stats
Milestone

Comments

@korypostma
Copy link

Describe your issue.

``

Expected result:
linregress = LinregressResult(slope=-0.08117408906882585, intercept=212.77543859649109, rvalue=-0.14521616256311687, pvalue=0.3777456213641338, stderr=0.09092294431311998, intercept_stderr=182.21245386044583) theilslopes = (0.0, 52.0, -0.2, 0.1388888888888889) siegelslopes= (-0.07448275862068965, 201.11448275862068) linregress = LinregressResult(slope=0.7799595141700404, intercept=-1526.3721997300943, rvalue=0.648915022489458, pvalue=7.885128179731413e-06, stderr=0.1503448757858782, intercept_stderr=301.2958824556652) theilslopes = (0.6923076923076923, -1350.3846153846152, 0.41379310344827586, 1.2307692307692308) siegelslopes= (0.980072463768116, -1927.0851449275362) linregress = LinregressResult(slope=-1.130971659919028, intercept=2291.15951417004, rvalue=-0.7989778176023127, pvalue=1.0820717823130094e-09, stderr=0.13994295315037716, intercept_stderr=280.4501007599699) theilslopes = (-1.0434782608695652, 2123.1304347826085, -1.3333333333333333, -0.5) siegelslopes= (-1.0947293447293447, 2218.258547008547) linregress = LinregressResult(slope=0.02186234817813765, intercept=7.469905533063432, rvalue=0.07626971562528158, pvalue=0.6444538083463693, stderr=0.04698691899214995, intercept_stderr=94.16327059776272) theilslopes = (0.0, 52.0, 0.0, 0.0) siegelslopes= (0.0, 52.0)

Current result:
linregress = LinregressResult(slope=-0.08117408906882585, intercept=212.77543859649109, rvalue=-0.14521616256311687, pvalue=0.3777456213641338, stderr=0.09092294431311998, intercept_stderr=182.21245386044583) theilslopes = (5.5, -10970.0, 1.0909090909090908, 11.272727272727273) siegelslopes= (-0.027777777777777776, 107.36111111111111) linregress = LinregressResult(slope=0.7799595141700404, intercept=-1526.3721997300943, rvalue=0.648915022489458, pvalue=7.885128179731413e-06, stderr=0.1503448757858782, intercept_stderr=301.2958824556652) theilslopes = (1.6666666666666667, -3303.0, 1.2857142857142858, 2.1333333333333333) siegelslopes= (0.04999999999999999, -62.39999999999998) linregress = LinregressResult(slope=-1.130971659919028, intercept=2291.15951417004, rvalue=-0.7989778176023127, pvalue=1.0820717823130094e-09, stderr=0.13994295315037716, intercept_stderr=280.4501007599699) theilslopes = (11.333333333333334, -22680.0, 8.73076923076923, 14.125) siegelslopes= (-0.1539473684210526, 339.12499999999994) linregress = LinregressResult(slope=0.02186234817813765, intercept=7.469905533063432, rvalue=0.07626971562528158, pvalue=0.6444538083463693, stderr=0.04698691899214995, intercept_stderr=94.16327059776272) theilslopes = (0.3333333333333333, -616.0, 0.07142857142857142, 0.625) siegelslopes= (0.0, 52.0)

Reproducing Code Example

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def plot_theil(xin, yin, png_path='out_theilslopes.png'):
    xin = xin.astype(np.int16)
    yin = yin.astype(np.uint8)
    #yin = yin.astype(float)
    lsq_res = stats.linregress(xin, yin)
    res = stats.theilslopes(yin, xin)
    sie_res = stats.siegelslopes(yin, xin)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(xin, yin, 'b.')
    ax.plot(xin, res[1] + res[0] * xin, 'r-')
    ax.plot(xin, sie_res[1] + sie_res[0] * xin, 'b-')
    ax.plot(xin, lsq_res[1] + lsq_res[0] * xin, 'g-')

    print('linregress = ', lsq_res)
    print('theilslopes = ', res)
    print('siegelslopes= ', sie_res)

    #plt.show()
    plt.savefig(png_path, bbox_inches='tight')

xin = np.array(range(1985,2023+1))
yin = np.array([46,46,51,54,49,50,52,53,52,54,56,53,56,57,57,55,53,52,51,49,55,54,38,44,37,40,31,38,45,50,53,58,57,53,56,55,51,48,45])
plot_theil(xin, yin, 'out_theilslopes.png')

yin = np.array([34,41,38,37,35,41,25,38,17,29,14,13,14,13,21,15,27,25,36,32,42,33,38,33,31,29,43,53,53,53,53,53,53,53,53,53,53,53,53])
plot_theil(xin, yin, 'out_theilslopes_1.png')
yin = np.array([55,50,40,50,49,35,37,42,34,41,32,13,11,13,17,38,36,33,33,32,38,37,37,40,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8])
plot_theil(xin, yin, 'out_theilslopes_2.png')
yin = np.array([52,52,52,52,52,52,52,51,54,54,44,40,55,52,52,52,52,52,50,51,46,46,47,61,50,49,54,52,52,52,52,52,52,52,52,52,52,52,52])
plot_theil(xin, yin, 'out_theilslopes_3.png')

Error message

No error

SciPy/NumPy/Python version and system information

import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.8.1 1.22.4 sys.version_info(major=3, minor=10, micro=4, releaselevel='final', serial=0)
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
  NOT AVAILABLE
openblas_clapack_info:
  NOT AVAILABLE
flame_info:
  NOT AVAILABLE
accelerate_info:
  NOT AVAILABLE
atlas_3_10_threads_info:
  NOT AVAILABLE
atlas_3_10_info:
  NOT AVAILABLE
atlas_threads_info:
  NOT AVAILABLE
atlas_info:
  NOT AVAILABLE
lapack_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas']
    library_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\lib']
    language = f77
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
atlas_3_10_blas_threads_info:
  NOT AVAILABLE
atlas_3_10_blas_info:
  NOT AVAILABLE
atlas_blas_threads_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
blas_info:
    libraries = ['cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\lib']
    include_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\include']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    libraries = ['cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\lib']
    include_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\include']
    language = f77
lapack_opt_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas', 'cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\lib']
    language = f77
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/kpostma/Anaconda3/envs/rcmap\\Library\\include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2,AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL
    not found =
@korypostma korypostma added the defect A clear bug or issue that prevents SciPy from being installed or used as expected label Dec 11, 2023
@j-bowhay j-bowhay added this to the 1.13.0 milestone Dec 12, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.stats
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants