<jupyter_start><jupyter_text>RCS Data Analysis - Curve Fitting<jupyter_code># We will be using SciPY Python library for mathematics,science and engineering
# https://www.scipy.org/ comes with Anaconda
# Most important external Python library NumPy library underneath: http://www.numpy.org/
# Data from the book
# Building Machine Learning Systems with Python
# by <NAME> and <NAME>
# published by PACKT Publishing

# It is made available under the MIT License
import scipy as sp
data = sp.genfromtxt('web_traffic.tsv', delimiter='\t') #default is comma, so we need to specify tab as delimiter
type(data)
data[:12]
# what is our data structure like
data.shape
# We need to clean and preprocess our data a bit:
# split is with index notation (comes from NumPy)
x =data[:, 0] # we grab everything from first dimension(rows) and only the first(0th) column
y = data[:,1]
# how many Nan(invalid entries ) are there ? Too many and we would have a problem
sp.sum(sp.isnan(y))
y[:10]
sp.isnan(y)[:10]
~sp.isnan(y)[:10]
truthtable = sp.isnan(y)
truthtable[:10]
mylist = []
for el in truthtable:
    mylist.append(not el)
mylist[:10]
mlist = [not el for el in truthtable]
mlist[:10]
y[:10]
not sp.isnan(y)
any([True,False,True])
all([True,False,False,True])
not True
True == 1
~True
~1
# we use a negative mask to select only those values which are not NAN
x = x[~sp.isnan(y)]
y = y[~sp.isnan(y)]
len(x),len(y) #sanity check
# Visualizating our data
import matplotlib.pyplot as plt
735 // 24
# plot the (x,y) points with dots of size 10
plt.scatter(x, y, s=10)
plt.title("Web site traffic over the last month")
plt.xlabel("Time(hours)")
plt.ylabel("Hits/h")
plt.xticks([w*7*24 for w in range(10)],
          [f'week {w}' for w in range(10)])
plt.autoscale(tight=True)
plt.grid(True, linestyle='-', color='0.80')
plt.show()<jupyter_output><empty_output><jupyter_text>Question when will our traffic reach 50,000 hits an hour ? (We'll need a new server by then..)  We will need to extrapolate with some model  All models are wrong, but some models are usefulhttps://en.wikipedia.org/wiki/All_models_are_wrong<jupyter_code>#### Error function
# f is our model here
# x is our independent variable
# y is the truth

def error(f, x, y):
    return sp.sum((f(x)-y)**2)

# NumPy thus SciPy strength is vectorized functions, we do all calculations at once or near once, good for larger datasets!<jupyter_output><empty_output><jupyter_text>Could our data be a simple linear progression (straight line )?<jupyter_code>#SciPY has polyfit function
fp1, residuals, rank, sv, rcond = sp.polyfit(x, y, 1, full=True)
print(f"Model Parameters {fp1}")
# 1st degree f(x)= ax+b
# this means f(x) = 2.59619213 * x + 989.02487106 is our linear function
#Error
print(residuals)
# poly1d() to create a model function from the model parameters
f1 = sp.poly1d(fp1)
err1=error(f1, x, y)
err1
plt.scatter(x, y, s=10)
plt.title("Web site traffic over the last month")
plt.xlabel("Time(hours)")
plt.ylabel("Hits/h")
plt.xticks([w*7*24 for w in range(10)],
          [f'week {w}' for w in range(10)])
plt.autoscale(tight=True)
plt.grid(True, linestyle='-', color='0.80')

fx = sp.linspace(0,x[-1], 1000) # generate X-values for plotting
plt.plot(fx, f1(fx), linewidth=3, color="green")
plt.legend(["d={f1.order}"], loc="upper left")

plt.show()<jupyter_output><empty_output><jupyter_text>Hmm first 3-4 weeks seem OK, but then what happens? Is our model really linear?<jupyter_code>## Absolute errors not useful but good for comparing to other values
err1
# Lets try a 2nd order function(square)
f2p = sp.polyfit(x, y, 2)
f2p
# f2(x)=ax^2 + bx + c  in our case 1.05322215e-02, -5.26545650e+00,  1.97476082e+03
f2 = sp.poly1d(f2p) # we need a function from what we just generated
err2 = error(f2, x, y)
err2
err2**0.5
plt.scatter(x, y, s=10)
plt.title("Web site traffic over the last month")
plt.xlabel("Time(hours)")
plt.ylabel("Hits/h")
plt.xticks([w*7*24 for w in range(10)],
          [f'week {w}' for w in range(10)])
plt.autoscale(tight=True)
plt.grid(True, linestyle='-', color='0.80')

fx = sp.linspace(0,x[-1], 1000) # generate X-values for plotting
plt.plot(fx, f1(fx), linewidth=3, color="green")
#plt.legend(["d={f1.order}"], loc="upper left")

plt.plot(fx, f2(fx), linewidth=3, color="yellow")
plt.legend(["d={f1.order}","d={f2.order}"], loc="upper left")

plt.show()
colors = ['g', 'k', 'b', 'm', 'r']
linestyles = ['-', '-.', '--', ':', '-']
def plot_models(x, y, models, fname=None, mx=None, ymax=None, xmin=None, maxhours=None):
    ''' plot input data '''

    plt.figure(num=None, figsize=(10, 8))
    plt.clf()
    plt.scatter(x, y, s=10)
    plt.title("Web traffic over the last month")
    plt.xlabel("Time")
    plt.ylabel("Hits/hour")
    plt.xticks(
        [w * 7 * 24 for w in range(10)], [f'week {w}' for w in range(10)])

    if models:
        if mx is None:
            if maxhours is None:
                mx = sp.linspace(0, x[-1], 1000)
            else:
                mx = sp.linspace(0, maxhours, 1000)
        for model, style, color in zip(models, linestyles, colors):
            # print "Model:",model
            # print "Coeffs:",model.coeffs
            plt.plot(mx, model(mx), linestyle=style, linewidth=2, c=color)

        plt.legend([f"d={m.order}" for m in models], loc="upper left")

    plt.autoscale(tight=True)
    plt.ylim(ymin=0)
    if ymax:
        plt.ylim(ymax=ymax)
    if xmin:
        plt.xlim(xmin=xmin)
    plt.grid(True, linestyle='-', color='0.75')
    if fname != None:
        plt.savefig(fname)
plot_models(x, y, None, "nomodels.png")
models=[f1,f2]
plot_models(x, y, models, "models.png")
f3 = sp.poly1d(sp.polyfit(x, y, 3))
f10 = sp.poly1d(sp.polyfit(x, y, 10))
f15 = sp.poly1d(sp.polyfit(x, y, 15))
models+=[f3,f10]
models+=[f15]
models
plot_models(x,y, models, "mod5.png")
for m in models:
    print(f"{m} \nError: {error(m, x, y)}\n")<jupyter_output>2.596 x + 989 
Error: 317389767.339778

         2
0.01053 x - 5.265 x + 1975 
Error: 179983507.8781792

          3           2
3.05e-05 x - 0.02359 x + 4.942 x + 1333 
Error: 139350144.03172523

           10             9             8             7             6
-3.74e-22 x  + 1.365e-18 x - 2.143e-15 x + 1.899e-12 x - 1.046e-09 x
              5             4           3          2
 + 3.709e-07 x - 8.456e-05 x + 0.01192 x - 0.9416 x + 33.37 x + 1264 
Error: 121942326.36347426

           15            14             13             12            11
-1.75e-33 x  + 9.35e-30 x  - 2.238e-26 x  + 3.173e-23 x  - 2.97e-20 x 
              10             9             8             7
 + 1.938e-17 x  - 9.095e-15 x + 3.127e-12 x - 7.964e-10 x
              6             5            4          3         2
 + 1.509e-07 x - 2.107e-05 x + 0.002094 x - 0.1366 x + 5.027 x - 78.61 x + 1841 
Error: 112482232.88272907<jupyter_text>Inflection point (human look)<jupyter_code>inflection = int(3.5*7*24) # calculate the inflection point in hours
xa = x[:inflection] # data before the inflection point
ya = y[:inflection]
xb = x[inflection:] # data after
yb = y[inflection:]
fa = sp.poly1d(sp.polyfit(xa, ya, 1))
fb = sp.poly1d(sp.polyfit(xb, yb, 1))
fa_error = error(fa, xa, ya)
fb_error = error(fb, xb, yb)

for m in [fa,fb]:
    print(f"{m} \nError: {error(m, x, y)}\n")
## Trusting Error 
plot_models(x,y, [fa,fb], "modab.png")
plot_models(x,y, models)
mx = sp.linspace(0, 1200, 1000)
plot_models(x,y, models, maxhours=1000)<jupyter_output><empty_output><jupyter_text>Bad case of overfitting!! on d10 and d15<jupyter_code>plot_models(x,y, [fa,fb], maxhours=1000)
# Solving equations
from scipy.optimize import fsolve
print(fb)
print(fb - 500000)
reached_max = fsolve(fb - 50000, x0=800) / (7 * 24) # x0 is the starting estimate for speedup
print(f"50,000 hits/hour expected in week {reached_max[0]}")
reached_max
# Solve documentation https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.fsolve.html<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15081967213114755, 'fraction_numerical': 0.07364438839848676, 'mean_word_length': 3.3721058434399116, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 26, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '209586', 'license': 'MIT', 'n_tokens_mistral': 3437, 'n_tokens_neox': 3053, 'n_words': 866, 'path': 'scikit-learn/Curve Fitting - Regression_in_class_02.2020.ipynb', 'repo_name': 'ValRCS/RCS_ML_01_20', 'size': '7948'}