<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Problem Set 12
First the exercises
Step1:  Let us load up a sample dataset.
Step6:  Now construct a KNN classifier
Step7:  Calculate accuracy on this very small subset.
Step8:  Let's time these different methods to see if the "faster_preds" is actually faster
Step11:  Okay now, let us try the clustering algorithm.
Step12:  Let us load the credit card dataset and extract a small dataframe of numerical features to test on.
Step14:  Now let us write our transformation function.
Step15:  Now let us build some simple loss functions for 1d labels.
Step17:  Now let us define the find split function.
Step18:  One hot encode our dataset
Step20:  Test this to see if it is reasonable
Step21:  Test this out.
Step22:  The naive option
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
import keras
from keras.datasets import mnist
Explanation: Problem Set 12
First the exercises:
* Let $\mu=\frac{1}{|S|}\sum_{x_i\in S} x_i$ let us expand 
\begin{align}
\sum_{x_i\in S} ||x_i-\mu||^2 &=\sum_{x_i\in S}(x_i-\mu)^T(x_i-\mu)\
    &= |S|\mu^T\mu+\sum_{x_i\in S}\left( x_i^Tx_i-2\mu^T x_i \right) \
    &= \frac{1}{|S|}\left(\sum_{(x_i,x_j)\in S\times S} x_i^T x_j\right) + \sum_{x_i\in S} \left( x_i^Tx_i-\frac{2}{|S|}\left(\sum_{x_j\in S} x_j^T x_i\right)\right)\
      &= \sum_{x_i\in S} x_i^Tx_i-\frac{1}{|S|}\sum_{(x_i,x_j)\in S\times S} x_j^T x_i\
      &= \frac{1}{2}\left(\sum_{x_i\in S} x_i^Tx_i-\frac{2}{|S|}\sum_{(x_i,x_j)\in S\times S} x_j^T x_i+\sum_{x_j\in S} x_j^Tx_j \right)\
      &= \frac{1}{2|S|}\left(\sum_{(x_i,x_j)\in S\times S} x_i^Tx_i-2\sum_{(x_i,x_j)\in S\times S} x_j^T x_i+\sum_{(x_i,x_j)\in S\times S} x_j^Tx_j \right)\
      &= \frac{1}{2|S|}\sum_{(x_i,x_j)\in S\times S} (x_i-x_j)^T(x_i-x_j)\
      &= \frac{1}{2|S|}\sum_{(x_i,x_j)\in S\times S} ||x_i-x_j||^2
\end{align}
as desired.
* So the $K$-means algorithm consists of iterations of two steps, we will show that either the algorithm has stabilized or that each of these steps decreases
[ T=\sum_{c=1}^K \sum_{x_i\in S_c}||x_i-\mu_c||^2,] where $S_c$ is the $c$th cluster and $\mu_c$ is the previously defined mean over that cluster. The sequence defined by these sums is therefore monotonically decreasing and bounded below so it will eventually approach the maximal lower bound. 
The value of $\mu$ that minimizes $\sum_{x_i\in S_c}||x_i-\mu||^2$ is $\frac{1}{|S_c|}\sum_{x_i\in S_c} x_i$ (we can check this by setting the derivative with respect to $\mu$ to zero). So updating the mean estimates will never increase $T$. If we do not update the mean estimates than the cluster assignments will not change on the next step. 
The next step maps samples to their closest mean which can only decrease the sum $T$.
Python Lab
Now let us load our standard libraries.
End of explanation
(x_train, y_train), (x_test, y_test) = mnist.load_data()
x_train.shape
import matplotlib.pyplot as plt
%matplotlib inline
randix = np.random.randint(0,60000)
plt.imshow(x_train[randix])
print("Label is {}.".format(y_train[randix]))
x_train_f = x_train.reshape(60000,-1)
x_train_f.shape
x_test_f = x_test.reshape(-1, 28**2)
x_test_f.shape
from sklearn.preprocessing import OneHotEncoder as OHE
ohe = OHE(sparse = False)
y_train_ohe = ohe.fit_transform(y_train.reshape(-1,1))
y_test_ohe = ohe.fit_transform(y_test.reshape(-1,1))
np.argmax(y_train_ohe[randix]) == y_train[randix]
Explanation: Let us load up a sample dataset.
End of explanation
from scipy.spatial.distance import cdist
from sklearn.neighbors import KDTree
class KNNClassifier(object):
    def fit(self,x,y,k=1,fun=lambda x: np.mean(x,axis=0)):
        Fits a KNN regressor.
        Args:
            x (numpy array) Array of samples indexed along first axis.
            y (numpy array) Array of corresponding labels.
            k (int) the number of neighbors
            fun (function numpy array --> desired output) Function to be applied to k-nearest
              neighbors for predictions
        
        self.x = x[:]
        self.y = y[:]
        self.k = k
        self.f = fun
        self.tree = KDTree(self.x)
        
    def predict_one(self, sample):
        
        Run prediction on sample
        
        Args:
            new_x (numpy array) sample
        
        
        dists = cdist(sample.reshape(1,-1),self.x)
        ix = np.argpartition(dists,self.k-1)[0,0:self.k]
        return self.f(self.y[ix])
    
    def predict(self, samples):
        
        Run predictions on list.
        Args:
            samples (numpy array) samples
        
        return np.array([self.predict_one(x) for x in samples])
    
    def faster_predict(self,samples):
        
        Run faster predictions on list.
        Args:
            samples (numpy array) samples
        
        _, ixs = self.tree.query(samples, k=self.k)
        #print(ixs)
        return np.array([self.f(self.y[ix]) for ix in ixs])
        
classifier = KNNClassifier()
classifier.fit(x_train_f, y_train_ohe, k=1)
preds=classifier.predict(x_test_f[:500])
Explanation: Now construct a KNN classifier
End of explanation
np.mean(np.argmax(preds,axis=1)==y_test[:500])
faster_preds = classifier.faster_predict(x_test_f[:500])
np.mean(np.argmax(faster_preds,axis=1)==y_test[:500])
Explanation: Calculate accuracy on this very small subset.
End of explanation
from timeit import default_timer as timer
start = timer()
classifier.predict(x_test_f[:500])
end = timer()
print(end-start)
start = timer()
classifier.faster_predict(x_test_f[:500])
end = timer()
print(end-start)
Explanation: Let's time these different methods to see if the "faster_preds" is actually faster:
End of explanation
def cluster_means(x,cluster_assignments,k):
    
    Return the new cluster means and the within cluster squared distance given the cluster assignments
    
    cluster_counter = np.zeros((k,1))
    cluster_means = np.zeros((k, x.shape[1]))
    for cluster, pt in zip(cluster_assignments, x):
        #print(x)
        cluster_means[cluster] += pt
        cluster_counter[cluster]+=1
    cluster_means = cluster_means/cluster_counter
    
    wcss = 0.
    for cluster, pt in zip(cluster_assignments, x):
        wcss+=np.sum((pt-cluster_means[cluster])**2)
        
    return cluster_means, wcss
class KMeansCluster(object):
    
    #Fit a clustering object on a dataset x consisting of samples on each row
    #by the K-means algorithm into k clusters
    def fit(self,x,k):
        
        Fit k-means clusterer
        
        Args:
            x (numpy array) samples
            k (int) number of clusters
        
        num_samples, num_features = x.shape[0], x.shape[1]
        #Randomly assign clusters
        cluster_assignments = np.random.randint(0,k,num_samples)
        
        #initialize
        cluster_mus = np.zeros((k,num_features))
        
        #update
        new_cluster_mus, wcss = cluster_means(x,cluster_assignments,k)
        count = 1
        while (cluster_mus!=new_cluster_mus).any() and count < 10**3:
            count += 1
            print("Iteration {:3d}, WCSS = {:10f}".format(count,wcss),end="\r")
            cluster_mus = new_cluster_mus
            #calculate distances
            distances = cdist(x,cluster_mus, metric = 'sqeuclidean')
            np.argmin(distances, axis = 1, out = cluster_assignments)
            new_cluster_mus, wcss = cluster_means(x,cluster_assignments,k)
        self.cluster_means = cluster_means
        self.cluster_assignments = cluster_assignments
        self.x = x[:]
        self.wcss = wcss
clusterer = KMeansCluster()
clusterer.fit(x_train_f,10)
clusterer2 = KMeansCluster()
clusterer2.fit(x_train_f,10)
from sklearn.metrics import confusion_matrix
confusion_matrix(y_train, clusterer2.cluster_assignments)
cluster_samples = clusterer2.x[clusterer2.cluster_assignments == 0]
plt.imshow(cluster_samples[0].reshape(28,28))
plt.imshow(cluster_samples[1].reshape(28,28))
plt.imshow(cluster_samples[23].reshape(28,28))
plt.imshow(cluster_samples[50].reshape(28,28))
np.mean(classifier.faster_predict(cluster_samples),axis=0)
Explanation: Okay now, let us try the clustering algorithm.
End of explanation
big_df = pd.read_csv("UCI_Credit_Card.csv")
big_df.head()
len(big_df)
len(big_df.dropna())
df = big_df.drop(labels = ['ID'], axis = 1)
labels = df['default.payment.next.month']
df.drop('default.payment.next.month', axis = 1, inplace = True)
num_samples = 25000
train_x, train_y = df[0:num_samples], labels[0:num_samples]
test_x, test_y = df[num_samples:], labels[num_samples:]
test_x.head()
train_y.head()
Explanation: Let us load the credit card dataset and extract a small dataframe of numerical features to test on.
End of explanation
class bin_transformer(object):
    
    def __init__(self, df, num_quantiles = 2):
        #identify list of quantiles
        self.quantiles = df.quantile(np.linspace(1./num_quantiles, 1.-1./num_quantiles,num_quantiles-1))
            
    
    def transform(self, df):
        
        Args:
            df (pandas dataframe) : dataframe to transform
        
        Returns:
            new (pandas dataframe) : new dataframe where for every feature of the original there will be 
                num_quantiles-1 features corresponding to whether or not the original values where greater
                than or equal to the corresponding quantile.
            fns (dictionary (string,float)) returns dictionary of quantiles
                
        
        new = pd.DataFrame()
        fns = {}
        for col_name in df.axes[1]:
            for ix, q in self.quantiles.iterrows():
                quart = q[col_name]
                new[col_name+str(ix)] = (df[col_name] >= quart)
                fn = quart
                fns[col_name+str(ix)] = [col_name, fn]
        return new, fns
transformer = bin_transformer(train_x,2)
train_x_t, tr_fns = transformer.transform(train_x)
test_x_t, test_fns = transformer.transform(test_x)
train_x_t.head()
Explanation: Now let us write our transformation function.
End of explanation
def bdd_cross_entropy(pred, label):
    return np.mean(-np.sum(label*np.log(pred+10**(-8)),axis=1))
def MSE(pred,label):
    return np.mean(np.sum((pred-label)**2, axis=1))
def acc(pred,label):
    return np.mean(np.argmax(pred,axis=1)==np.argmax(label, axis=1))
def SSE(x,y):
    return np.sum((x-y)**2)
def gini(x,y):
    return 1-np.sum(np.mean(y,axis=0)**2)
Explanation: Now let us build some simple loss functions for 1d labels.
End of explanation
def find_split(x, y, loss, verbose = False):
    
    Args:
        x (dataframe) : dataframe of boolean values
        y (dataframe (1 column)) : dataframe of labeled values
        loss (function: (yvalue, dataframe of labels)-->float) : calculates loss for prediction of yvalue
            for a dataframe of true values.
        verbose (bool) : whether or not to include debugging info
    
    min_ax = None
    N = x.shape[0]
    base_loss = loss(np.mean(y,axis=0),y)
    min_loss = base_loss
    for col_name in x.axes[1]:
        mask = x[col_name]
        num_pos = np.sum(mask)
        num_neg = N - num_pos
        if num_neg*num_pos == 0:
            continue
        pos_y = np.mean(y[mask], axis = 0)
        neg_y = np.mean(y[~mask], axis = 0)
        l = (num_pos*loss(pos_y, y[mask]) + num_neg*loss(neg_y, y[~mask]))/N
        if verbose:
            print("Column {0} split has improved loss {1}".format(col_name, base_loss-l))
        if l < min_loss:
            min_loss = l
            min_ax = col_name
    return min_ax, min_loss, base_loss-min_loss
        
Explanation: Now let us define the find split function.
End of explanation
ohe = OHE(sparse = False)
train_y_ohe = ohe.fit_transform(train_y.values.reshape(-1,1))
train_y_ohe[0:5],train_y.values[0:5]
test_y_ohe = ohe.transform(test_y.values.reshape(-1,1))
Explanation: One hot encode our dataset
End of explanation
find_split(train_x_t, train_y_ohe, bdd_cross_entropy, verbose = False)
np.mean(train_y_ohe[train_x_t['LIMIT_BAL0.5']],axis=0)
np.mean(train_y_ohe[~train_x_t['LIMIT_BAL0.5']],axis = 0)
np.mean(train_y_ohe,axis=0)
#Slow but simple
class decision_tree(object):
    def __init__(self):
        self.f = None
    
    def fit(self, x,y,depth=5,loss=MSE, minsize = 1, quintiles = 2, verbose = False):
        #Construct default function
        mu = np.mean(y, axis=0)
        
        self.f = lambda a: mu 
        
        # Check our stopping criteria
        if(x.shape[0]<=minsize or depth == 0):
            return 
       
        # transform our data
        tr = bin_transformer(x, quintiles)
        tr_x, fns = tr.transform(x)
        split, split_loss, improvement = find_split(tr_x,y,loss)
        
        if verbose:
            print("Improvement: {}".format(improvement))
        #if no good split was found return
        if split == None:
            return
        
        # Build test function
        col_to_split = fns[split][0]
        splitter = lambda a: (a[col_to_split] >= fns[split][1])
        mask = tr_x[split]
        left = decision_tree()
        right = decision_tree()
        left.fit(x[~mask],y[~mask],depth-1,loss, minsize, quintiles)
        right.fit(x[mask],y[mask],depth-1,loss, minsize, quintiles)
        def g(z):
            if(splitter(z)):
                return right.f(z)
            else:
                return left.f(z)
        self.f = g
    def predict(self, x):
        
        Used for bulk prediction
        
        num_samples = x.shape[0]
        
        return np.array([self.f(x.iloc[ix,:]) for ix in range(num_samples)])
        
Explanation: Test this to see if it is reasonable:
End of explanation
dt = decision_tree()
dt.fit(train_x, train_y_ohe, loss = MSE, minsize = 1, depth = 6, quintiles = 50)
dt.predict(test_x.iloc[0:3,:]), test_y_ohe[0:3]
preds = dt.predict(train_x)
np.mean(np.argmax(preds, axis=1)==train_y)
Explanation: Test this out.
End of explanation
1-np.mean(test_y)
class gradient_boosting_trees(object):
    
    def fit(self, x, y, depth = 2, quintiles = 10, num_trees  = 10):
        self.forest = [None]*num_trees
        cur_y = y[:]
        
        for ix in range(num_trees):
            self.forest[ix] = decision_tree()
            self.forest[ix].fit(x, cur_y, loss=MSE, depth = depth, quintiles = quintiles, minsize = 1)
            
            preds =  self.forest[ix].predict(x)
            cur_y = cur_y - preds
    
    def predict(self,x):
        s = 0.
        preds = [tree.predict(x) for tree in self.forest]
        for t in preds:
            s+=t
        return s
    
    
forest = gradient_boosting_trees()
train_y_ohe = ohe.fit_transform(train_y.values.reshape(-1,1))
forest.fit(train_x, train_y_ohe, depth = 20, num_trees = 5, quintiles = 20)
forest.predict(test_x.iloc[0:3,:]), test_y_ohe[0:3]
for_preds = forest.predict(train_x)
for_preds[0:5,:]
train_y_ohe[0:3]
np.mean(np.argmax(for_preds, axis=1)==train_y)
for_preds = forest.predict(test_x)
np.mean(np.argmax(for_preds, axis=1)==test_y)
from sklearn import tree
sktree = tree.DecisionTreeClassifier(max_depth=20)
sktree.fit(train_x, train_y_ohe)
Explanation: The naive option:
End of explanation
<END_TASK>