I have two data arrays, xdata of shape (40,) and ydata of shape (40, 721, 1440) (time, lat, lon). My final goal is to compute the regression slope between these two datasets and obtain the errors from the slope to show the error distribution, where x is average sea surface temperatures for 40 years and y is an atmospheric variable. I have done this using one method where I calculated the covariance and then took the square root of this array, which if I understand correct, leaves me with the standard error. To verify this approach, I found the scipy.stat.linregress function and wanted to use this as it returns the standard error while calculating the slope. Although, I am running into errors when coding this.
To get my arrays the same size:
Y = ydata.stack(allpoints = ['lat','lon'])
X = xdata.values[:, None] * np.ones(Y.shape)
Here, I have stacked my ydata into a 2D array rather than a 3D array, leaving with me with the dimensions of (40, 1038240). Then I create a temporary array to of the xdata, with an extra dimensions filled with ones to get the two arrays of the same shape. After this I pass it through the function:
test = stats.linregress(X,Y)
And I am left with a value error saying:
ValueError: too many values to unpack (expected 4)
UPDATE:
I was able to get the package I was interested in using working:
from scipy import stats
ny = 721
nx = 1440
n = nx * ny
cape_2d = cape_ds.values.reshape(40,n)
reg_results = np.empty((5,n))
for i in range(n):
reg_results[:,i] = stats.linregress(b,cape_2d[:,i])
slope, intercept, r_val, p_val, std_err = reg_results.reshape((5,ny,nx))
err_shape = std_err.reshape(721*1440)
plt.plot(err_shape)
At this point, I think I may have the correct answer, but Im worried that my standard errors are too high...Im not sure what to for it to look like, I was just hoping to get a normal distribution.