weighted decay code

Here's some code to generate those basis functions from yesterday. As you might imagine, you need to increase the decay rate as you ask for more basis functions. I didn't least squares it, just applied the usual trick with inner product to get the true parameter values for the basis functions. If you were anything other than the most truly lazy of gods creations you would derive an upper triangular form for c iteratively. You might also have the good sense to perform any of a million things that would make this more stable, such as using x^k/k! as the elements we multiply the coefficients against for instance. Or maybe exploiting the alpha parameters in the gamma function so this scales right.

My own failings aside you would be wrong to criticize this on the condition number of G in comparison to any polynomial based method. We're integrating over infinite distance, not the finite distances a polynomial asks of us, it's a different problem. If you call the fit bad in the final plot I'll ask you to have some grace as this is recovery of the coefficients for the true parameters in the infinite series, not the line of best fit.

import numpy as np
import scipy.special as sp

#Generate matrix of inner products of all functions up to order n-1

#of form x^(k)e^-theta*x over 0,inf
n = 10
x = np.arange(n)
X = x[:,np.newaxis] + x[np.newaxis,:] + 1
theta = 1
Gl = sp.gammaln(X)-(X*np.log(2*theta))
G = np.exp(Gl)

#Use matrix of inner products to determine coefficients for

#orthonormal basis
u,s,v = np.linalg.svd(G)
c = (u/np.sqrt(s))

#Evaluate othonormal basis functions
dt = 0.001
t = np.arange(dt,30,dt)
V = np.vander(t,N=n,increasing=True)
w = np.exp(-t*theta)
P = V*w[:,np.newaxis]
y = P @ (c)

#Plot basis functions
import matplotlib.pyplot as plt
fig, axs = plt.subplots(2, 2)
axs[0, 0].plot(t,y)
axs[0, 0].set_title('Basis functions')
print('integral', np.sum(y**2,axis=0)*dt)

#Demonstrate fitting property
f = np.sin(t)
b = y.T @ f * dt

#Plot fitted data
axs[1, 0].plot(t[:int(len(t)/3)],(y@b)[:int(len(t)/3)])
axs[1, 0].plot(t[:int(len(t)/3)],f[:int(len(t)/3)])
axs[1, 0].set_title('Local fit')
axs[1, 1].plot(t,(y@b))
axs[1, 1].plot(t,f)

axs[1,1].set_title('\n Full forecast')