<file_sep>import tensorflow as tf
from tensorflow.examples.tutorials.mnist import input_data
import argparse

LOGDIR = "/tmp/mnist_tutorial/"
mnist = input_data.read_data_sets("/tmp/data/", one_hot=True)

def weight_variable(shape):
  '''Create a weight variable with initialization'''
  return tf.Variable(tf.truncated_normal(shape, stddev=0.1), name="W")

def bias_variable(shape):
  '''Create a bias variable with initialization'''
  return tf.Variable(tf.constant(0.1, shape=shape), name="B")

def conv2d(x, weights):
  '''Perform 2d convolution with specified weights'''
  return tf.nn.conv2d(input=x, filter=weights, strides=[1,1,1,1], padding='SAME')

def max_pool(x):
  '''Perform max pooling'''
  return tf.nn.max_pool(x, ksize=[1,2,2,1], strides=[1,2,2,1], padding='SAME')

def conv_layer(input, kernel_size, in_channels, out_channels, batch_norm=False, name='conv'):
  '''Convolutional layer with 2d convolution, relu activation, and max pooling'''
  with tf.name_scope(name):
    w = weight_variable([kernel_size, kernel_size, in_channels, out_channels])
    b = bias_variable([out_channels])
    if not batch_norm:
        activation = tf.nn.relu(conv2d(input, w) + b)
    else:
        activation = tf.nn.relu(tf.nn.batch_normalization(conv2d(input, w)+b, 0,1,0.1,10,0.0001))
    tf.summary.histogram("Weights", w)
    tf.summary.histogram("Biases", b)
    tf.summary.histogram("Activations", activation)
    pool = max_pool(activation)
    return pool

def fc_layer(input, in_channels, out_channels, name='fc'):
  '''Fully connected layer'''
  with tf.name_scope(name):
    w = weight_variable([in_channels, out_channels])
    b = bias_variable([out_channels])
    out = tf.matmul(input, w) + b
    tf.summary.histogram("Weights", w)
    tf.summary.histogram("Biases", b)
    tf.summary.histogram("Activations", out)
    return out

def train():
  parser = argparse.ArgumentParser()
  parser.add_argument('--test', default=False, type=bool)
  parser.add_argument('--batch_norm', default=False, type=bool)
  args = parser.parse_args()

  tf.reset_default_graph()
  sess = tf.Session()
  # Input flattened image
  x = tf.placeholder(tf.float32, shape=[None, 784], name='x')
  # Class label
  y = tf.placeholder(tf.float32, shape=[None, 10], name='y')

  x_image = tf.reshape(x, [-1,28,28,1])

  # Model
  conv_out1 = conv_layer(x_image, 5, 1, 32, args.batch_norm, "conv1")
  conv_out2 = conv_layer(conv_out1, 5, 32, 64, args.batch_norm, "conv2")
  flattened_out = tf.reshape(conv_out2, [-1, 7*7*64])
  fc_out1 = fc_layer(flattened_out, 7*7*64, 1024, "fc1")
  relu = tf.nn.relu(fc_out1)
  dropout = tf.layers.dropout(inputs=relu, rate=0.4)
  logits = fc_layer(dropout, 1024, 10, "fc2")

  loss = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=logits, labels=y))
  tf.summary.scalar("Loss", loss)

  optimizer = tf.train.AdamOptimizer(0.001)
  train_step = optimizer.minimize(loss)

  correct_prediction = tf.equal(tf.argmax(logits,1), tf.argmax(y,1))
  accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
  train_acc = tf.summary.scalar('Accuracy', accuracy)
  validation_acc = tf.summary.scalar('Validation_accuracy', accuracy)

  merged_summary = tf.summary.merge_all()

  sess.run(tf.global_variables_initializer())
  writer = tf.summary.FileWriter('/tmp/mnist_tutorial')
  writer.add_graph(sess.graph)

  for epoch in range(2001):
    batch = mnist.train.next_batch(50)
    if epoch % 100 == 0:
      [summary, train_accuracy] = sess.run([train_acc, accuracy], \
                                feed_dict={x: batch[0], y: batch[1]})
      writer.add_summary(summary, epoch)
      print('Epoch %d, train accuracy: %f'%(epoch, train_accuracy))
      [validation_summary, val_accuracy] = sess.run([validation_acc, accuracy], \
                feed_dict={x: mnist.validation.images, y:mnist.validation.labels})
      writer.add_summary(validation_summary, epoch)
      print('Epoch %d, validation accuracy: %f'%(epoch, val_accuracy))
    sess.run(train_step, feed_dict={x: batch[0], y: batch[1]})
    if args.test:
        if epoch == 801:
            break
  if args.test:
    print('Test accuracy: %g'%accuracy.eval(session=sess, \
                feed_dict={x: mnist.test.images, y:mnist.test.labels}))
  writer.close()

train()
<file_sep>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mlxtend.preprocessing import standardize
from sklearn.cross_validation import train_test_split
import random
import argparse
from process import *

parser = argparse.ArgumentParser()
parser.add_argument('--num_epochs', default=50, type=int)
parser.add_argument('--num_steps', default=300, type=int)
parser.add_argument('--eval_steps', default=30, type=int, help='compute the accuracy of the current classifier on the set held out for the epoch every 30 steps')
parser.add_argument('--batch_size', default=1, type=int)
parser.add_argument('--m', default=1, type=float, help='parameter for learning rate')
parser.add_argument('--n', default=0.1, type=float, help='parameter for learning rate')
parser.add_argument('--lamda', default=0.001, type=float, help='regularization parameter')
parser.add_argument('--num_held_out', default=50, type=int)
parser.add_argument('--test', default=False, type=bool)
parser.add_argument('--path', default='yuchecw2_prediction.csv', type=str)
args = parser.parse_args()

train_data = pd.read_csv("train.data.csv")
label_list = [i for i in train_data['class']]
label = np.empty(len(label_list))
for i in range(label.shape[0]):
    if (label_list[i] == ' >50K'):
        label[i] = 1
    else:
        label[i] = -1

feature = pd.DataFrame.as_matrix(train_data)[:, :-1]
# Preprocess the data
feature = preprocess(feature)
# Standardize columns in train_feature
feature = standardize(feature)
# Train-validation split
train_feature, val_feature, train_label, val_label = train_test_split(
                feature, label, test_size=0.1)

plt.figure()
for lamda in [0.001, 0.01, 0.1, 1]:
    # Initialize a, b
    a = np.random.randn(train_feature.shape[1], 1)  # shape=(14, 1)
    b = 0
    accuracy_list = []
    mag_list = []
    loss_list = []
    for epoch in range(args.num_epochs):
        actual_train_feature, held_out_feature, actual_train_label, held_out_label = train_test_split(train_feature, train_label, test_size=args.num_held_out/train_feature.shape[0], random_state=epoch)
        actual_train_size = actual_train_feature.shape[0]
        for step in range(args.num_steps):
            lr = compute_lr(args.m, args.n, epoch)
            batch_num = random.sample(range(actual_train_size), args.batch_size)
            a = update_a(a, b, actual_train_feature, actual_train_label, batch_num, lr, lamda)
            b = update_b(b, a, actual_train_feature, actual_train_label, batch_num, lr, lamda)
            if (step%args.eval_steps == 0):
                eval_dict = evaluate(held_out_feature, held_out_label, a, b, lamda)
                accuracy_list.append(eval_dict['accuracy'])
                mag_list.append(eval_dict['mag'])
                loss_list.append(eval_dict['loss'])

    if (args.test):
        test_data = pd.read_csv("test.data.csv")
        test_feature = pd.DataFrame.as_matrix(test_data)
        test_feature = preprocess(test_feature)
        test_feature = standardize(test_feature)
        prediction = predict(test_feature, a, b)
        prediction[prediction['Label'] == '>'] = ">50K"
        prediction[prediction['Label'] == '<'] = "<=50K"
        prediction.to_csv(args.path)

    step_list = range(len(accuracy_list))

    for step in step_list:
        step *= args.num_steps


    label = 'λ = {}'.format(lamda)
    plt.plot(step_list, mag_list, label=label, linewidth=0.8)
    plt.ylim(0, 100)

plt.xlabel('step number')
plt.ylabel('magnitude')
plt.title('Magnitude vs step number', fontsize=12, fontweight='bold')
plt.legend(loc='lower right')
plt.show()
<file_sep>import os
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.cross_validation import train_test_split
from sklearn.ensemble import RandomForestClassifier
import matplotlib.pyplot as plt
import argparse

def class_confusion_matrix(predict, actual):
    confusion_matrix = np.zeros((14,14)) # len(name_label_dict.keys()) = 14
    for i in range(len(predict)):
        confusion_matrix[predict[i], actual[i]] += 1
    return confusion_matrix

# Can also be used to obtain test_data, test_label
def get_data_label(train_data_dict, data_dict, name_label_dict):
    train_data =  np.vstack(train_data_dict['Brush_teeth'])
    train_label = [name_label_dict['Brush_teeth']]*train_data.shape[0]
    for key in data_dict.keys():
        if key is not 'Brush_teeth':
            current_train_data = np.vstack(train_data_dict[key])
            train_data = np.vstack([train_data, current_train_data])
            train_label += [name_label_dict[key]]*current_train_data.shape[0]
    return train_data, train_label

def build_name_label_dict(data_dict):
    name_label_dict = {}
    for i, key in enumerate(data_dict.keys()):
        name_label_dict[key] = i
    return name_label_dict

def split_data_dict(data_dict, test_size=0.20, random_state=42):
    train_data_dict = {}
    test_data_dict = {}
    for label in data_dict.keys():
        train_data, test_data = train_test_split(data_dict[label], test_size=test_size, random_state=random_state)
        train_data_dict[label] = train_data
        test_data_dict[label] = test_data
    return train_data_dict, test_data_dict

# (label, data(shape: (n_samples, n_features))) pair
def build_data_dict(folders, dim, overlap):
    data_dict = {}
    for folder in folders:
        if folder in data_dict.keys():
            data_dict[folder] += [preprocess('HMP_Dataset/'+folder+'/'+file, dim, overlap) for file in os.listdir('HMP_Dataset/'+folder) if not file.startswith('.')]
        else:
            key_exists = False
            for key in data_dict.keys():
                if folder.startswith(key):
                    key_exists = True
                    data_dict[key] += [preprocess('HMP_Dataset/'+folder+'/'+file, dim, overlap) for file in os.listdir('HMP_Dataset/'+folder) if not file.startswith('.')]
            if not key_exists:
                data_dict[folder] = [preprocess('HMP_Dataset/'+folder+'/'+file, dim, overlap) for file in os.listdir('HMP_Dataset/'+folder) if not file.startswith('.')]
    return data_dict

# Build frequency histogram with length n_clusters
def build_histogram(labels, n_clusters):
    histogram = [0]*n_clusters
    total = 0
    for label in labels:
        histogram[label] += 1
        total += 1
    X = np.array(histogram).reshape(1,-1)
    X = X / total
    return X

# Build a histogram of shape (n_clusters,)
def prediction(kmeans, train_data_dict, n_clusters, name_label_dict):
    predict_hist = []
    predict_label = []
    for key in train_data_dict.keys():
        for data in train_data_dict[key]:
            hist = build_histogram(kmeans.predict(data), n_clusters)
            predict_hist.append(hist)
            predict_label.append(name_label_dict[key])
    predict_hist = np.array(predict_hist).squeeze(1)
    return predict_hist, predict_label

def preprocess(file, dim, overlap):
    signal = np.loadtxt(file)
    signal_patches = cut_signal(signal, dim, overlap)
    return np.vstack(signal_patches)

def cut_signal(signal, dim, overlap):
    i=0
    signal_patches = []
    while(i+dim < signal.shape[0]):
        signal_patches.append(signal[i:i+dim, :].reshape(1, -1))
        i += (int)(dim*(1-overlap)) # no overlap
    if i < signal.shape[0]:
        signal_patches.append(signal[-dim:, :].reshape(1, -1))
    return signal_patches

def main(dim, n_clusters, overlap, fig):
    folders = os.listdir('HMP_Dataset')
    folders.pop(0) # Remove '.DS_Store'

    data_dict = build_data_dict(folders, dim, overlap)
    name_label_dict = build_name_label_dict(data_dict)

    train_data_dict, test_data_dict = split_data_dict(data_dict)
    train_data, train_label = get_data_label(train_data_dict, data_dict, name_label_dict) #train_data[i,:] -> train_label[i]
    test_data, test_label = get_data_label(test_data_dict, data_dict, name_label_dict)

    kmeans = KMeans(n_clusters=n_clusters, init='k-means++', n_init=5).fit(train_data)
    predict_hist, predict_label = prediction(kmeans, train_data_dict, n_clusters, name_label_dict)
    clf = RandomForestClassifier(n_estimators=20, max_depth=10, random_state=0)
    clf.fit(predict_hist, predict_label)

    #Test
    test_predict_hist, test_actual_label = prediction(kmeans, test_data_dict, n_clusters, name_label_dict)
    test_predict_label = clf.predict(test_predict_hist)
    n_correct = len(test_predict_label[test_predict_label==test_actual_label])
    accuracy = n_correct / len(test_predict_label)
    print('Accuracy: {}'.format(accuracy))
    confusion_matrix = class_confusion_matrix(test_predict_label, test_actual_label)

    for key_num, key in enumerate(train_data_dict.keys()):
        sample = np.zeros((1,dim*3))
        for i in range(len(train_data_dict[key])):
            for j in range(train_data_dict[key][i].shape[0]):
                sample = np.vstack((sample, train_data_dict[key][i][j,:].reshape(1,dim*3)))

        hist = build_histogram(kmeans.predict(sample[1:,:]), n_clusters)
        ax = fig.add_subplot(3,5,key_num+1)
        ax.bar(np.arange(1, n_clusters+1)-0.4, hist.squeeze(0).tolist(), width=0.8)
        ax.set_xlabel('Cluster number', fontsize=8)
        ax.set_ylabel('Frequency', fontsize=8)
        ax.set_title('Histogram for {}'.format(key), fontsize=8)
        ax.grid()
    print(confusion_matrix)

fig = plt.figure(figsize=(13,7))
confusion_matrix = main(dim=5, n_clusters=50, overlap=0.5, fig=fig)
fig.tight_layout()
plt.show()
# np.savetxt("confusion_matrix.csv", confusion_matrix, delimiter='')
<file_sep>import numpy as np
import pandas as pd
from sklearn.naive_bayes import BernoulliNB
from sklearn.naive_bayes import GaussianNB
import random
import math

def preprocess(data, test):
    box_size = 20
    original_size = 28
    if not test:
        streched_data = np.empty((data.shape[0], box_size*box_size+1))
    else:
        # No label column
        streched_data = np.empty((data.shape[0], box_size*box_size))
    for k in range(data.shape[0]):
        # Label
        if not test:
            streched_data[k, 0] = data[k, 0]
        # Find the margin
        bottom = find_bottom(data, k, original_size)
        top = find_top(data, k, original_size)
        num_rows = bottom-top+1
        right = find_right(data, k, original_size, test)
        left = find_left(data, k, original_size, test)
        num_cols = right-left+1
        # Find the bounding box
        bounding_box = np.empty((num_rows, num_cols))
        for j in range(top, bottom+1):
            for i in range(left, right+1):
                if not test:
                    bounding_box[j-top, i-left] = data[k, original_size*j+i+1]
                else:
                    bounding_box[j-top, i-left] = data[k, original_size*j+i]
        # plt.subplot(1,2,1)
        # plt.matshow(bounding_box)
        # plt.title("#%d bounding box"%(k))
        # Resize the bounding box
        streched_bounding_box = resize_matrix(bounding_box, box_size, box_size)
        for i in range(box_size):
            for j in range(box_size):
                if not test:
                    streched_data[k, box_size*i+j+1] = streched_bounding_box[i, j]
                else:
                    streched_data[k, box_size*i+j] = streched_bounding_box[i,j]
    return streched_data

def resize_matrix(mat, n_row_out, n_col_out):
    n_row_in = mat.shape[0]
    n_col_in = mat.shape[1]
    mat_out = np.empty((n_row_out, n_col_out))
    row_ratio = n_row_in/n_row_out
    col_ratio = n_col_in/n_col_out
    for i in range(n_row_out):
        for j in range(n_col_out):
            mat_out[i, j] = mat[math.floor(i*row_ratio), math.floor(j*col_ratio)]
    return mat_out

def find_top(data, num, original_size):
    for i in range(1, len(data[num, :])):
        if data[num, i].item() is not 0:
            return math.floor(i/original_size)

def find_bottom(data, num, original_size):
    for i in range(len(data[num, :]) - 1, 0, -1):
        if data[num, i].item() is not 0:
            return math.floor(i/original_size)

def find_left(data, num, original_size, test):
    for left in range(original_size):
        for j in range(original_size):
            if not test:
                if data[num, original_size*j+left+1].item() is not 0:
                    return left
            else:
                if data[num, original_size*j+left].item() is not 0:
                    return left

def find_right(data, num, original_size, test):
    for right in range(original_size-1,-1,-1):
        for j in range(original_size):
            if not test:
                if data[num, original_size*j+right+1].item() is not 0:
                    return right
            else:
                if data[num, original_size*j+right].item() is not 0:
                    return right
<file_sep># Part 1d
library(caret)
library(e1071)
library(klaR)
# Read CSV into R
data = read.csv(file="pima-indians-diabetes.csv", header=TRUE, sep=",")
attributes = c(3, 4, 6, 8) #blood pressure, skin thickness, BMI, age
accuracy_vec = seq(1, 10, by=1)
for(seed in seq(1, 10)) {
  # Test train split
  set.seed(seed)
  in_train = createDataPartition(data[,9], p=0.8, list=FALSE)
  train_data = data[in_train, ]
  test_data = data[-in_train, ]
  test_label = test_data[,9]
  # Train the model
  feature = train_data[, attributes]
  label = as.factor(train_data[,9])
  model = svmlight(feature, label, pathsvm="/Users/jeffreywang/Desktop/AML/svm_light")
  prediction = predict(model, test_data[, attributes])
  predict_label = prediction$class
  correct = test_label == predict_label
  
  accuracy = length(correct[correct==TRUE])/length(correct)
  accuracy_vec[seed] = accuracy
}
print(accuracy_vec)
print(mean(accuracy_vec))
<file_sep>setwd('/Users/jeffreywang/Desktop/AML/HW7')
library(gdata)
data <- read.csv('blogData_train.csv')
library(glmnet)
library("rlist")
library(data.table)
xmat = as.matrix(data[,-281])
ymat = as.matrix(data[,281])
# Poisson generalized model with lasso
model <- cv.glmnet(xmat, ymat, family='poisson',nfolds=10, alpha=1)
# A plot of cross-validated deviance against the regularization variable
plot(model)
# A scatter plot of true values vs predicted values for training data
predicted = predict(model, newx=xmat, s="lambda.1se")
plot(predicted, ymat, main="Training data\nPredicted v.s. True", xlab="Predicted value", ylab="True value")
reg <- lm(ymat~predicted) # R squared = 0.1398
abline(reg, col='red')
# A scatter plot of true values vs predicted values for testing data
test_data_files = list.files(path="BlogFeedback")
# Load test data from csv files
test_data <- read.csv(paste("BlogFeedback/", test_data_files[1], sep=""))
for (i in 2:60) {
  test_data = rbindlist(list(test_data, read.csv(paste("BlogFeedback/", test_data_files[i], sep=""))))
}
# Split test data to x_test and y_test
par(mfrow=c(1,1))
x_test = as.matrix(test_data[,-281])
y_test = as.matrix(test_data[,281])
predicted_test = predict(model, x_test, s="lambda.1se")
plot(predicted_test, y_test, main="Testing data\nPredicted v.s. True", xlab="Predicted value", ylab="True value")
reg_test <- lm(y_test~predicted_test)
abline(reg_test, col='red')

<file_sep>import numpy as np
from PIL import Image
from sklearn.cluster import KMeans
from collections import Counter
import matplotlib.pyplot as plt

def load_image(filename):
    img = Image.open(filename)
    img_arr = np.asarray(img, dtype='int32')
    img_ravel = np.concatenate((img_arr[:,:,0].ravel().reshape(-1,1), img_arr[:,:,1].ravel().reshape(-1,1), img_arr[:,:,2].ravel().reshape(-1,1)), axis=1)
    return img_arr, img_ravel, img_scaled # (H*W, C)

# Perform EM Algorithm to segement image
def image_segmentation(image_name, filename, n_segments, random_seed):
    # Initialization
    orig_img, unscaled_img, img = load_image(filename)
    nrows, ncols, _ = orig_img.shape
    kmeans = KMeans(n_clusters=n_segments, random_state=random_seed).fit(img)
    blob_weight = np.zeros(n_segments)   # pi, i.e. weight of each blob
    counter = Counter(kmeans.labels_)
    for i in range(n_segments):
        blob_weight[i] = counter[i] / img.shape[0]
    blob_center = kmeans.cluster_centers_ # mu, i.e. mean of each normal distribution
    w = np.zeros((img.shape[0], n_segments))
    Q_prev = 0
    num_epochs = 30
    epsilon = 0.1
    # EM Algorithm
    for epoch in range(num_epochs):
        # E step
        for i in range(img.shape[0]):
            total = 0
            for j in range(n_segments):
                w[i,j] = np.exp(-0.5*(np.dot(img[i,:]-blob_center[j,:], img[i,:]-blob_center[j,:])))*blob_weight[j]
                total += w[i,j]
            w[i,:] = w[i,:] / total
        # M step
        for j in range(n_segments):
            total = np.zeros((1,3))
            for i in range(img.shape[0]):
                total += img[i,]*w[i,j]
            total = total/(np.sum(w[:,j]))
            blob_center[j] = total
            blob_weight[j] = np.sum(w[:,j]) / img.shape[0]
        # Calculate Q
        Q = 0
        for i in range(img.shape[0]):
            for j in range(n_segments):
                Q += (-0.5*(np.dot(img[i,:]-blob_center[j,:], img[i,:]-blob_center[j,:])) + np.log(blob_weight[j]))*w[i,j]
        # Check if converge or not
        if np.abs(Q - Q_prev) < epsilon:
            break
        Q_prev = Q

    # plt.figure()
    # plt.plot(list(range(len(Q_list))), Q_list, 'b-')
    # plt.xlabel('Epoch')
    # plt.ylabel('Q')
    # plt.title('Q v.s. Epoch {}-{}'.format(image_name, n_segments))
    # plt.show()

    cluster_index = [np.argsort(w[i,:])[-1] for i in range(img.shape[0])]
    segmented_img_ravel = np.zeros(img.shape)
    for i in range(img.shape[0]):
        segmented_img_ravel[i,:] = blob_center[cluster_index[i],:]
    segmented_img = np.zeros(orig_img.shape)
    for i in range(img.shape[0]):
        segmented_img[i//ncols, i%ncols, :] = segmented_img_ravel[i,:]
    segmented_img = Image.fromarray(segmented_img.astype(np.uint8),'RGB')
    segmented_img.save('{}_segment-{}_seed-{}.jpg'.format(image_name, n_segments, random_seed))

def main():
    image_name = ['RobertMixed03','smallsunset','smallstrelitzia', 'tree']
    filename = [i+'.jpg' for i in image_name]
    n_segments_list = [10,20,50]

    for i in range(4):
        for n_segments in n_segments_list:
            image_segmentation(image_name[i], filename[i], n_segments, 0)
            print('Image: {}, num_segments: {}'.format(image_name[i], n_segments))

    for random_seed in [1,2,3,4,5]:
        image_segmentation(image_name[3], filename[3], 20, random_seed)
        print('Tree image with random seed: {}'.format(random_seed))

main()
<file_sep>import numpy as np
import gzip
import struct
import matplotlib.pyplot as plt
from mlxtend.data import loadlocal_mnist

def extract_images(filename, num_images):
    with gzip.open(filename) as f:
        zero, data_type, dims = struct.unpack('>HBB', f.read(4))
        shape = tuple(struct.unpack('>I', f.read(4))[0] for d in range(dims))
        images = np.fromstring(f.read(), dtype=np.uint8).reshape(shape)
        images = np.divide(images, 255.0)
    return images[:num_images, :, :]

def binarize(images):
    bi_images = np.zeros(images.shape)
    for k in range(images.shape[0]):
        for i in range(images.shape[1]):
            for j in range(images.shape[2]):
                if images[k,i,j] >= 0.5:
                    bi_images[k,i,j] = 1.0
                else:
                    bi_images[k,i,j] = -1.0
    return bi_images

def create_noisy(images):
    noisy_images = np.copy(images)
    size = images.shape[1] * images.shape[2]
    flip_size = (int)(size*0.02)
    for k in range(images.shape[0]):
        choice = np.random.permutation(np.arange(size))[:flip_size].tolist()
        for i in range(size):
            if i in choice:
                noisy_images[k, i//images.shape[1], i%images.shape[1]] \
                            = -images[k, i//images.shape[1], i%images.shape[1]]
    return noisy_images

def denoise(noisy_images, theta_hh=0.2):
    # theta_hh = 0.2
    theta_hx = 0.2
    epsilon = 0.001
    num_epochs = 20
    # Initialize diff list
    diff = [[] for _ in range(noisy_images.shape[0])]
    for i in diff:
        i.append(0)
    length = noisy_images.shape[1]
    images = np.copy(noisy_images)
    for k in range(images.shape[0]):
        # Initialize edge weights pi
        pi = np.random.rand(length, length)
        prev_pi = np.copy(pi)
        for epoch in range(num_epochs):
            exponent = np.zeros((length, length))
            for i in range(images.shape[1]):
                for j in range(images.shape[2]):
                    if i is not 0: # Not on the top edge
                        exponent[i,j] += theta_hh*(2*pi[i-1,j]-1) + theta_hx*noisy_images[k,i-1,j]
                    if i is not images.shape[1]-1: # Not on the bottom edge
                        exponent[i,j] += theta_hh*(2*pi[i+1,j]-1) + theta_hx*noisy_images[k,i+1,j]
                    if j is not 0: # Not on the left edge
                        exponent[i,j] += theta_hh*(2*pi[i,j-1]-1) + theta_hx*noisy_images[k,i,j-1]
                    if j is not images.shape[1]-1: # Not on the right edge
                        exponent[i,j] += theta_hh*(2*pi[i,j+1]-1) + theta_hx*noisy_images[k,i,j+1]
                    # Update edge weights
                    pi[i,j] = np.exp(exponent[i,j]) / (np.exp(exponent[i,j]) + np.exp(-exponent[i,j]))
                    if pi[i,j] < 0.5:
                        images[k,i,j] = -1.0
                    else:
                        images[k,i,j] = 1.0

            # diff[k].append(np.linalg.norm(pi-prev_pi,2))
            diff[k].append(np.sum(np.power(pi-prev_pi,2)))
            prev_pi = np.copy(pi)
            if diff[k][-1] < epsilon:
                break
    return images

def accuracy(binary_images,denoise_images):
    accuracy_list = np.zeros((binary_images.shape[0],1))
    for k in range(binary_images.shape[0]):
        n_incorrect = np.count_nonzero(binary_images[k,:,:]-denoise_images[k,:,:])
        accuracy_list[k] = 1 - (n_incorrect / (binary_images.shape[1]*binary_images.shape[2]))
    return accuracy_list

def confusion(binary_images, denoise_images):
    true_positive_list = np.zeros((binary_images.shape[0],1))
    false_positive_list = np.zeros((binary_images.shape[0],1))
    for k in range(binary_images.shape[0]):
        true_positive = 0
        false_positive = 0
        for i in range(binary_images.shape[1]):
            for j in range(binary_images.shape[2]):
                if denoise_images[k,i,j] == 1.0:
                    if binary_images[k,i,j] == 1.0:
                        true_positive += 1
                    else:
                        false_positive += 1
        true_positive_list[k] = true_positive / (binary_images.shape[1]**2)
        false_positive_list[k] = false_positive / (binary_images.shape[1]**2)
    return np.mean(true_positive_list), np.mean(false_positive_list)

def main():
    # img_filename = 'train-images-idx3-ubyte.gz'
    # images = extract_images(img_filename, 500)
    images, labels = loadlocal_mnist(images_path='train-images-idx3-ubyte', labels_path='train-labels-idx1-ubyte')
    images = images.reshape(-1,28,28)[:500,:,:]
    labels = labels[:500]
    label_img_dict = {}
    for i in range(labels.shape[0]):
        if labels[i] not in label_img_dict.keys():
            label_img_dict[labels[i]] = images[i,:,:].reshape(-1,28,28)
        else:
            label_img_dict[labels[i]] = np.concatenate((label_img_dict[labels[i]], images[i,:,:].reshape(-1,28,28)), axis=0)

    # Sample images for each digit
    plt.figure(figsize=(4.5,15))
    for num in range(10):
        orig_img = label_img_dict[num]
        binary_img = binarize(orig_img)
        plt.subplot(10,3,num*3+1)
        plt.imshow(binary_img[0,:,:], cmap='gray')
        noisy_img = create_noisy(binary_img)
        plt.subplot(10,3,num*3+2)
        plt.imshow(noisy_img[0,:,:], cmap='gray')
        denoise_img = denoise(noisy_img)
        plt.subplot(10,3,num*3+3)
        plt.imshow(denoise_img[0,:,:], cmap='gray')
    plt.savefig('samples.png')

    binary_images = binarize(images)
    noisy_images = create_noisy(binary_images)
    denoise_images = denoise(noisy_images)
    accuracy_list = accuracy(binary_images, denoise_images)
    avg_accuracy = sum(accuracy_list[:500]) / 500
    print('Average accuracy on the first 500 images: {}'.format(avg_accuracy))

    plt.figure()
    plt.scatter(list(range(accuracy_list.shape[0])), accuracy_list)
    plt.xlabel('Image Number')
    plt.ylabel('Accuracy')
    plt.ylim(0.9, 1.0)
    plt.title('Fraction of correct pixels')
    plt.savefig('accuracy.png')

    # Most accurate
    max_idx = np.argmax(accuracy_list)
    plt.figure()
    plt.imshow(binary_images[max_idx,:,:])
    plt.title('Most accurate binary image')
    plt.savefig('most_accurate_binary_image.png')
    plt.figure()
    plt.imshow(noisy_images[max_idx,:,:])
    plt.title('Most accurate noisy image')
    plt.savefig('most_accurate_noisy_image.png')
    plt.figure()
    plt.imshow(denoise_images[max_idx,:,:])
    plt.title('Most accurate denoised image')
    plt.savefig('most_accurate_denoised_image.png')

    # Least accurate
    min_idx = np.argmin(accuracy_list)
    plt.figure()
    plt.imshow(binary_images[min_idx,:,:])
    plt.title('Least accurate binary image')
    plt.savefig('least_accurate_binary_image.png')
    plt.figure()
    plt.imshow(noisy_images[min_idx,:,:])
    plt.title('Least accurate noisy image')
    plt.savefig('least_accurate_noisy_image.png')
    plt.figure()
    plt.imshow(denoise_images[min_idx,:,:])
    plt.title('Least accurate denoised image')
    plt.savefig('least_accurate_denoised_image.png')

    # ROC
    denoise_images_list = []
    accuracies = []
    true_positive_list = []
    false_positive_list = []
    for theta_hh in [-1,0,0.2,1,2]:
        denoise_images_list.append(denoise(noisy_images, theta_hh))
        accuracies.append(accuracy(binary_images, denoise_images))
        true_positive, false_positive = confusion(binary_images, denoise_images_list[-1])
        true_positive_list.append(true_positive)
        false_positive_list.append(false_positive)

    txt_list = [-1,0,0.2,1,2]
    fig, ax = plt.subplots()
    ax.scatter(false_positive_list, true_positive_list)
    for i, txt in enumerate(txt_list):
        ax.annotate(txt, (false_positive_list[i], true_positive_list[i]))
    ax.set_xlabel('False Positive Rate')
    ax.set_ylabel('True Positive Rate')
    plt.title('Receiver Operating Curve')
    plt.savefig('roc.png')

main()
<file_sep>import numpy as np
from sklearn.decomposition import PCA
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from scipy.spatial.distance import euclidean
# from skbio.stats.ordination import pcoa
from adjustText import adjust_text

def unpickle(file):
    import pickle
    with open(file, 'rb') as fo:
        dict = pickle.load(fo, encoding='bytes')
    return dict

meta_data = unpickle('cifar-10-batches-py/batches.meta')
data_file = ['cifar-10-batches-py/data_batch_%d'%i for i in range(1,6)]

# Key: label, value: a list of images with that label
label_data_dict = {i:[] for i in range(10)}
for file_num in range(1,6):
    original_dict = unpickle(data_file[file_num-1])
    for i, key in enumerate(original_dict[b'labels']):
        label_data_dict[key].append(original_dict[b'data'][i,:])

# Stack all images of a category into a single 2D numpy array
# Key: label, value: a matrix of images with that label
label_matrix_dict = {}
for i in range(10):
    label_matrix_dict[i] = np.empty((len(label_data_dict[i]), 3072))
    for j in range(len(label_data_dict[i])):
        label_matrix_dict[i][j,:] = label_data_dict[i][j]

mean_image_list = []
pca_list = [PCA(n_components=20) for i in range(10)]
recon_image_list = []
mse_list = []
for i in range(10):
    mean_image_list.append(np.mean(label_matrix_dict[i], axis=0))
    pca_list[i].fit(label_matrix_dict[i])
    recon_image_list.append(pca_list[i].inverse_transform(pca_list[i].transform(label_matrix_dict[i])))
    mse_list.append(mean_squared_error(label_matrix_dict[i], recon_image_list[i]))

plt.figure(1)
plt.bar(list(range(10)), mse_list)
plt.xlabel('Label')
plt.ylabel('MSE Loss')
plt.title('Error of each category')
plt.show()

distance_matrix = np.empty((10, 10))
for i in range(10):
    for j in range(10):
        distance_matrix[i,j] = euclidean(mean_image_list[i], mean_image_list[j])

# Principle Coordinate Analysis
A = np.identity(10) - 1/10 * np.ones((10,10))
W = -1/2 * np.dot(np.dot(A, distance_matrix), A.T)
eigen_value, eigen_vector = np.linalg.eig(W)
y = np.dot(eigen_vector[:,0:2], np.sqrt(np.diag(eigen_value[0:2])))

labels = np.array([str(x)[2:-1] for x in meta_data[b'label_names']])
fig, ax = plt.subplots()
plt.plot(y[:,0], y[:,1], 'bo')
texts = [plt.text(y[i,0], y[i,1], labels[i], ha='center', va='center') for i in range(labels.shape[0])]
adjust_text(texts)
plt.title('2D Map obtained from PCoA')
plt.show()

# ordination_result = pcoa(distance_matrix)
# ordination_result.proportion_explained
<file_sep>Part 1
Run hw1a.R, hw1b.R, hw1d.R by
$ Rscript hw1a.R
to obtain the accuracy of each prediction

Part2
Run hw2a.py by
$ python hw2a.py
To change models and preprocess images, specify it by adding command line arguments.
Options:
--model_name   ['Gaussian', 'Bernoulli', 'Random Forest'] 	Select a model for training.
--preprocess   bool						Preprocess the image by finding its bounding box and crop it.
--num_epochs   int						Number of epochs where each epoch has different train-validatio splits.
--test	       bool						Test or not. If true, prediction labels will be written into a csv file and mean images will be shown.
--path         str						Specify the path where the csv file containing prediction labels is stored.
--n_estimators int						Number of estimators in the random forest model.
--max_depth    int						The max depth for each tree in the random forest model.

<file_sep>import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error
from sklearn.decomposition import PCA
import argparse

parser = argparse.ArgumentParser()
parser.add_argument('--data_file', default='hw3-data/iris.csv', type=str, help='data for calculating mean and covmat')
parser.add_argument('--data_file_noisy', default='hw3-data/dataV.csv', type=str, help='data for resconstruction')
parser.add_argument('--save_reconstruct', default=False, type=bool)
args = parser.parse_args()

data = pd.read_csv(args.data_file)
data = pd.DataFrame.as_matrix(data)
data_noisy = pd.read_csv(args.data_file_noisy)
data_noisy = pd.DataFrame.as_matrix(data_noisy)


num_principal_components = list(range(5))
# num_principal_components = [2] # For reconstructing dataII with two principle components
for n_components in num_principal_components:
    pca_noiseless = PCA(n_components=n_components)
    pca_noiseless.fit(data)
    recon_noiseless = pca_noiseless.inverse_transform(pca_noiseless.transform(data_noisy))
    mse_noiseless = mean_squared_error(data, recon_noiseless)

    pca_noise = PCA(n_components=n_components)
    pca_noise.fit(data_noisy)
    recon_noise = pca_noise.inverse_transform(pca_noise.transform(data_noisy))
    if (args.save_reconstruct):
        pd.DataFrame(recon_noise).to_csv('yuchecw2-recon.csv')
    mse_noise = mean_squared_error(data, recon_noise)
    print('n_components: {}, mse_noiseless: {}, mse_noise: {}'.format(n_components, mse_noiseless, mse_noise))
<file_sep># Read CSV into R
# Part 1A
data = read.csv(file="Desktop/pima-indians-diabetes.csv")
accuracy_vec = seq(1, 10, by=1)
train_ratio = 0.8
attribute = 1
for (iter in seq(1, 10, by=1)) {
  # Test-train spilts
  train_test_list = test_train_split(data, train_ratio, iter)
  train_data = train_test_list[[1]]
  test_data = train_test_list[[2]]
  # Calculate prior probabilities
  outcome = train_data[,9]
  n_positive = length(outcome[outcome==1])
  n_negative = length(outcome[outcome==0])
  P_positive = n_positive/length(outcome)
  P_negative = n_negative/length(outcome)
  # Calculate likelihood
  feature = train_data[,1]
  positive_feature_mean = mean(feature[outcome==1])
  positive_feature_var = var(feature[outcome==1])
  negative_feature_mean = mean(feature[outcome==0])
  negative_feature_var = var(feature[outcome==0])
  confusion_matrix = estimate(test_data, positive_feature_mean, positive_feature_var, 
                              negative_feature_mean, negative_feature_var, P_positive, P_negative)
  n_correct = confusion_matrix[1,1] + confusion_matrix[2,2]
  accuracy = n_correct / nrow(test_data)
  accuracy_vec[iter] = accuracy
}
print(accuracy_vec)

# Evaluate on test data and return the confusion matrix
estimate = function(test_data, positive_mean, positive_var, negative_mean, negative_var, P_positive, P_negative) {
  true_positive = 0
  false_positive = 0
  true_negative = 0
  false_negative = 0
  logP_positive = log(P_positive)
  logP_negative = log(P_negative)
  for(index in seq(1, nrow(test_data), by=1)) {
    positive_log_p = evaluate_log_p(test_data[index,1], positive_mean, positive_var)
    negative_log_p = evaluate_log_p(test_data[index,1], negative_mean, negative_var)
    if (logP_positive + positive_log_p > logP_negative + negative_log_p) {
      if (test_data[index, 9] == 1) {
        true_positive = true_positive + 1
      } else {
        false_positive = false_positive + 1
      }
    } else {
      if (test_data[index, 9] == 0) {
        true_negative = true_negative + 1
      } else {
        false_negative = false_negative + 1
      }
    }
  }
  confusion_matrix = matrix(
    c(true_positive, false_positive, false_negative, true_negative),
    nrow = 2,
    ncol = 2,
    byrow = TRUE
  )
  return(confusion_matrix)
}

# Split the data
test_train_split = function(data, train_ratio, seed) {
  set.seed(seed)
  sample = sample(seq_len(nrow(data)), size=floor(train_ratio*nrow(data)))
  train_data = data[sample,]
  test_data = data[-sample,]
  return(list(train_data, test_data))
}

# Evalute log probability
evaluate_log_p = function(feature, feature_mean, feature_var) {
  log_p = - ((feature - feature_mean)^2)/(2*feature_var) 
  log_p = log_p + log(1/sqrt(2*pi*feature_var))
  return(log_p)
}
<file_sep>import numpy as np
import pandas as pd
import math

# Input: feature matrix of shape (43958, 14)
def preprocess(train_feature):
    workclass = {' Private': 0, ' Self-emp-not-inc': 1, ' Self-emp-inc': 2, ' Federal-gov': 3,
                ' Local-gov': 4, ' State-gov': 5, ' Without-pay': 6, ' Never-worked': 7};
    education = {' Bachelors': 0, ' Some-college': 1, ' 11th': 2, ' HS-grad': 3, ' Prof-school': 4,
            ' Assoc-acdm': 5, ' Assoc-voc': 6, ' 9th': 7, ' 7th-8th': 8, ' 12th': 9, ' Masters': 10,
            ' 1st-4th': 11, ' 10th': 12, ' Doctorate': 13, ' 5th-6th': 14, ' Preschool': 15};
    martial_status = {' Married-civ-spouse': 0, ' Divorced': 1, ' Never-married': 2, ' Separated': 3,
                    ' Widowed': 4, ' Married-spouse-absent': 5, ' Married-AF-spouse': 6};
    occupation = {' Tech-support': 0, ' Craft-repair': 1, ' Other-service': 2, ' Sales': 3, ' Exec-managerial': 4,
            ' Prof-specialty': 5, ' Handlers-cleaners': 6, ' Machine-op-inspct': 7, ' Adm-clerical': 8,
            ' Farming-fishing': 9, ' Transport-moving': 10, ' Priv-house-serv': 11, ' Protective-serv': 12,
            ' Armed-Forces': 13};
    relationship = { ' Wife': 0, ' Own-child': 1, ' Husband': 2, ' Not-in-family': 3, ' Other-relative': 4,
                ' Unmarried': 5};
    race = {' White': 0, ' Asian-Pac-Islander': 1, ' Amer-Indian-Eskimo': 2, ' Other': 3, ' Black': 4};
    sex = {' Female': 0, ' Male': 1};
    native_country ={' United-States': 0, ' Cambodia': 1, ' England': 2, ' Puerto-Rico': 3, ' Canada': 4,
                ' Germany': 5, ' Outlying-US(Guam-USVI-etc)': 6, ' India': 7, ' Japan': 8,
                ' Greece': 9, ' South': 10, ' China': 11, ' Cuba': 12, ' Iran': 13, ' Honduras': 14,
                ' Philippines': 15, ' Italy': 16, ' Poland': 17, ' Jamaica': 18, ' Vietnam': 19,
                ' Mexico': 20, ' Portugal': 21, ' Ireland': 22, ' France': 23, ' Dominican-Republic': 24,
                ' Laos': 25, ' Ecuador': 26, ' Taiwan': 27, ' Haiti': 28, ' Columbia': 29, ' Hungary': 30,
                ' Guatemala': 31, ' Nicaragua': 32, ' Scotland': 33, ' Thailand': 34, ' Yugoslavia': 35,
                ' El-Salvador': 36, ' Trinadad&Tobago': 37, ' Peru': 38, ' Hong': 39, ' Holand-Netherlands': 40};
    missing_workclass = []
    has_workclass = []
    missing_education = []
    has_education = []
    missing_martial_status = []
    has_martial_status = []
    missing_martial_status = []
    has_martial_status = []
    missing_occupation = []
    has_occupation = []
    missing_relationship = []
    has_relationship = []
    missing_race = []
    has_race = []
    missing_sex = []
    has_sex = []
    missing_native_country = []
    has_native_country = []

    for i in range(train_feature.shape[0]):
        # workclass
        if not train_feature[i, 1] == ' ?':
            train_feature[i, 1] = workclass[train_feature[i, 1]]
            has_workclass.append(i)
        else:
            missing_workclass.append(i)
        # education
        if not train_feature[i, 3] == ' ?':
            train_feature[i, 3] = education[train_feature[i, 3]]
            has_education.append(i)
        else:
            missing_education.append(i)
        # martial_status
        if not train_feature[i, 5] == ' ?':
            train_feature[i, 5] = martial_status[train_feature[i, 5]]
            has_martial_status.append(i)
        else:
            missing_martial_status.append(i)
        # occupation
        if not train_feature[i, 6] == ' ?':
            train_feature[i, 6] = occupation[train_feature[i, 6]]
            has_martial_status.append(i)
        else:
            missing_occupation.append(i)
        # relationship
        if not train_feature[i, 7] == ' ?':
            train_feature[i, 7] = relationship[train_feature[i, 7]]
            has_relationship.append(i)
        else:
            missing_relationship.append(i)
        # race
        if not train_feature[i, 8] == ' ?':
            train_feature[i, 8] = race[train_feature[i, 8]]
            has_race.append(i)
        else:
            missing_race.append(i)
        # sex
        if not train_feature[i, 9] == ' ?':
            train_feature[i, 9] = sex[train_feature[i, 9]]
            has_sex.append(i)
        else:
            missing_sex.append(i)
        # native_country
        if not train_feature[i, 13] == ' ?':
            train_feature[i, 13] = native_country[train_feature[i, 13]]
            has_native_country.append(i)
        else:
            missing_native_country.append(i)

    substitue_val = {}
    substitue_val['workclass'] = round(train_feature[:, 1][has_workclass].mean())
    for j in missing_workclass:
        train_feature[j,1] = substitue_val['workclass']
    substitue_val['education'] = round(train_feature[:, 3][has_education].mean())
    for j in missing_education:
        train_feature[j,3] = substitue_val['education']
    substitue_val['martial_status'] = round(train_feature[:, 5][has_martial_status].mean())
    for j in missing_martial_status:
        train_feature[j,5] = substitue_val['martial_statu']
    substitue_val['occupation'] = round(train_feature[:, 6][has_occupation].mean())
    for j in missing_occupation:
        train_feature[j,6] = substitue_val['occupation']
    substitue_val['relationship'] = round(train_feature[:, 7][has_relationship].mean())
    for j in missing_relationship:
        train_feature[j,7] = substitue_val['relationship']
    substitue_val['race'] = round(train_feature[:, 8][has_race].mean())
    for j in missing_race:
        train_feature[j,8] = substitue_val['race']
    substitue_val['sex'] = round(train_feature[:, 9][has_sex].mean())
    for j in missing_sex:
        train_feature[j,9] = substitue_val['sex']
    substitue_val['native_country'] = round(train_feature[:, 13][has_native_country].mean())
    for j in missing_native_country:
        train_feature[j,13] = substitue_val['native_country']

    return train_feature

# Input
# a: numpy array of size (14, 1)
# feature: numpy array of size (14, 1)
# b: double
# Output
# gamma: double
def compute_gamma(a, feature, b):
    gamma = np.matmul(np.transpose(a), feature).item() + b
    return gamma

# Input
# gamma: numpy array of size (14, 1)
def cost_function(gamma, label):
    cost = max(0, 1-label*gamma)
    return cost

# Input
# batch_num: a list consists of indices of training feature, supposed to be length of 1
# train_feature: numpy array of shape (43958, 14)
# train_label: numpy array of shape (43958, )
def hinge_loss(batch_num, train_feature, train_label, a, b):
    loss = 0
    for i in batch_num:
        feature_i = np.expand_dims(train_feature[i,:], axis=1)  # shape=(14,1)
        gamma = compute_gamma(a, feature_i, b)
        loss += cost_function(gamma, train_label[i])
    loss /= len(batch_num)
    return loss

# Input
# lamda: regularization parameter
# a: numpy array of shape (14, 1)
def regularization_loss(lamda, a):
    return 0.5*lamda*(np.matmul(np.transpose(a), a).item())

# This is the loss we're going to minimize using stochastic gradient descent
def total_loss(batch_num, train_feature, train_label, a, b, lamda):
    return hinge_loss(batch_num, train_feature, train_label, a, b) + regularization_loss(lamda, a)

# Output
# u: numpy array of shape (15, 1)
def obtain_u(a, b):
    b_arr = np.expand_dims(np.array([b]), axis=1)
    u = np.concatenate((a, b_arr), axis=0)
    return u

def compute_lr(m, n, epoch):
    return m/(n+epoch)

# Input
# current_a: numpy array of shape (14, 1)
# current_b: double
def update_a(current_a, current_b, train_feature, train_label, batch_num, lr, lamda):
    grad = np.zeros(current_a.shape)
    for i in batch_num:
        if (cost_function(compute_gamma(current_a, train_feature[i,:], current_b), train_label[i]) == 0):
            grad += lamda*current_a
        else:
            grad = grad + (lamda*current_a - np.expand_dims(train_label[i]*train_feature[i, :], axis=1))
    grad /= len(batch_num)
    return current_a - lr*grad

def update_b(current_b, current_a, train_feature, train_label, batch_num, lr, lamda):
    grad = 0
    for i in batch_num:
        if (cost_function(compute_gamma(current_a, train_feature[i, :], current_b), train_label[i]) == 0):
            continue
        else:
            grad += (-train_label[i])
    grad /= len(batch_num)
    return current_b - lr*grad

def evaluate(held_out_feature, held_out_label, a, b, lamda):
    n_correct = 0
    for i in range(held_out_feature.shape[0]):
        gamma_i = compute_gamma(a, held_out_feature[i,:], b)
        if gamma_i*held_out_label[i] > 0:
            n_correct += 1
    accuracy = n_correct / held_out_feature.shape[0]
    mag = math.sqrt(np.matmul(a.T, a).item() + b**2)
    loss = total_loss(range(held_out_feature.shape[0]), held_out_feature, held_out_label, a, b, lamda)
    return {'accuracy': accuracy, 'mag': mag, 'loss': loss}

def predict(test_feature, a, b):
    test_label = np.empty((test_feature.shape[0]), dtype=str)
    for i in range(test_feature.shape[0]):
        gamma = compute_gamma(a, np.expand_dims(test_feature[i,:], axis=1), b)
        if (gamma > 0):
            test_label[i] = ">50K"
        else:
            test_label[i] = "<=50K"
    return pd.DataFrame(data={'Label': test_label}, dtype=str)
<file_sep>''' Autoencoder for MNIST dataset'''
from __future__ import division, print_function, absolute_import

import tensorflow as tf
from tensorflow.examples.tutorials.mnist import input_data
import numpy as np
import matplotlib.pyplot as plt
import argparse

def weight_variable(shape, name):
    return tf.Variable(tf.truncated_normal(shape, stddev=0.1), name=name)

def bias_variable(shape, name):
    return tf.Variable(tf.constant(0.1, shape=shape), name=name)

def encoder(x, in_channels, hid_channels_1, hid_channels_2, name='encoder'):
    with tf.name_scope(name):
        w1 = weight_variable([in_channels, hid_channels_1], name='w1')
        b1 = bias_variable([hid_channels_1], name='b1')
        hid_out_1 = tf.nn.sigmoid(tf.matmul(x, w1) + b1)
        w2 = weight_variable([hid_channels_1, hid_channels_2], name='w2')
        b2 = bias_variable([hid_channels_2], name='b2')
        hid_out_2 = tf.nn.sigmoid(tf.matmul(hid_out_1, w2) + b2)
        return hid_out_2

def decoder(x, hid_channels_2, hid_channels_1, out_channels, name='decoder'):
    with tf.name_scope(name):
        w1 = weight_variable([hid_channels_2, hid_channels_1], name='w1')
        b1 = bias_variable([hid_channels_1], name='b1')
        hid_out_1 = tf.nn.sigmoid(tf.matmul(x, w1) + b1)
        w2 = weight_variable([hid_channels_1, out_channels], name='w2')
        b2 = bias_variable([out_channels], name='b2')
        hid_out_2 = tf.nn.sigmoid(tf.matmul(hid_out_1, w2) + b2)
        return hid_out_2

def train():
    mnist = input_data.read_data_sets('/tmp/data', one_hot=True)
    img_test = mnist.test.images
    label_test = mnist.test.labels
    label_img_dict = {}
    for i in range(label_test.shape[0]):
        num = np.argmax(label_test[i,:])
        if num not in label_img_dict.keys():
            label_img_dict[num] = img_test[i,:]
        else:
            label_img_dict[num] = np.vstack((label_img_dict[num], img_test[i,:]))

    # Network parameters
    in_channels = 784 # = out_channels
    hid_channels_1 = 256
    hid_channels_2 = 32
    # Training parameters
    learning_rate = 0.01
    num_epochs = 100000
    display_epoch = 1000
    batch_size = 1

    x = tf.placeholder(tf.float32, shape=[None, in_channels], name='x') # Input image a.k.a groud truth
    y = tf.placeholder(tf.float32, shape=[None, hid_channels_2], name='y') # Encoder output

    # Auto-encoder Model
    encoder_out = encoder(x, in_channels, hid_channels_1, hid_channels_2) # Latent representation
    decoder_out = decoder(y, hid_channels_2, hid_channels_1, in_channels) # Reconstructed image

    # Calculate loss
    loss = tf.losses.mean_squared_error(x, decoder_out)
    optimizer = tf.train.RMSPropOptimizer(learning_rate).minimize(loss)

    # Initialize the variables
    init = tf.global_variables_initializer()

    with tf.Session() as sess:
        sess.run(init)
        print('Start session!')
        # Training
        loss_list = []
        for epoch in range(num_epochs+1):
            batch = mnist.train.next_batch(batch_size)
            en_out = sess.run(encoder_out, feed_dict={x: batch[0]})
            _, mse_loss = sess.run([optimizer, loss], feed_dict={x: batch[0], y: en_out})
            loss_list.append(mse_loss)
            if epoch % display_epoch == 0:
                print('Epoch %d, loss: %.3f'%(epoch, mse_loss))
        plt.figure()
        plt.plot(list(range(len(loss_list))), loss_list)
        plt.xlabel('Epoch number')
        plt.ylabel('Loss')
        plt.title('Loss v.s. Epoch')
        plt.show()
        # Testing
        for are_same_digits in [True, False]:
            plt.figure(figsize=(20,18))
            for digit in range(10):
                digit1, digit2 = None, None
                if are_same_digits: # Same digits
                    digit_arr = label_img_dict[digit]
                    # Randomly select two same digits
                    same_digits = digit_arr[np.random.choice(digit_arr.shape[0], 2, replace=False), :]
                    digit1 = same_digits[0,:].reshape(1,-1)
                    digit2 = same_digits[1,:].reshape(1,-1)
                else:   # Different digits
                    # Randomly select two different digits
                    permutation_arr = np.random.permutation(10)
                    num1 = permutation_arr[0]
                    num2 = permutation_arr[1]
                    digit1 = label_img_dict[num1][np.random.choice(label_img_dict[num1].shape[0], 1, replace=False), :]
                    digit2 = label_img_dict[num2][np.random.choice(label_img_dict[num2].shape[0], 1, replace=False), :]
                # Calculate corresponding codes
                code1 = sess.run(encoder_out, feed_dict={x: digit1})
                code2 = sess.run(encoder_out, feed_dict={x: digit2})
                reconstructed_img = [None]*9
                reconstructed_img[0], reconstructed_img[8] = digit1, digit2
                # Compute 7 evenly spaced linear interpolates
                for t in range(1, 8):
                    code = ((8-t)/8)*code1 + (t/8)*code2
                    reconstructed_img[t] = sess.run(decoder_out, feed_dict={y: code})
                for idx in range(9):
                    plt.subplot(10,9,9*digit+idx+1)
                    plt.imshow(reconstructed_img[idx].reshape(28,28), cmap='gray')
            if are_same_digits:
                plt.savefig('same.png')
            else:
                plt.savefig('different.png')

train()
<file_sep># Part 1B
# Read CSV into R
data = read.csv(file="Desktop/pima-indians-diabetes.csv")
accuracy_vec = seq(1, 10, by=1)
train_ratio = 0.8
attributes = c(3, 4, 6, 8) #blood pressure, skin thickness, BMI, age
param = matrix(, nrow=8, ncol=4)

for (iter in seq(1, 10, by=1)) {
  # Test-train spilts
  train_test_list = test_train_split(data, train_ratio, iter)
  train_data = train_test_list[[1]]
  test_data = train_test_list[[2]]
  # Calculate prior probabilities
  outcome = train_data[,9]
  n_positive = length(outcome[outcome==1])
  n_negative = length(outcome[outcome==0])
  P_positive = n_positive/length(outcome)
  P_negative = n_negative/length(outcome)
  # Calculate likelihood
  # Calculate normal distribution parameter for each feature
  for (attribute in attributes) {
    feature = train_data[,attribute]
    # Ignore zero entries aka missing values
    # I don't change 0 entries to NA, ignore those entries in this step
    remain_outcome = outcome[feature!=0]
    feature = feature[feature!=0]
    param[attribute, 1] = mean(feature[remain_outcome==1])   # positive_feature_mean
    param[attribute, 2] = var(feature[remain_outcome==1])    # positive_feature_var
    param[attribute, 3] = mean(feature[remain_outcome==0])   # negative_feature_mean
    param[attribute, 4] = var(feature[remain_outcome==0])    # negative_feature_var
  }
  confusion_matrix = estimate(test_data, param, attribute, P_positive, P_negative)
  n_correct = confusion_matrix[1,1] + confusion_matrix[2,2]
  accuracy = n_correct / nrow(test_data)
  accuracy_vec[iter] = accuracy
}
print(accuracy_vec)
print(mean(accuracy_vec))

# Evaluate on test data and return the confusion matrix
estimate = function(test_data, param, attribute, P_positive, P_negative) {
  true_positive = 0
  false_positive = 0
  true_negative = 0
  false_negative = 0
  logP_positive = log(P_positive)
  logP_negative = log(P_negative)
  for(index in seq(1, nrow(test_data), by=1)) {
    positive_log_p = 0
    negative_log_p = 0
    # Calculate log likelihood
    for (attribute in attributes) {
      if (test_data[index, attribute] == 0) {
        # Ignore this attribute
        next
      } else {
        positive_log_p = positive_log_p + evaluate_log_p(test_data[index,attribute], param[attribute, 1], param[attribute, 2])
        negative_log_p = negative_log_p + evaluate_log_p(test_data[index,attribute], param[attribute, 3], param[attribute, 4])
      }
    }
    if (logP_positive + positive_log_p > logP_negative + negative_log_p) {
      if (test_data[index, 9] == 1) {
        true_positive = true_positive + 1
      } else {
        false_positive = false_positive + 1
      }
    } else {
      if (test_data[index, 9] == 0) {
        true_negative = true_negative + 1
      } else {
        false_negative = false_negative + 1
      }
    }
  }
  confusion_matrix = matrix(
    c(true_positive, false_positive, false_negative, true_negative),
    nrow = 2,
    ncol = 2,
    byrow = TRUE
  )
  return(confusion_matrix)
}

# Split the data
test_train_split = function(data, train_ratio, seed) {
  set.seed(seed)
  sample = sample(seq_len(nrow(data)), size=floor(train_ratio*nrow(data)))
  train_data = data[sample,]
  test_data = data[-sample,]
  return(list(train_data, test_data))
}

# Evalute log probability
evaluate_log_p = function(feature, feature_mean, feature_var) {
  log_p = - ((feature - feature_mean)^2)/(2*feature_var) 
  log_p = log_p + log(1/sqrt(2*pi*feature_var))
  return(log_p)
}
<file_sep>library(caret)
library(e1071)
library(klaR)
# Read CSV into R
data = read.csv(file="Desktop/AML/pima-indians-diabetes.csv", header=TRUE, sep=",")
accuracy_vec = seq(1, 10, by=1)
for(seed in seq(1, 10)) {
  # Test train split
  set.seed(seed)
  in_train = createDataPartition(data[,9], p=0.8, list=FALSE)
  train_data = data[in_train, ]
  test_data = data[-in_train, ]
  # Train the model
  feature = train_data[,-9]
  label = as.factor(train_data[,9])
  model = train(feature, label, 'nb', trControl=trainControl(method='cv', number=10))
  prediction = predict(model$finalModel, test_data[,-9])
  label_prediction = prediction$class
  # Calculate confusion matrix
  true_positive = 0
  false_positive = 0
  true_negative = 0
  false_negative = 0
  for(i in seq(1, nrow(test_data))) {
    if (label_prediction[i] == "Yes") {
      if (test_data[i,9] == 1) {
        true_positive = true_positive + 1
      } else {
        false_positive = false_positive + 1
      }
    } else {
      if (test_data[i, 9] == 0) {
        true_negative = true_negative + 1
      } else {
        false_negative = false_negative + 1
      }
    } 
  }
  confusion_matrix = matrix(
    c(true_positive, false_positive, false_negative, true_negative),
    nrow = 2,
    ncol = 2,
    byrow = TRUE
  )
  # print(confusion_matrix)
  n_correct = true_positive + true_negative
  accuracy = n_correct/nrow(test_data)
  accuracy_vec[seed] = accuracy
}
print(accuracy_vec)
print(mean(accuracy_vec))
<file_sep>import numpy as np
import numpy.linalg as la
import pandas as pd
import seaborn as sns
import scipy.stats as stats
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from statsmodels.stats.outliers_influence import OLSInfluence
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Load the data
data = np.loadtxt('data.txt')
df = pd.DataFrame(data)

# Linear regression
X = data[:,:13]
y = data[:,13].reshape(-1,1)
# Method 1
lm = sm.OLS(y, sm.add_constant(X)).fit()
# Method 2
reg = LinearRegression().fit(X, y)
r_squared = reg.score(X,y)

# Generate influence plot
fig, ax = plt.subplots(figsize=(9,6))
fig = sm.graphics.influence_plot(lm, alpha=0.001, ax=ax, criterion='cooks')
fig.show()

####################
# Calculate leverage
####################
beta_hat = np.dot(la.inv(np.dot(X.T, X)), np.dot(X.T, y))
hat_matrix = X@la.inv(np.dot(X.T, X))@(X.T)
# plt.bar(list(range(hat_matrix.shape[0])),np.diag(hat_matrix))
# plt.boxplot(np.diag(hat_matrix))
leverage = np.diag(hat_matrix)
# Remove leverage outliers and record their indices
outliers_num_leverage= []
for i in range(leverage.shape[0]):
    if (leverage[i] > np.percentile(leverage, 75)+1.5*stats.iqr(leverage)):
        outliers_num_leverage.append(i)
leverage_remove_outliers = leverage[leverage <= np.percentile(leverage, 75)+1.5*stats.iqr(leverage)]

# Calculate residuals and mean square error
y_predicted = X@beta_hat # Fitted value
e = y - X@beta_hat
N = y.shape[0]
mean_squared_error = np.dot(e.T, e)/N
# r_squred = np.var(X@beta_hat)/np.var(y)

###########################################
# Calculate the cook distance for each data
###########################################
cook_distance = np.zeros((y.shape[0],1))
for i in range(y.shape[0]):
    X_remove = np.delete(X, i, 0)
    y_remove = np.delete(y, i, 0)
    beta_i = np.dot(la.inv(np.dot(X_remove.T, X_remove)), np.dot(X_remove.T, y_remove))
    y_p = X_remove@la.inv(np.dot(X_remove.T, X_remove))@(X_remove.T)@y_remove
    cook_distance[i] = np.dot((y_p-X_remove@beta_i).T, y_p-X_remove@beta_i)/(N*mean_squared_error)
# plt.bar(list(range(y.shape[0])), cook_distance)
# plt.boxplot(cook_distance)
outliers_num_cook_distance = []
for i in range(cook_distance.shape[0]):
    if (cook_distance[i] > np.percentile(cook_distance, 75)+1.5*stats.iqr(cook_distance)):
        outliers_num_cook_distance.append(i)
cook_distance_remove_outliers = cook_distance[cook_distance <= np.percentile(cook_distance, 75)+1.5*stats.iqr(cook_distance)]

###################################
# Calculate standardized residuals
###################################
s = np.zeros((y.shape[0], 1))
for i in range(y.shape[0]):
    s[i] = e[i]/np.sqrt(mean_squared_error*(1-leverage[i]))
# plt.bar(list(range(y.shape[0])), np.abs(s))
# plt.boxplot(s)
outliers_num_s = []
for i in range(s.shape[0]):
    if (s[i] > np.percentile(s, 75)+1.5*stats.iqr(s)):
        outliers_num_s.append(i)

# Remove outliers
by_eye = True
if by_eye:
    outliers_num = [364, 365, 368, 370, 371, 372, 380, 405, 410, 418]
else:
    outliers_num = np.union1d(outliers_num_leverage, outliers_num_cook_distance)
    outliers_num = np.union1d(outliers_num, outliers_num_s)

lm_remove = sm.OLS(y_remove, sm.add_constant(X_remove)).fit()
X_remove = np.delete(X, outliers_num, 0)
y_remove = np.delete(y, outliers_num, 0)
s_remove = np.delete(s, outliers_num, 0)

# Page1
n_interval = 50
lmbda_list = np.linspace(-2, 2, n_interval)
llf_list = np.zeros(n_interval)
for i in range(lmbda_list.shape[0]):
    llf_list[i] = stats.boxcox_llf(lmbda_list[i], y_remove)
fig, ax = plt.subplots(figsize=(8,4))
plt.plot(lmbda_list, llf_list)
plt.xlabel('Lambda value')
plt.ylabel('Log likelihood')
plt.title('Box-Cox Transformation Curve')
plt.show()

# Page 2
# Generate diagnostic plots before removing outliers
before = True # Change this variable to False to generate diagnostic plot after removing outliers
if before:
    model_residuals = lm.resid
    model_norm_residuals = lm.get_influence().resid_studentized_internal
    model_leverage = lm.get_influence().hat_matrix_diag
    model_cooks = lm.get_influence().cooks_distance[0]

    diagnostic_plot = plt.figure()
    plt.scatter(model_leverage, model_norm_residuals, alpha=0.5)
    sns.regplot(model_leverage, model_norm_residuals, scatter=False, ci=False, lowess=True,
                        line_kws={'color':'red', 'lw':1, 'alpha':0.8}, label='Cook\'s distance')
    diagnostic_plot.axes[0].set_xlim(0, max(model_leverage)+0.01)
    diagnostic_plot.axes[0].set_ylim(-3,5)
    plt.title('Residuals vs Leverage\nBefore removing outliers')
    plt.xlabel('Leverage')
    plt.ylabel('Standardized Residuals')
    plt.legend()

    leverage_top_3 = np.flip(np.argsort(model_cooks), 0)[:3]
    for i in leverage_top_3:
        diagnostic_plot.axes[0].annotate(i, xy=(model_leverage[i], model_norm_residuals[i]))
    plt.show()

# Boxcox transformation
y_boxcox, maxlog = stats.boxcox(y_remove)

# Influence plots after removing outliers and boxcox transformation
lm_remove = sm.OLS(y_boxcox, sm.add_constant(X_remove)).fit()
fig, ax = plt.subplots(figsize=(9,6))
fig = sm.graphics.influence_plot(lm_remove, alpha=0.001, ax=ax, criterion='cooks')
fig.show()

# Page 3
# Standardized residuals vs Fitted values without any transformation
sns.residplot(y_predicted, s)
plt.xlabel('Fitted values')
plt.ylabel('Standardized residuals')
plt.title('Without any transformation\nResiduals vs Fitted')
plt.show()
# Standardized residuals vs Fitted values removing all outliers and boxcox transformation
beta_hat_prime = np.dot(la.inv(np.dot(X_remove.T, X_remove)), np.dot(X_remove.T, y_remove))
y_predicted_prime = X_remove@beta_hat_prime # Fitted house price
sns.residplot(y_boxcox, np.delete(s, outliers_num))
plt.xlabel('Fitted values')
plt.ylabel('Standardized residuals')
plt.title('Removing all outliers and apply boxcox transformation\nResiduals vs Fitted')
plt.show()

# Page 4
# Fitted house price vs True house price
plt.scatter(y_predicted_prime, y_remove)
model = LinearRegression().fit(y_predicted_prime, y_remove)
plt.plot(y_predicted_prime, model.predict(y_predicted_prime), 'r')
print(model.score(y_predicted_prime, y_remove))
plt.xlabel('Fitted house price')
plt.ylabel('True hous price')
plt.title('After removing outliers and boxcox transformation')
plt.show()
<file_sep>import numpy as np
import pandas as pd
from sklearn.naive_bayes import BernoulliNB
from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
import matplotlib.pyplot as plt
import random
import math
import argparse
from helper import *

parser = argparse.ArgumentParser()
parser.add_argument('--model_name', default='Gaussian', type=str)
parser.add_argument('--preprocess', default=False, type=bool)
parser.add_argument('--num_epochs', default=10, type=int)
parser.add_argument('--test', default=True, type=bool)
parser.add_argument('--path', default='yuchecw2_1.csv', type=str)
parser.add_argument('--n_estimators', default=30, type=int)
parser.add_argument('--max_depth', default=16, type=int)
args = parser.parse_args()

train_data = pd.read_csv("MNIST/train.csv")
test_data = pd.read_csv("MNIST/test.csv")
train_data = pd.DataFrame.as_matrix(train_data)
test_data = pd.DataFrame.as_matrix(test_data)
if args.preprocess:
    train_data = preprocess(train_data, False)
if not args.test:
    ratio = 0.2
    num_data = train_data.shape[0]
    num_val = (int)(num_data*ratio)
    accuracy_list = []

    for epoch in range(args.num_epochs):
        random.seed(epoch)
        val_samples = random.sample(range(num_data), num_val)
        mask_array = np.ones((num_data, 1), dtype=bool)
        for i in val_samples:
            mask_array[i] = False

        label = np.empty((num_data-num_val, 1))
        feature = np.empty((num_data-num_val, train_data.shape[1]-1))
        idx = 0
        for i in range(num_data):
            if (mask_array[i]):
                label[idx] = train_data[i,0]
                feature[idx, :] = train_data[i,1:]
                idx = idx + 1

        val_label = train_data[val_samples, 0]
        val_feature = train_data[val_samples, 1:]

        # Create a model and train it
        if args.model_name == 'Gaussian':
            model = GaussianNB()
        elif args.model_name == 'Bernoulli':
            # Thresholding
            feature[feature>127] = 255
            feature[feature<=127] = 0
            val_feature[val_feature>127] = 255
            val_feature[val_feature<=127] = 0
            model = BernoulliNB()
        elif args.model_name == 'Random_Forest':
            model = RandomForestClassifier(n_estimators=args.n_estimators, max_depth=args.max_depth)

        label = np.ravel(label)
        fit_model = model.fit(feature, label)
        label_prediction = fit_model.predict(val_feature)
        num_correct = (val_label == label_prediction).sum()
        accuracy = num_correct/num_val
        accuracy_list.append(accuracy)
    print('Training and cross validation: uses %s model, preprocessed: %d, accuracy: %f'%(args.model_name, args.preprocess, sum(accuracy_list)/len(accuracy_list)))
else:
    label = train_data[:,0]
    feature = train_data[:,1:]
    # Create a model and train it
    if args.model_name == 'Gaussian':
        model = GaussianNB()
    elif args.model_name == 'Bernoulli':
        feature[feature>127] = 255
        feature[feature<=127] = 0
        model = BernoulliNB()
    elif args.model_name == 'Random_Forest':
        model = RandomForestClassifier(n_estimators=args.n_estimators, max_depth=args.max_depth)

    label = np.ravel(label)
    fit_model = model.fit(feature, label)

    if args.preprocess:
        test_data = preprocess(test_data, True)
    test_prediction = fit_model.predict(test_data)
    df = pd.DataFrame(test_prediction)
    df_test_data = pd.DataFrame(test_data)
    df.to_csv(args.path)

    # Calculate class means
    if args.preprocess:
        original_size = 20 # box_size
    else:
        original_size = 28

    num_classes = 10
    mean = [None]*num_classes
    for label in range(num_classes):
        mean[label] = df_test_data[df[0]==label].mean()

    mean_matrix = []
    for label in range(num_classes):
        a = np.empty((original_size, original_size))
        for j in range(original_size):
            for i in range(original_size):
                if (mean[label][original_size*j+i] < 127):
                    a[j, i] = 255
                else:
                    a[j, i] = 0
        mean_matrix.append(a)

    plt.figure()
    for label in range(10):
        plt.subplot(1, 10, label+1)
        plt.imshow(mean_matrix[label], cmap='Greys')
        plt.axis('off')
    plt.show()
