pmslavin · June 28, 2022 21:48
diff --git a/cholesky.py b/cholesky.py
 """
   Second-degree polynomial fit of f(xseries) => yseries
   by Choleskification and Ly=b, L^x=y
 """
 import numpy as np
 import matplotlib.pyplot as plt

 """
 # Test series: y=x^2
 xseries = [0, 0.5, 1.0, 1.5, 2.0, 2.5]
 yseries = [0, 0.25, 1.0, 2.25, 4.0, 6.25]
 """
 xseries = (0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5)
 yseries = (11, 4, 4, 6, 9, 12, 17, 18, 14, 18, 20, 20)
 # Good fit to smooth polynomial-like func...
 #yseries = [np.sin(x/12.0*np.pi) for x in range(12)]
 # Bad fit to periodicity...
 #yseries = [np.sin(x/6.0*np.pi) for x in range(12)]


 def solve(xs, ys):
    n = len(xs)

    sx1 = sum(xs)
    sx2 = sum(i*i for i in xs)
    sx3 = sum(i*i*i for i in xs)
    sx4 = sum(i*i*i*i for i in xs)

    sy1   = sum(yseries)
    sx1y1 = sum(x*y for (x, y) in zip(xs, ys))
    sx2y1 = sum(x*x*y for (x, y) in zip(xs, ys))

    # Coefficient array representing degree 2 LS fit
    A = np.array([ [  n, sx1, sx2],
                   [sx1, sx2, sx3],
                   [sx2, sx3, sx4]
                 ])

    # Ax = b; build b
    b = np.array([sy1, sx1y1, sx2y1])
    # Factor L into LL^
    L = np.linalg.cholesky(A)

    # Ly = b
    y1 = sy1/L[0][0]
    y2 = (b[1] - L[1][0]*y1)/L[1][1]
    y3 = (b[2] - L[2][0]*y1 - L[2][1]*y2)/L[2][2]
    y = np.array([y1, y2, y3])

    # L.Tx = y
    x3 = y3/L.T[2][2]
    x2 = (y2 - L.T[1][2]*x3)/L.T[1][1]
    x1 = (y1 - L.T[0][2]*x3 - L.T[0][1]*x2)/L.T[0][0]
    return np.array([x1, x2, x3])


 if __name__ == "__main__":
    print xseries
    print yseries
    x = solve(xseries, yseries)
    print x
    f  = lambda i: x[0] + x[1]*i + x[2]*i*i
    f1 = lambda i: x[1] + 2.0*x[2]*i
    fs = [f(v) for v in xseries]

    plt.plot(xseries, yseries)
    plt.plot(xseries, fs)
    plt.show()
	"""
	Second-degree polynomial fit of f(xseries) => yseries
	by Choleskification and Ly=b, L^x=y
	"""
	import numpy as np
	import matplotlib.pyplot as plt

	"""
	# Test series: y=x^2
	xseries = [0, 0.5, 1.0, 1.5, 2.0, 2.5]
	yseries = [0, 0.25, 1.0, 2.25, 4.0, 6.25]
	"""
	xseries = (0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5)
	yseries = (11, 4, 4, 6, 9, 12, 17, 18, 14, 18, 20, 20)
	# Good fit to smooth polynomial-like func...
	#yseries = [np.sin(x/12.0*np.pi) for x in range(12)]
	# Bad fit to periodicity...
	#yseries = [np.sin(x/6.0*np.pi) for x in range(12)]


	def solve(xs, ys):
	n = len(xs)

	sx1 = sum(xs)
	sx2 = sum(i*i for i in xs)
	sx3 = sum(iii for i in xs)
	sx4 = sum(iii*i for i in xs)

	sy1 = sum(yseries)
	sx1y1 = sum(x*y for (x, y) in zip(xs, ys))
	sx2y1 = sum(xxy for (x, y) in zip(xs, ys))

	# Coefficient array representing degree 2 LS fit
	A = np.array([ [ n, sx1, sx2],
	[sx1, sx2, sx3],
	[sx2, sx3, sx4]
	])

	# Ax = b; build b
	b = np.array([sy1, sx1y1, sx2y1])
	# Factor L into LL^
	L = np.linalg.cholesky(A)

	# Ly = b
	y1 = sy1/L[0][0]
	y2 = (b[1] - L[1][0]*y1)/L[1][1]
	y3 = (b[2] - L[2][0]y1 - L[2][1]y2)/L[2][2]
	y = np.array([y1, y2, y3])

	# L.Tx = y
	x3 = y3/L.T[2][2]
	x2 = (y2 - L.T[1][2]*x3)/L.T[1][1]
	x1 = (y1 - L.T[0][2]x3 - L.T[0][1]x2)/L.T[0][0]
	return np.array([x1, x2, x3])


	if __name__ == "__main__":
	print xseries
	print yseries
	x = solve(xseries, yseries)
	print x
	f = lambda i: x[0] + x[1]i + x[2]i*i
	f1 = lambda i: x[1] + 2.0x[2]i
	fs = [f(v) for v in xseries]

	plt.plot(xseries, yseries)
	plt.plot(xseries, fs)
	plt.show()