I'm working with an ensemble of precipitation data and I want to post-process. One technique that I want to try is EMOS (Ensemble Model Output Statistics). There is an R package called EnsembleMOS that does the job. However, since most of my other scripts are Python, I want to translate the ensembleMOScsg0 part of that package to Python.
The problem that I'm facing is that scipy.optimize.minimize is returning the initial guess. I've found many different questions with a similar issue, and the answers were almost always related to some kind of problem on the function that is minimized.
My function:
def crps(train,obs,pars):
a = pars[0]
c = pars[1]
d = pars[2]
shift = pars[3]
# print(a,c,d,pars[6])
media = a
crps = 0
for i in range(len(train)):
var = np.var(train.iloc[i])
for j in range(len(train.columns)-1):
media += pars[j + 4]*train.iloc[i,j]
dp = c + d*var
shape = (media**2)/(dp)
scale = dp/media
gamma = ss.gamma.cdf
beta = ss.beta.rvs
x = obs.iloc[i]
z = (x + shift**2)/scale
cc = (shift**2)/scale
crps += scale*z*(2*gamma(z,a=shape)-1) \
- scale*cc*(gamma(cc,a=shape))**2 \
+ media*(1+2*gamma(cc,a=shape)*gamma(cc,a=shape+1)-gamma(cc,a=shape))**2 \
- 2*gamma(z,a=shape+1) \
- media*(1-gamma(2*cc,a=shape)) \
* beta(0.5,shape) \
/np.pi
return crps
The objective of the minimization is to find the parameters of a linear regression equation and its variance, that's why I separated the first 4 pars (a the intercept, c and d from the variance and the shift part of the gamma distribution) from the other 12 (the b's, since I have 12 ensemble members, I need 12 b's).
The minimization part is:
pars = [0.0001,0.0001,0.0001,0.0001] + [0.0001]*12
bnds = ((.0001,20),) * 4
bnds = bnds + ((.0001,20),) * 12
result = minimize(partial(crps,train,obs),pars,method='L-BFGS-B',bounds=bnds)
The process ends with a success flag, but it always converges to the initial values. Now I'm thinking about using rpy2 and running this package directly with Python.
More information about this procedure can be found on Baran and Nemoda (2016).