I am trying to find a numerical package which will fit a natural which minimizes weighted least squares.

There is a package in scipy which does what I want for unnatural splines.

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
import random

x = np.arange(0,5,1.0/2)
xs = np.arange(0,5,1.0/500)

y = np.sin(x+1)
for i in range(len(y)):
    y[i] += .2*random.random() - .1

knots = np.array([1,2,3,4])
tck = interpolate.splrep(x,y,s=1,k=3,t=knots,task=-1)
ynew = interpolate.splev(xs,tck,der=0)

plt.figure()
plt.plot(xs,ynew,x,y,'x')