%%hylang

(* 2 (reduce + (range 1 7) 0))

import matplotlib
matplotlib.use('nbagg')
%matplotlib inline

%%hylang

(import [functools [reduce]]
        [pprint [pprint]]
        [cytoolz [itertoolz]]
        [numpy :as np]
        [matplotlib :as mpl]
        [matplotlib.pyplot :as plt]
        [seaborn :as sns])

%%hylang

(def palette_name "husl")
(def colors (sns.color_palette palette_name 8))
(colors.reverse)
(def cmap (mpl.colors.LinearSegmentedColormap.from_list palette_name colors))

%%hylang

(def data (apply np.genfromtxt ["filip.csv"] {"dtype" float "delimiter" "," "names" True}))

%%hylang

(def x (get data "x"))
(def y (get data "y"))
(def [x-min x-max] [(x.min) (x.max)])

%%hylang

x

%%hylang

y

%%hylang

(apply plt.scatter [x y])
(plt.show)

%%hylang

(def ones (np.ones_like x))
(def size (* np.pi (** (* 20 ones) 2)))

%%hylang

(def figsize {"figsize" [18 18]})
(def fontsize {"fontsize" 24})

%%hylang

(def plot-options {"s" size "c" x "alpha" 0.5 "cmap" cmap})

%%hylang

(apply plt.figure [] figsize)
(apply plt.title ["NIST Data Set for Linear Regression Demo"] fontsize)
(apply plt.xlabel ["$x$ values"] fontsize)
(apply plt.ylabel ["$y$ values"] fontsize)
(apply plt.scatter [x y] plot-options)
(plt.show)

%%hylang

(def coeffs (np.polyfit x y 10))

coeffs

%%hylang

(def [coeffs residuals rank singular-values rcond]
     (apply np.polyfit [x y 10] {"full" true}))

%%hylang

(def model (np.poly1d coeffs))

%%hylang

[(model -9) (model -7) (model -6) (model -5) (model -4) (model -3)]

%%hylang

(def y-predicted (model x))

%%hylang

(def y-mean (/ (np.sum y) y.size))
(def sstot (np.sum (** (- y y-mean) 2)))
(def ssreg (np.sum (** (- y-predicted y-mean) 2)))
(def ssres (np.sum (** (- y y-predicted) 2)))

%%hylang

(def rsquared (- 1 (/ ssres sstot)))

rsquared

%%hylang

(- 0.99672741617239946  0.996727416185620)

%%hylang

(defclass PolynomialLinearModel []
  "A pretty sweet utility Python-Lisp class for creating
  a polynomial curve fitting model")

(defun PolynomialLinearModel.--init-- [self x y degree]
  (setv self.x x)
  (setv self.y y)
  (setv self.degree degree)
  (setv self.results None)
  (setv self.model None)
  (setv [self.coeffs self.residuals self.rank self.singular-values self.rcond]
        [None None None None None])
  (self.polyfit)
  None)

(defun PolynomialLinearModel.get-y-mean [self]
  "Get the mean value of the observed data"
  (/ (np.sum self.y) self.y.size))

(defun PolynomialLinearModel.get-ss-tot [self]
  "Get total sum of the squares"
  (np.sum (** (- self.y (self.get-y-mean)) 2)))

(defun PolynomialLinearModel.get-ss-reg [self]
  "Get the regression sum of squares"
  (np.sum (** (- self.y-predicted (self.get-y-mean)) 2)))

(defun PolynomialLinearModel.get-ss-res [self]
  "Get the sum of squares of residuals"
  (np.sum (** (- self.y self.y-predicted) 2)))

(defun PolynomialLinearModel.get-r-squared [self]
  "Get the R^2 value for the polynomial fit"
  (- 1 (/ (self.get-ss-res) (self.get-ss-tot))))

(defun PolynomialLinearModel.polyfit [self]
  "Do all the business"
  (setv [self.coeffs self.residuals self.rank self.singular-values self.rcond]
        (apply np.polyfit [self.x self.y self.degree] {"full" true}))
  (setv self.model (np.poly1d self.coeffs))
  (setv self.y-predicted (self.model self.x))
  (setv self.r-squared (self.get-r-squared))
  (setv self.results {
        "coeffs" (self.coeffs.tolist)
        "residuals" (self.residuals.tolist)
        "rank" self.rank
        "singular-values" (self.singular-values.tolist)
        "rcond" self.rcond
        "r-squared" self.r-squared}))

(defun PolynomialLinearModel.--str-- [self]
  "Provide a string representation of the data"
  (str self.results))

(defun PolynomialLinearModel.--repr-- [self]
  "Provide a representation of the data"
  (self.--str--))

(defun PolynomialLinearModel.predict [self xs]
  "Given a set of input values, produce outputs using the model"
  (self.model xs))

%%hylang

(def model (PolynomialLinearModel x y 10))

%%hylang

(model.predict -9)

%%hylang

(pprint model.results)

%%hylang

(def fitted-xs (np.linspace x-min x-max 200))
(def fitted-ys (model.predict fitted-xs))

%%hylang

(def [figure axes] (apply plt.subplots [] figsize))
(plt.hold True)
(apply axes.set_title ["NIST Data Set and Polynomial Fit"] fontsize)
(apply axes.set_xlabel ["$x$ values"] fontsize)
(apply axes.set_ylabel ["$y$ values"] fontsize)
(apply axes.scatter [x y] plot-options)
(axes.plot fitted-xs fitted-ys)
(plt.show)