Fitting a power law for extrapolation

Question

I am trying to figure out how to interpolate a function, and I want to be able to extrapolate a small way beyond the interpolation range. There is some background theory which leads me to expect that the behavior for large values of the dependent variable x will asymptote to −A·x⁻ⁿ behaviour (maybe with some offsets). I want to find the best fit power law, but only for large values of x. Is there a good way of implementing this in Python?

I tried to do least squares for the log-linearized version of my power law, but it gave me an inappropriate fit - which I suspect is because of a coding error rather than a math mistake. I was thinking of doing this by introducing an offset into the power law, so that the log equation would be ln(−y) = ln(A) − n ln(x−x₀). What I have so far is:

X_asymptote=[]
for i in aX_vals1:
    if i>20*lab_inc(T-1):
        X_asymptote.append(i)
VT1_optimized_asymptote=[]
for i in X_asymptote:
    a=aX_vals1.index(i)
    VT1_optimized_asymptote.append(VT1_optimized_val[a])
#Converting to logs. Note the minus sign! You'll need to put this back in in the unravelling.

X_asymptote_logs=[]
for i in X_asymptote:
    X_asymptote_logs.append(np.log(i))
X_asymptote_logs_array=np.asarray(X_asymptote_logs)

VT1_optimized_asymptote_logs=[]
for i in VT1_optimized_asymptote:
    VT1_optimized_asymptote_logs.append(np.log(-i))
    #VT1_optimized_asymptote_logs.append(-1*np.log(i))
VT1_optimized_asymptote_logs_array=np.asarray(VT1_optimized_asymptote)


#Linear fit
fitfunc=lambda p,x:p[0]+p[1]*x
errfunc=lambda p,x,y:(y-fitfunc(p,x)) 

out=scipy.optimize.leastsq(errfunc,[-4,5],args=(X_asymptote_logs_array,VT1_optimized_asymptote_logs_array),full_output=1)
pfinal=out[0]
covar=out[1]

index=pfinal[1]
amp=np.exp(pfinal[0])

xnew=np.linspace(X_asymptote_logs[0],aX_vals1_array[len(aX_vals1_array)-1],100)
ynew=VT1_optimized(xnew)
plt.title("Optimized value function, linear extrapolation")
plt.plot(aX_vals1_array, VT1_optimized_val_array, 'o', xnew, ynew, '-')
plt.ylabel('Value function with optimized asset allocation')
plt.xlabel('ratio X/Y')
plt.show()

Could anyone give me a hand with this?

Any help is much appreciated!

Update: I have a pragmatic solution to the problem. I just fitted the power law for the whole data range. It seems to have worked remarkably well, and the code now does want I want, but not in the way I was hoping to understand with this comment. The code I now have is:

#Linear fit
fitfunc=lambda p,x:p[0]+p[1]*(x)
errfunc=lambda p,x,y:(y-fitfunc(p,x)) 

out=scipy.optimize.leastsq(errfunc,[-4,5],args=(np.asarray(X_asymptote_logs),np.asarray(VT1_optimized_asymptote_logs)),full_output=1)
pfinal=out[0]
covar=out[1]

index=pfinal[1]
amp=np.exp(pfinal[0])
print(index,amp)

def VPL(X):
    return -(amp*(X)**(index))

def VT1_optimized1(X):
    if aX_vals1_array[len(aX_vals1_array)-1]>=X:
        VT1_op=scipy.interpolate.UnivariateSpline(aX_vals1_array,VT1_optimized_val_array,k=3,s=0,check_finite=True)
        return VT1_op(X)
    if X>aX_vals1_array[len(aX_vals1_array)-1]:
        VT1_asymptote=scipy.interpolate.UnivariateSpline(aX_vals1_array,VT1_optimized_val_array,k=5,s=0,check_finite=True)
        return VPL(X)


xnew=np.linspace(aX_vals1_array[0],aX_vals1_array[len(aX_vals1_array)-1]*2,100)

y0=[]
for i in xnew:
    y0.append(VT1_optimized1(i))
ynew=np.asarray(y0)
plt.title("Optimized value function, linear extrapolation")
plt.plot(np.asarray(aX_vals1), np.asarray(VT1_optimized_val), 'o', xnew, ynew, '-')
plt.ylabel('Value function with optimized asset allocation')
plt.xlabel('X')
plt.show()

#This seems to work almost too well! Let's check the asset allocation.




# Plotting
# xnew=np.linspace(aX_vals1_array[0],aX_vals1_array[len(aX_vals1_array)-1],100)
# ynew=VT1_optimized(xnew)
# plt.title("Optimized value function, linear extrapolation")
# plt.plot(aX_vals1_array, VT1_optimized_val_array, 'o', xnew, ynew, '-')
# plt.ylabel('Value function with optimized asset allocation')
# plt.xlabel('ratio X/Y')
# plt.show()


# xnew=np.linspace(aX_vals1_array[0],(aX_vals1_array[len(aX_vals1_array)-1]),100)
# y0=[]
# for i in xnew:
#     y0.append(VT1_optimized1(i))
# ynew=np.asarray(y0)
# plt.title("Optimized value function, lin. extrapolation")
# plt.plot(aX_vals1_array, VT1_optimized_val_array, 'o', xnew, ynew, '-')
# plt.ylabel('asset allocation')
# plt.xlabel('X')
# plt.show()
# exit()

This question is more about math than any particular language. — Axalix, Dec 29 '17 at 21:10
I'm not really concerned about the math. I could well be making a math error! But I'm hoping to be able to implement the partial fit in python to see just that! — лаетековский чайник, Dec 29 '17 at 21:13
Actually, after edits I think I see your point. I still lean to thinking it's a coding question, because I just want to see what a reasonable way of doing the partial fit is. I guess the mathematical question is how you sew these two things together, whereas I think I'm just asking how you implement a given sewing. — лаетековский чайник, Dec 29 '17 at 21:20
I don't understand why you don't use something like [this](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.interpolate.InterpolatedUnivariateSpline.html) for example. — Georgy, Dec 29 '17 at 21:30
Have you tried using better initial guesses for the fit parameters? — mkrieger1, Dec 29 '17 at 21:36
Hi Georgy! I think just an interpolated spline won't work, because an inverse power-law has an infinite Taylor Series, so the extrapolation would break down very quickly, as indeed happens for my code. The -4 parameter is the theoretical value, so I think that's a good estimate. The offset I admit to just guessing! — лаетековский чайник, Dec 29 '17 at 22:27
Yeah, some data would be nice. Tip, use `zip` to avoid half of the loops and the `_.index()` part that does not really increase readability. — mikuszefski, Jan 12 '18 at 07:24

Fitting a power law for extrapolation

0 Answers0