0

I have difficulty performing scipy.optimize.minimize with Martin Eastwood's interpolation formula—

z=(x^w1/(x^w2+y^w3))*w4 * 17 (we get 16 instead of 17
while x[3], x[4], x[16], x[18] reside in the formula)

My data set (17/12/12 preml.ge)

x=np.array([33,43,28,26,28,30,26,24,15,21,23,28,19,18,19,22,15,19,18,15])
y=np.array([15,24,17,16,21,25,22,21,13,20,23,29,25,24,26,32,24,31,32,30])
z=np.array([36,42,29,24,27,29,23,27,24,23,22,20,25,16,17,15,18, 9,15,10])
data=np.array([x, y, z])

Ten years ago, Martin Eastwood (an enthusiast blogger) found:

w1=1.122777, w2=1.072388, w3=1.127248, w4=2.499973
where RMSE=3.657522858 for my problem.

What I want to know is which approach I could use to get these w–parameters, like those, for above dependent estimation.
I read these answers, but the method seems not easy to trace to me. I need your help.
Added. A further question, how can we estimate w–parameters intended for each set {x_i, y_i, z_i} instead of whole {x, y, z} as above?

1 Answers1

0

Using least-squares is better because the method sees variation of each sample individually instead of only the final sum.

import scipy.optimize
import numpy as np
import matplotlib.pyplot
x=np.array([33,43,28,26,28,30,26,24,15,21,23,28,19,18,19,22,15,19,18,15])
y=np.array([15,24,17,16,21,25,22,21,13,20,23,29,25,24,26,32,24,31,32,30])
z=np.array([36,42,29,24,27,29,23,27,24,23,22,20,25,16,17,15,18, 9,15,10])

pred = lambda w: (x**w[0]/(x**w[1]+y**w[2]))*w[3
w_given = 1.122777, 1.072388, 1.127248, 2.499973]
w,_ = scipy.optimize.leastsq(lambda w: (z - pred(w)), (1,1,1,1))
w_guided,_ = scipy.optimize.leastsq(lambda w: (z - pred(w)), w_given)

Let's visualize

plt.plot(z, pred(w), '.')
# 17 introduced here arbitrarily
plt.plot(z, pred(w_given)*17, '+')
plt.plot(z, pred(w_guided), '+')
plt.plot(np.sort(z), np.sort(z), '--')
plt.legend(['dumb guess', 'given w (scaled)', 'init with given w', 'target'])

Checking if the fit gave better results than the initial guess (sanity check)

(np.mean((z - pred(w))**2), 
 np.mean((z - pred(w_guided))**2), 
 np.mean((z - pred(w_given)*17)**2), 
 np.mean((z - pred(w_given)*16)**2))

(10.987132120174204,
 10.987132121290418,
 15.064715846376691,
 17.341093598858798)
Bob
  • 13,867
  • 1
  • 5
  • 27