cs231n Assignment#1 svm

Posted on 2017-03-03 | | Visitors

Multiclass Support Vector Machine exercise

Complete and hand in this completed worksheet (including its outputs and any supporting code outside of the worksheet) with your assignment submission. For more details see the assignments page on the course website.

In this exercise you will:

implement a fully-vectorized loss function for the SVM
implement the fully-vectorized expression for its analytic gradient
check your implementation using numerical gradient
use a validation set to tune the learning rate and regularization strength
optimize the loss function with SGD
visualize the final learned weights

# Run some setup code for this notebook.
import random
import numpy as np
from cs231n.data_utils import load_CIFAR10
import matplotlib.pyplot as plt
# This is a bit of magic to make matplotlib figures appear inline in the
# notebook rather than in a new window.
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 8.0) # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'
# Some more magic so that the notebook will reload external python modules;
# see http://stackoverflow.com/questions/1907993/autoreload-of-modules-in-ipython
%load_ext autoreload
%autoreload 2

CIFAR-10 Data Loading and Preprocessing

# Load the raw CIFAR-10 data.
cifar10_dir = 'cs231n/datasets/cifar-10-batches-py'
X_train, y_train, X_test, y_test = load_CIFAR10(cifar10_dir)
# As a sanity check, we print out the size of the training and test data.
print 'Training data shape: ', X_train.shape
print 'Training labels shape: ', y_train.shape
print 'Test data shape: ', X_test.shape
print 'Test labels shape: ', y_test.shape

Training data shape:  (50000L, 32L, 32L, 3L)
Training labels shape:  (50000L,)
Test data shape:  (10000L, 32L, 32L, 3L)
Test labels shape:  (10000L,)

# Visualize some examples from the dataset.
# We show a few examples of training images from each class.
classes = ['plane', 'car', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship', 'truck']
num_classes = len(classes)
samples_per_class = 7
for y, cls in enumerate(classes):
    idxs = np.flatnonzero(y_train == y)
    idxs = np.random.choice(idxs, samples_per_class, replace=False)
    for i, idx in enumerate(idxs):
        plt_idx = i * num_classes + y + 1
        plt.subplot(samples_per_class, num_classes, plt_idx)
        plt.imshow(X_train[idx].astype('uint8'))
        plt.axis('off')
        if i == 0:
            plt.title(cls)
plt.show()

1	?np.random.choice

# Split the data into train, val, and test sets. In addition we will
# create a small development set as a subset of the training data;
# we can use this for development so our code runs faster.
num_training = 49000
num_validation = 1000
num_test = 1000
num_dev = 500
# Our validation set will be num_validation points from the original
# training set.
mask = range(num_training, num_training + num_validation)
X_val = X_train[mask]
y_val = y_train[mask]
# Our training set will be the first num_train points from the original
# training set.
mask = range(num_training)
X_train = X_train[mask]
y_train = y_train[mask]
# We will also make a development set, which is a small subset of
# the training set.
mask = np.random.choice(num_training, num_dev, replace=False)
X_dev = X_train[mask]
y_dev = y_train[mask]
# We use the first num_test points of the original test set as our
# test set.
mask = range(num_test)
X_test = X_test[mask]
y_test = y_test[mask]
print 'Train data shape: ', X_train.shape
print 'Train labels shape: ', y_train.shape
print 'Validation data shape: ', X_val.shape
print 'Validation labels shape: ', y_val.shape
print 'Test data shape: ', X_test.shape
print 'Test labels shape: ', y_test.shape

Train data shape:  (49000L, 32L, 32L, 3L)
Train labels shape:  (49000L,)
Validation data shape:  (1000L, 32L, 32L, 3L)
Validation labels shape:  (1000L,)
Test data shape:  (1000L, 32L, 32L, 3L)
Test labels shape:  (1000L,)

# Preprocessing: reshape the image data into rows
X_train = np.reshape(X_train, (X_train.shape[0], -1))
X_val = np.reshape(X_val, (X_val.shape[0], -1))
X_test = np.reshape(X_test, (X_test.shape[0], -1))
X_dev = np.reshape(X_dev, (X_dev.shape[0], -1))
# As a sanity check, print out the shapes of the data
print 'Training data shape: ', X_train.shape
print 'Validation data shape: ', X_val.shape
print 'Test data shape: ', X_test.shape
print 'dev data shape: ', X_dev.shape

Training data shape:  (49000L, 3072L)
Validation data shape:  (1000L, 3072L)
Test data shape:  (1000L, 3072L)
dev data shape:  (500L, 3072L)

# Preprocessing: subtract the mean image
# first: compute the image mean based on the training data
mean_image = np.mean(X_train, axis=0)
print mean_image[:10] # print a few of the elements
plt.figure(figsize=(4,4))
plt.imshow(mean_image.reshape((32,32,3)).astype('uint8')) # visualize the mean image
plt.show()

[ 130.64189796  135.98173469  132.47391837  130.05569388  135.34804082
  131.75402041  130.96055102  136.14328571  132.47636735  131.48467347]

# second: subtract the mean image from train and test data
X_train -= mean_image
X_val -= mean_image
X_test -= mean_image
X_dev -= mean_image

# third: append (at the last) the bias dimension of ones (i.e. bias trick) so that our SVM
# only has to worry about optimizing a single weight matrix W.
X_train = np.hstack([X_train, np.ones((X_train.shape[0], 1))])
X_val = np.hstack([X_val, np.ones((X_val.shape[0], 1))])
X_test = np.hstack([X_test, np.ones((X_test.shape[0], 1))])
X_dev = np.hstack([X_dev, np.ones((X_dev.shape[0], 1))])
print X_train.shape, X_val.shape, X_test.shape, X_dev.shape

(49000L, 3073L) (1000L, 3073L) (1000L, 3073L) (500L, 3073L)

SVM Classifier

Your code for this section will all be written inside cs231n/classifiers/linear_svm.py.

As you can see, we have prefilled the function compute_loss_naive which uses for loops to evaluate the multiclass SVM loss function.

# Evaluate the naive implementation of the loss we provided for you:
from cs231n.classifiers.linear_svm import svm_loss_naive
import time
# generate a random SVM weight matrix of small numbers
W = np.random.randn(3073, 10) * 0.0001 
loss, grad = svm_loss_naive(W, X_dev, y_dev, 0.00001)
print 'loss: %f' % (loss, )

loss: 8.831645

The grad returned from the function above is right now all zero. Derive and implement the gradient for the SVM cost function and implement it inline inside the function svm_loss_naive. You will find it helpful to interleave your new code inside the existing function.

To check that you have correctly implemented the gradient correctly, you can numerically estimate the gradient of the loss function and compare the numeric estimate to the gradient that you computed. We have provided code that does this for you:

# Once you've implemented the gradient, recompute it with the code below
# and gradient check it with the function we provided for you
# Compute the loss and its gradient at W.
loss, grad = svm_loss_naive(W, X_dev, y_dev, 0.0)
# Numerically compute the gradient along several randomly chosen dimensions, and
# compare them with your analytically computed gradient. The numbers should match
# almost exactly along all dimensions.
from cs231n.gradient_check import grad_check_sparse
f = lambda w: svm_loss_naive(w, X_dev, y_dev, 0.0)[0]
grad_numerical = grad_check_sparse(f, W, grad)
# do the gradient check once again with regularization turned on
# you didn't forget the regularization gradient did you?
loss, grad = svm_loss_naive(W, X_dev, y_dev, 1e2)
f = lambda w: svm_loss_naive(w, X_dev, y_dev, 1e2)[0]
grad_numerical = grad_check_sparse(f, W, grad)

numerical: -13.865929 analytic: -13.865929, relative error: 1.283977e-12
numerical: 7.842142 analytic: 7.735021, relative error: 6.876784e-03
numerical: 3.464393 analytic: 3.464393, relative error: 9.040092e-11
numerical: -23.034911 analytic: -23.034911, relative error: 6.876266e-12
numerical: -0.185311 analytic: -0.185311, relative error: 2.538774e-10
numerical: 25.825504 analytic: 25.825504, relative error: 1.336035e-11
numerical: 4.457836 analytic: 4.457836, relative error: 1.015819e-10
numerical: 3.184691 analytic: 3.184691, relative error: 8.849109e-11
numerical: 10.428446 analytic: 10.374317, relative error: 2.601982e-03
numerical: 12.479957 analytic: 12.479957, relative error: 6.825191e-12
numerical: 12.237949 analytic: 12.326308, relative error: 3.597051e-03
numerical: 4.377103 analytic: 4.377103, relative error: 3.904758e-11
numerical: -1.951930 analytic: -1.951930, relative error: 1.432276e-10
numerical: 33.752503 analytic: 33.752503, relative error: 4.254520e-12
numerical: 11.367149 analytic: 11.367149, relative error: 1.682727e-11
numerical: 16.461879 analytic: 16.461879, relative error: 4.766805e-12
numerical: 3.814562 analytic: 3.814562, relative error: 1.087469e-10
numerical: 13.931226 analytic: 13.931226, relative error: 9.578349e-12
numerical: -27.291095 analytic: -27.395406, relative error: 1.907445e-03
numerical: -7.610407 analytic: -7.610407, relative error: 1.015282e-12

Inline Question 1:

It is possible that once in a while a dimension in the gradcheck will not match exactly. What could such a discrepancy be caused by? Is it a reason for concern? What is a simple example in one dimension where a gradient check could fail? Hint: the SVM loss function is not strictly speaking differentiable

Your Answer: Maybe the SVM loss function is not differentiable on that dimension

?np.max

1	np.sum(np.maximum(0, X_dev.dot(W) - X_dev.dot(W)[np.arange(len(y_dev)), [y_dev]].T + 1))

4915.822409730994

# Next implement the function svm_loss_vectorized; for now only compute the loss;
# we will implement the gradient in a moment.
tic = time.time()
loss_naive, grad_naive = svm_loss_naive(W, X_dev, y_dev, 0.00001)
toc = time.time()
print 'Naive loss: %e computed in %fs' % (loss_naive, toc - tic)
from cs231n.classifiers.linear_svm import svm_loss_vectorized
tic = time.time()
loss_vectorized, _ = svm_loss_vectorized(W, X_dev, y_dev, 0.00001)
toc = time.time()
print 'Vectorized loss: %e computed in %fs' % (loss_vectorized, toc - tic)
# The losses should match but your vectorized implementation should be much faster.
print 'difference: %f' % (loss_naive - loss_vectorized)

Naive loss: 8.831645e+00 computed in 0.071000s
Vectorized loss: 8.831645e+00 computed in 0.000000s
difference: 0.000000

# Complete the implementation of svm_loss_vectorized, and compute the gradient
# of the loss function in a vectorized way.
# The naive implementation and the vectorized implementation should match, but
# the vectorized version should still be much faster.
tic = time.time()
_, grad_naive = svm_loss_naive(W, X_dev, y_dev, 0.00001)
toc = time.time()
print 'Naive loss and gradient: computed in %fs' % (toc - tic)
tic = time.time()
_, grad_vectorized = svm_loss_vectorized(W, X_dev, y_dev, 0.00001)
toc = time.time()
print 'Vectorized loss and gradient: computed in %fs' % (toc - tic)
# The loss is a single number, so it is easy to compare the values computed
# by the two implementations. The gradient on the other hand is a matrix, so
# we use the Frobenius norm to compare them.
difference = np.linalg.norm(grad_naive - grad_vectorized, ord='fro')
print 'difference: %f' % difference

Naive loss and gradient: computed in 0.084000s
Vectorized loss and gradient: computed in 0.005000s
difference: 0.000000

Stochastic Gradient Descent

We now have vectorized and efficient expressions for the loss, the gradient and our gradient matches the numerical gradient. We are therefore ready to do SGD to minimize the loss.

# In the file linear_classifier.py, implement SGD in the function
# LinearClassifier.train() and then run it with the code below.
from cs231n.classifiers import LinearSVM
svm = LinearSVM()
tic = time.time()
loss_hist = svm.train(X_train, y_train, learning_rate=1e-7, reg=5e4,
                      num_iters=1500, verbose=True)
toc = time.time()
print 'That took %fs' % (toc - tic)

iteration 0 / 1500: loss 791.772037
iteration 100 / 1500: loss 286.021346
iteration 200 / 1500: loss 107.673095
iteration 300 / 1500: loss 41.812791
iteration 400 / 1500: loss 18.665578
iteration 500 / 1500: loss 10.614984
iteration 600 / 1500: loss 6.664814
iteration 700 / 1500: loss 6.509693
iteration 800 / 1500: loss 5.792204
iteration 900 / 1500: loss 4.986855
iteration 1000 / 1500: loss 5.914691
iteration 1100 / 1500: loss 5.058078
iteration 1200 / 1500: loss 5.491475
iteration 1300 / 1500: loss 5.609450
iteration 1400 / 1500: loss 5.376595
That took 5.454000s

# A useful debugging strategy is to plot the loss as a function of
# iteration number:
plt.plot(loss_hist)
plt.xlabel('Iteration number')
plt.ylabel('Loss value')
plt.show()

# Write the LinearSVM.predict function and evaluate the performance on both the
# training and validation set
y_train_pred = svm.predict(X_train)
print 'training accuracy: %f' % (np.mean(y_train == y_train_pred), )
y_val_pred = svm.predict(X_val)
print 'validation accuracy: %f' % (np.mean(y_val == y_val_pred), )

training accuracy: 0.364980
validation accuracy: 0.378000

# Use the validation set to tune hyperparameters (regularization strength and
# learning rate). You should experiment with different ranges for the learning
# rates and regularization strengths; if you are careful you should be able to
# get a classification accuracy of about 0.4 on the validation set.
learning_rates = [1e-8, 1e-7, 2e-7]
regularization_strengths = [1e4, 2e4, 3e4, 4e4, 5e4, 6e4, 7e4, 8e4, 1e5]
# results is dictionary mapping tuples of the form
# (learning_rate, regularization_strength) to tuples of the form
# (training_accuracy, validation_accuracy). The accuracy is simply the fraction
# of data points that are correctly classified.
results = {}
best_val = -1   # The highest validation accuracy that we have seen so far.
best_svm = None # The LinearSVM object that achieved the highest validation rate.
################################################################################
# TODO:                                                                        #
# Write code that chooses the best hyperparameters by tuning on the validation #
# set. For each combination of hyperparameters, train a linear SVM on the      #
# training set, compute its accuracy on the training and validation sets, and  #
# store these numbers in the results dictionary. In addition, store the best   #
# validation accuracy in best_val and the LinearSVM object that achieves this  #
# accuracy in best_svm.                                                        #
#                                                                              #
# Hint: You should use a small value for num_iters as you develop your         #
# validation code so that the SVMs don't take much time to train; once you are #
# confident that your validation code works, you should rerun the validation   #
# code with a larger value for num_iters.                                      #
################################################################################
for learning_rate in learning_rates:
    for regularization_strength in regularization_strengths:
        svm = LinearSVM()
        loss_hist = svm.train(
            X_train, y_train, learning_rate, \
            regularization_strength, num_iters=1500, batch_size=200)
        y_train_pred = svm.predict(X_train)
        y_val_pred = svm.predict(X_val)
        training_accuracy = np.mean(y_train == y_train_pred)
        validation_accuracy = np.mean(y_val == y_val_pred)
        results[(learning_rate, regularization_strength)] = \
        (training_accuracy, validation_accuracy)
        
        if validation_accuracy > best_val:
            best_val = validation_accuracy
            best_svm = svm
################################################################################
#                              END OF YOUR CODE                                #
################################################################################
    
# Print out results.
for lr, reg in sorted(results):
    train_accuracy, val_accuracy = results[(lr, reg)]
    print 'lr %e reg %e train accuracy: %f val accuracy: %f' % (
                lr, reg, train_accuracy, val_accuracy)
    
print 'best validation accuracy achieved during cross-validation: %f' % best_val

lr 1.000000e-08 reg 1.000000e+04 train accuracy: 0.221898 val accuracy: 0.247000
lr 1.000000e-08 reg 2.000000e+04 train accuracy: 0.233653 val accuracy: 0.258000
lr 1.000000e-08 reg 3.000000e+04 train accuracy: 0.234694 val accuracy: 0.225000
lr 1.000000e-08 reg 4.000000e+04 train accuracy: 0.255959 val accuracy: 0.249000
lr 1.000000e-08 reg 5.000000e+04 train accuracy: 0.259755 val accuracy: 0.273000
lr 1.000000e-08 reg 6.000000e+04 train accuracy: 0.267408 val accuracy: 0.269000
lr 1.000000e-08 reg 7.000000e+04 train accuracy: 0.269102 val accuracy: 0.287000
lr 1.000000e-08 reg 8.000000e+04 train accuracy: 0.277102 val accuracy: 0.285000
lr 1.000000e-08 reg 1.000000e+05 train accuracy: 0.295306 val accuracy: 0.301000
lr 1.000000e-07 reg 1.000000e+04 train accuracy: 0.369388 val accuracy: 0.374000
lr 1.000000e-07 reg 2.000000e+04 train accuracy: 0.380265 val accuracy: 0.390000
lr 1.000000e-07 reg 3.000000e+04 train accuracy: 0.375490 val accuracy: 0.378000
lr 1.000000e-07 reg 4.000000e+04 train accuracy: 0.375633 val accuracy: 0.385000
lr 1.000000e-07 reg 5.000000e+04 train accuracy: 0.369694 val accuracy: 0.375000
lr 1.000000e-07 reg 6.000000e+04 train accuracy: 0.372469 val accuracy: 0.383000
lr 1.000000e-07 reg 7.000000e+04 train accuracy: 0.356000 val accuracy: 0.370000
lr 1.000000e-07 reg 8.000000e+04 train accuracy: 0.352816 val accuracy: 0.355000
lr 1.000000e-07 reg 1.000000e+05 train accuracy: 0.356796 val accuracy: 0.377000
lr 2.000000e-07 reg 1.000000e+04 train accuracy: 0.393510 val accuracy: 0.395000
lr 2.000000e-07 reg 2.000000e+04 train accuracy: 0.377020 val accuracy: 0.382000
lr 2.000000e-07 reg 3.000000e+04 train accuracy: 0.363857 val accuracy: 0.373000
lr 2.000000e-07 reg 4.000000e+04 train accuracy: 0.368714 val accuracy: 0.372000
lr 2.000000e-07 reg 5.000000e+04 train accuracy: 0.361531 val accuracy: 0.364000
lr 2.000000e-07 reg 6.000000e+04 train accuracy: 0.354714 val accuracy: 0.368000
lr 2.000000e-07 reg 7.000000e+04 train accuracy: 0.348306 val accuracy: 0.365000
lr 2.000000e-07 reg 8.000000e+04 train accuracy: 0.358082 val accuracy: 0.378000
lr 2.000000e-07 reg 1.000000e+05 train accuracy: 0.347898 val accuracy: 0.358000
best validation accuracy achieved during cross-validation: 0.395000

# Visualize the cross-validation results
import math
x_scatter = [math.log10(x[0]) for x in results]
y_scatter = [math.log10(x[1]) for x in results]
# plot training accuracy
marker_size = 100
colors = [results[x][0] for x in results]
plt.subplot(3, 1, 1)
plt.scatter(x_scatter, y_scatter, marker_size, c=colors)
plt.colorbar()
plt.xlabel('log learning rate')
plt.ylabel('log regularization strength')
plt.title('CIFAR-10 training accuracy')
# plot validation accuracy
colors = [results[x][1] for x in results] # default size of markers is 20
plt.subplot(3, 1, 3)
plt.scatter(x_scatter, y_scatter, marker_size, c=colors)
plt.colorbar()
plt.xlabel('log learning rate')
plt.ylabel('log regularization strength')
plt.title('CIFAR-10 validation accuracy')
plt.show()

# Evaluate the best svm on test set
y_test_pred = best_svm.predict(X_test)
test_accuracy = np.mean(y_test == y_test_pred)
print 'linear SVM on raw pixels final test set accuracy: %f' % test_accuracy

linear SVM on raw pixels final test set accuracy: 0.383000

1 2	x = np.array([[[0], [1], [2]]]) np.squeeze(x)

array([0, 1, 2])

# Visualize the learned weights for each class.
# Depending on your choice of learning rate and regularization strength, these may
# or may not be nice to look at.
w = best_svm.W[:-1,:] # strip out the bias
w = w.reshape(32, 32, 3, 10)
w_min, w_max = np.min(w), np.max(w)
classes = ['plane', 'car', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship', 'truck']
for i in xrange(10):
  plt.subplot(2, 5, i + 1)
    
  # Rescale the weights to be between 0 and 255
  wimg = 255.0 * (w[:, :, :, i].squeeze() - w_min) / (w_max - w_min)
  plt.imshow(wimg.astype('uint8'))
  plt.axis('off')
  plt.title(classes[i])

Inline question 2:

Describe what your visualized SVM weights look like, and offer a brief explanation for why they look they way that they do.

Your answer: fill this in

Codes

linear_svm.py

import numpy as np
from random import shuffle
def svm_loss_naive(W, X, y, reg):
    """
    Structured SVM loss function, naive implementation (with loops).
    Inputs have dimension D, there are C classes, and we operate on minibatches
    of N examples.
    Inputs:
    - W: A numpy array of shape (D, C) containing weights.
    - X: A numpy array of shape (N, D) containing a minibatch of data.
    - y: A numpy array of shape (N,) containing training labels; y[i] = c means
      that X[i] has label c, where 0 <= c < C.
    - reg: (float) regularization strength
    Returns a tuple of:
    - loss as single float
    - gradient with respect to weights W; an array of same shape as W
    """
    dW = np.zeros(W.shape)  # initialize the gradient as zero
    # compute the loss and the gradient
    num_classes = W.shape[1]
    num_train = X.shape[0]
    loss = 0.0
    for i in xrange(num_train):
        scores = X[i].dot(W)
        correct_class_score = scores[y[i]]
        for j in xrange(num_classes):
            if j == y[i]:
                continue
            margin = scores[j] - correct_class_score + 1  # note delta = 1
            if margin > 0:
                dW[:, y[i]] -= X[i, :]
                dW[:, j] += X[i, :]
                loss += margin
    # Right now the loss is a sum over all training examples, but we want it
    # to be an average instead so we divide by num_train.
    loss /= num_train
    dW /= num_train
    # Add regularization to the loss.
    loss += 0.5 * reg * np.sum(W * W)
    dW += reg * W
    #############################################################################
    # TODO:                                                                     #
    # Compute the gradient of the loss function and store it dW.                #
    # Rather that first computing the loss and then computing the derivative,     #
    # it may be simpler to compute the derivative at the same time that the     #
    # loss is being computed. As a result you may need to modify some of the    #
    # code above to compute the gradient.                                       #
    #############################################################################
    return loss, dW
def svm_loss_vectorized(W, X, y, reg):
    """
    Structured SVM loss function, vectorized implementation.
    Inputs and outputs are the same as svm_loss_naive.
    """
    loss = 0.0
    dW = np.zeros(W.shape)  # initialize the gradient as zero
    #############################################################################
    # TODO:                                                                     #
    # Implement a vectorized version of the structured SVM loss, storing the    #
    # result in loss.                                                           #
    #############################################################################
    num_train = X.shape[0]
    delta = 1.0
    scores = X.dot(W)
    correct_class_score = scores[np.arange(num_train), y]
    margins = np.maximum(
        0, scores - correct_class_score[:, np.newaxis] + delta)
    margins[np.arange(num_train), y] = 0
    loss = np.sum(margins)
    loss /= num_train
    loss += 0.5 * reg * np.sum(W.T.dot(W))
    #############################################################################
    #                             END OF YOUR CODE                              #
    #############################################################################
    #############################################################################
    # TODO:                                                                     #
    # Implement a vectorized version of the gradient for the structured SVM     #
    # loss, storing the result in dW.                                           #
    #                                                                           #
    # Hint: Instead of computing the gradient from scratch, it may be easier    #
    # to reuse some of the intermediate values that you used to compute the     #
    # loss.                                                                     #
    #############################################################################
    X_mask = np.zeros(margins.shape)
    X_mask[margins > 0] = 1
    count = np.sum(X_mask, axis=1)
    X_mask[np.arange(num_train), y] = -count
    dW = X.T.dot(X_mask)
    dW /= num_train
    dW += np.multiply(W, reg)
    #############################################################################
    #                             END OF YOUR CODE                              #
    #############################################################################
    return loss, dW

linear_classifier

import numpy as np
from cs231n.classifiers.linear_svm import *
from cs231n.classifiers.softmax import *
class LinearClassifier(object):
    def __init__(self):
        self.W = None
    def train(self, X, y, learning_rate=1e-3, reg=1e-5, num_iters=100,
              batch_size=200, verbose=False):
        """
        Train this linear classifier using stochastic gradient descent.
        Inputs:
        - X: A numpy array of shape (N, D) containing training data; there are N
          training samples each of dimension D.
        - y: A numpy array of shape (N,) containing training labels; y[i] = c
          means that X[i] has label 0 <= c < C for C classes.
        - learning_rate: (float) learning rate for optimization.
        - reg: (float) regularization strength.
        - num_iters: (integer) number of steps to take when optimizing
        - batch_size: (integer) number of training examples to use at each step.
        - verbose: (boolean) If true, print progress during optimization.
        Outputs:
        A list containing the value of the loss function at each training iteration.
        """
        num_train, dim = X.shape
        num_classes = np.max(y) + 1  # assume y takes values 0...K-1 where K is number of classes
        if self.W is None:
            # lazily initialize W
            self.W = 0.001 * np.random.randn(dim, num_classes)
        # Run stochastic gradient descent to optimize W
        loss_history = []
        for it in xrange(num_iters):
            X_batch = None
            y_batch = None
            #########################################################################
            # TODO:                                                                 #
            # Sample batch_size elements from the training data and their           #
            # corresponding labels to use in this round of gradient descent.        #
            # Store the data in X_batch and their corresponding labels in           #
            # y_batch; after sampling X_batch should have shape (dim, batch_size)   #
            # and y_batch should have shape (batch_size,)                           #
            #                                                                       #
            # Hint: Use np.random.choice to generate indices. Sampling with         #
            # replacement is faster than sampling without replacement.              #
            #########################################################################
            mask = np.random.choice(num_train, batch_size, replace=True)
            X_batch = X[mask]
            y_batch = y[mask]
            #########################################################################
            #                       END OF YOUR CODE                                #
            #########################################################################
            # evaluate loss and gradient
            loss, grad = self.loss(X_batch, y_batch, reg)
            loss_history.append(loss)
            # perform parameter update
            #########################################################################
            # TODO:                                                                 #
            # Update the weights using the gradient and the learning rate.          #
            #########################################################################
            self.W = self.W - learning_rate * grad
            #########################################################################
            #                       END OF YOUR CODE                                #
            #########################################################################
            if verbose and it % 100 == 0:
                print 'iteration %d / %d: loss %f' % (it, num_iters, loss)
        return loss_history
    def predict(self, X):
        """
        Use the trained weights of this linear classifier to predict labels for
        data points.
        Inputs:
        - X: D x N array of training data. Each column is a D-dimensional point.
        Returns:
        - y_pred: Predicted labels for the data in X. y_pred is a 1-dimensional
          array of length N, and each element is an integer giving the predicted
          class.
        """
        X = X.T
        y_pred = np.zeros(X.shape[1])
        ###########################################################################
        # TODO:                                                                   #
        # Implement this method. Store the predicted labels in y_pred.            #
        ###########################################################################
        scores = X.T.dot(self.W)
        y_pred = np.argsort(scores, axis=1)[:, -1]
        ###########################################################################
        #                           END OF YOUR CODE                              #
        ###########################################################################
        return y_pred
    def loss(self, X_batch, y_batch, reg):
        """
        Compute the loss function and its derivative.
        Subclasses will override this.
        Inputs:
        - X_batch: A numpy array of shape (N, D) containing a minibatch of N
          data points; each point has dimension D.
        - y_batch: A numpy array of shape (N,) containing labels for the minibatch.
        - reg: (float) regularization strength.
        Returns: A tuple containing:
        - loss as a single float
        - gradient with respect to self.W; an array of the same shape as W
        """
        pass
class LinearSVM(LinearClassifier):
    """ A subclass that uses the Multiclass SVM loss function """
    def loss(self, X_batch, y_batch, reg):
        return svm_loss_vectorized(self.W, X_batch, y_batch, reg)
class Softmax(LinearClassifier):
    """ A subclass that uses the Softmax + Cross-entropy loss function """
    def loss(self, X_batch, y_batch, reg):
        return softmax_loss_vectorized(self.W, X_batch, y_batch, reg)

cs231n Assignment#1 kNN

Posted on 2017-03-02 | | Visitors

k-Nearest Neighbor (kNN) exercise

The kNN classifier consists of two stages:

During training, the classifier takes the training data and simply remembers it
During testing, kNN classifies every test image by comparing to all training images and transfering the labels of the k most similar training examples
The value of k is cross-validated

In this exercise you will implement these steps and understand the basic Image Classification pipeline, cross-validation, and gain proficiency in writing efficient, vectorized code.

# Run some setup code for this notebook.
import random
import numpy as np
from cs231n.data_utils import load_CIFAR10
import matplotlib.pyplot as plt
# This is a bit of magic to make matplotlib figures appear inline in the notebook
# rather than in a new window.
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 8.0) # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'
# Some more magic so that the notebook will reload external python modules;
# see http://stackoverflow.com/questions/1907993/autoreload-of-modules-in-ipython
%load_ext autoreload
%autoreload 2

# Load the raw CIFAR-10 data.
cifar10_dir = 'cs231n/datasets/cifar-10-batches-py'
X_train, y_train, X_test, y_test = load_CIFAR10(cifar10_dir)
# As a sanity check, we print out the size of the training and test data.
print 'Training data shape: ', X_train.shape
print 'Training labels shape: ', y_train.shape
print 'Test data shape: ', X_test.shape
print 'Test labels shape: ', y_test.shape

Training data shape:  (50000L, 32L, 32L, 3L)
Training labels shape:  (50000L,)
Test data shape:  (10000L, 32L, 32L, 3L)
Test labels shape:  (10000L,)

# Visualize some examples from the dataset.
# We show a few examples of training images from each class.
classes = ['plane', 'car', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship', 'truck']
num_classes = len(classes)
samples_per_class = 7
for y, cls in enumerate(classes):
    idxs = np.flatnonzero(y_train == y)
    idxs = np.random.choice(idxs, samples_per_class, replace=False)
    for i, idx in enumerate(idxs):
        plt_idx = i * num_classes + y + 1
        plt.subplot(samples_per_class, num_classes, plt_idx)
        plt.imshow(X_train[idx].astype('uint8'))
        plt.axis('off')
        if i == 0:
            plt.title(cls)
plt.show()

# Subsample the data for more efficient code execution in this exercise
num_training = 5000
mask = range(num_training)
X_train = X_train[mask]
y_train = y_train[mask]
num_test = 500
mask = range(num_test)
X_test = X_test[mask]
y_test = y_test[mask]

# Reshape the image data into rows
X_train = np.reshape(X_train, (X_train.shape[0], -1)) # Wow~
X_test = np.reshape(X_test, (X_test.shape[0], -1))
print X_train.shape, X_test.shape

(5000L, 3072L) (500L, 3072L)

from cs231n.classifiers import KNearestNeighbor
# Create a kNN classifier instance. 
# Remember that training a kNN classifier is a noop: 
# the Classifier simply remembers the data and does no further processing 
classifier = KNearestNeighbor()
classifier.train(X_train, y_train)

We would now like to classify the test data with the kNN classifier. Recall that we can break down this process into two steps:

First we must compute the distances between all test examples and all train examples.
Given these distances, for each test example we find the k nearest examples and have them vote for the label

Lets begin with computing the distance matrix between all training and test examples. For example, if there are Ntr training examples and Nte test examples, this stage should result in a Nte x Ntr matrix where each element (i,j) is the distance between the i-th test and j-th train example.

First, open cs231n/classifiers/k_nearest_neighbor.py and implement the function compute_distances_two_loops that uses a (very inefficient) double loop over all pairs of (test, train) examples and computes the distance matrix one element at a time.

# Open cs231n/classifiers/k_nearest_neighbor.py and implement
# compute_distances_two_loops.
# Test your implementation:
dists = classifier.compute_distances_two_loops(X_test)
print dists.shape

(500L, 5000L)

# We can visualize the distance matrix: each row is a single test example and
# its distances to training examples
plt.imshow(dists, interpolation='none')
plt.show()

Inline Question #1: Notice the structured patterns in the distance matrix, where some rows or columns are visible brighter. (Note that with the default color scheme black indicates low distances while white indicates high distances.)

What in the data is the cause behind the distinctly bright rows?
What causes the columns?

Your Answer: Maybe exists noises in test data set and train data set.

# Now implement the function predict_labels and run the code below:
# We use k = 1 (which is Nearest Neighbor).
y_test_pred = classifier.predict_labels(dists, k=1)
# Compute and print the fraction of correctly predicted examples
num_correct = np.sum(y_test_pred == y_test)
accuracy = float(num_correct) / num_test
print 'Got %d / %d correct => accuracy: %f' % (num_correct, num_test, accuracy)

Got 137 / 500 correct => accuracy: 0.274000

You should expect to see approximately 27% accuracy. Now lets try out a larger k, say k = 5:

y_test_pred = classifier.predict_labels(dists, k=5)
num_correct = np.sum(y_test_pred == y_test)
accuracy = float(num_correct) / num_test
print 'Got %d / %d correct => accuracy: %f' % (num_correct, num_test, accuracy)

Got 142 / 500 correct => accuracy: 0.284000

You should expect to see a slightly better performance than with k = 1.

# Now lets speed up distance matrix computation by using partial vectorization
# with one loop. Implement the function compute_distances_one_loop and run the
# code below:
dists_one = classifier.compute_distances_one_loop(X_test)
# To ensure that our vectorized implementation is correct, we make sure that it
# agrees with the naive implementation. There are many ways to decide whether
# two matrices are similar; one of the simplest is the Frobenius norm. In case
# you haven't seen it before, the Frobenius norm of two matrices is the square
# root of the squared sum of differences of all elements; in other words, reshape
# the matrices into vectors and compute the Euclidean distance between them.
difference = np.linalg.norm(dists - dists_one, ord='fro')
print 'Difference was: %f' % (difference, )
if difference < 0.001:
  print 'Good! The distance matrices are the same'
else:
  print 'Uh-oh! The distance matrices are different'

Difference was: 0.000000
Good! The distance matrices are the same

# Now implement the fully vectorized version inside compute_distances_no_loops
# and run the code
dists_two = classifier.compute_distances_no_loops(X_test)
# check that the distance matrix agrees with the one we computed before:
difference = np.linalg.norm(dists - dists_two, ord='fro')
print 'Difference was: %f' % (difference, )
if difference < 0.001:
  print 'Good! The distance matrices are the same'
else:
  print 'Uh-oh! The distance matrices are different'

Difference was: 0.000000
Good! The distance matrices are the same

# Let's compare how fast the implementations are
def time_function(f, *args):
  """
  Call a function f with args and return the time (in seconds) that it took to execute.
  """
  import time
  tic = time.time()
  f(*args)
  toc = time.time()
  return toc - tic
two_loop_time = time_function(classifier.compute_distances_two_loops, X_test)
print 'Two loop version took %f seconds' % two_loop_time
one_loop_time = time_function(classifier.compute_distances_one_loop, X_test)
print 'One loop version took %f seconds' % one_loop_time
no_loop_time = time_function(classifier.compute_distances_no_loops, X_test)
print 'No loop version took %f seconds' % no_loop_time
# you should see significantly faster performance with the fully vectorized implementation

Two loop version took 27.001000 seconds
One loop version took 59.630000 seconds
No loop version took 0.205000 seconds

Cross-validation

We have implemented the k-Nearest Neighbor classifier but we set the value k = 5 arbitrarily. We will now determine the best value of this hyperparameter with cross-validation.

num_folds = 5
k_choices = [1, 3, 5, 8, 10, 12, 15, 20, 50, 100]
X_train_folds = []
y_train_folds = []
################################################################################
# TODO:                                                                        #
# Split up the training data into folds. After splitting, X_train_folds and    #
# y_train_folds should each be lists of length num_folds, where                #
# y_train_folds[i] is the label vector for the points in X_train_folds[i].     #
# Hint: Look up the numpy array_split function.                                #
################################################################################
X_train_folds = np.array_split(X_train, num_folds)
y_train_folds = np.array_split(y_train, num_folds)
################################################################################
#                                 END OF YOUR CODE                             #
################################################################################
# A dictionary holding the accuracies for different values of k that we find
# when running cross-validation. After running cross-validation,
# k_to_accuracies[k] should be a list of length num_folds giving the different
# accuracy values that we found when using that value of k.
k_to_accuracies = {}
################################################################################
# TODO:                                                                        #
# Perform k-fold cross validation to find the best value of k. For each        #
# possible value of k, run the k-nearest-neighbor algorithm num_folds times,   #
# where in each case you use all but one of the folds as training data and the #
# last fold as a validation set. Store the accuracies for all fold and all     #
# values of k in the k_to_accuracies dictionary.                               #
################################################################################
for k in k_choices:
    k_to_accuracies[k] = []
    for fold in xrange(num_folds):
        train_X = np.append(
            X_train_folds[:fold], X_train_folds[fold+1:]).reshape(
            (X_train.shape[0] - X_train.shape[0]/num_folds, -1))
        train_y = np.append(
            y_train_folds[:fold], y_train_folds[fold+1:]).reshape(
            (y_train.shape[0] - y_train.shape[0]/num_folds, -1)).flatten()
        classifier.train(train_X, train_y)
        dists = classifier.compute_distances_no_loops(X_train_folds[fold])
        y_test_pred = classifier.predict_labels(dists, k)
        
        num_correct = np.sum(y_test_pred == y_train_folds[fold])
        accuracy = float(num_correct) / len(y_train_folds[fold])
        k_to_accuracies[k].append(accuracy)
        
################################################################################
#                                 END OF YOUR CODE                             #
################################################################################
# Print out the computed accuracies
for k in sorted(k_to_accuracies):
    for accuracy in k_to_accuracies[k]:
        print 'k = %d, accuracy = %f' % (k, accuracy)

k = 1, accuracy = 0.263000
k = 1, accuracy = 0.257000
k = 1, accuracy = 0.264000
k = 1, accuracy = 0.278000
k = 1, accuracy = 0.266000
k = 3, accuracy = 0.241000
k = 3, accuracy = 0.249000
k = 3, accuracy = 0.243000
k = 3, accuracy = 0.273000
k = 3, accuracy = 0.264000
k = 5, accuracy = 0.258000
k = 5, accuracy = 0.273000
k = 5, accuracy = 0.281000
k = 5, accuracy = 0.290000
k = 5, accuracy = 0.272000
k = 8, accuracy = 0.263000
k = 8, accuracy = 0.288000
k = 8, accuracy = 0.278000
k = 8, accuracy = 0.285000
k = 8, accuracy = 0.277000
k = 10, accuracy = 0.265000
k = 10, accuracy = 0.296000
k = 10, accuracy = 0.278000
k = 10, accuracy = 0.284000
k = 10, accuracy = 0.286000
k = 12, accuracy = 0.260000
k = 12, accuracy = 0.294000
k = 12, accuracy = 0.281000
k = 12, accuracy = 0.282000
k = 12, accuracy = 0.281000
k = 15, accuracy = 0.255000
k = 15, accuracy = 0.290000
k = 15, accuracy = 0.281000
k = 15, accuracy = 0.281000
k = 15, accuracy = 0.276000
k = 20, accuracy = 0.270000
k = 20, accuracy = 0.281000
k = 20, accuracy = 0.280000
k = 20, accuracy = 0.282000
k = 20, accuracy = 0.284000
k = 50, accuracy = 0.271000
k = 50, accuracy = 0.288000
k = 50, accuracy = 0.278000
k = 50, accuracy = 0.269000
k = 50, accuracy = 0.266000
k = 100, accuracy = 0.256000
k = 100, accuracy = 0.270000
k = 100, accuracy = 0.263000
k = 100, accuracy = 0.256000
k = 100, accuracy = 0.263000

# plot the raw observations
for k in k_choices:
  accuracies = k_to_accuracies[k]
  plt.scatter([k] * len(accuracies), accuracies)
# plot the trend line with error bars that correspond to standard deviation
accuracies_mean = np.array([np.mean(v) for k,v in sorted(k_to_accuracies.items())])
accuracies_std = np.array([np.std(v) for k,v in sorted(k_to_accuracies.items())])
plt.errorbar(k_choices, accuracies_mean, yerr=accuracies_std)
plt.title('Cross-validation on k')
plt.xlabel('k')
plt.ylabel('Cross-validation accuracy')
plt.show()

# Based on the cross-validation results above, choose the best value for k,   
# retrain the classifier using all the training data, and test it on the test
# data. You should be able to get above 28% accuracy on the test data.
best_k = 10
classifier = KNearestNeighbor()
classifier.train(X_train, y_train)
y_test_pred = classifier.predict(X_test, k=best_k)
# Compute and display the accuracy
num_correct = np.sum(y_test_pred == y_test)
accuracy = float(num_correct) / num_test
print 'Got %d / %d correct => accuracy: %f' % (num_correct, num_test, accuracy)

Got 139 / 500 correct => accuracy: 0.278000

k_nearest_neighbor.py

import numpy as np
import operator
class KNearestNeighbor(object):
    """ a kNN classifier with L2 distance """
    def __init__(self):
        pass
    def train(self, X, y):
        """
        Train the classifier. For k-nearest neighbors this is just
        memorizing the training data.
        Inputs:
        - X: A numpy array of shape (num_train, D) containing the training data
          consisting of num_train samples each of dimension D.
        - y: A numpy array of shape (N,) containing the training labels, where
             y[i] is the label for X[i].
        """
        self.X_train = X
        self.y_train = y
    def predict(self, X, k=1, num_loops=0):
        """
        Predict labels for test data using this classifier.
        Inputs:
        - X: A numpy array of shape (num_test, D) containing
             test data consisting of num_test samples each of dimension D.
        - k: The number of nearest neighbors that vote for
             the predicted labels.
        - num_loops: Determines which implementation to use
                     to compute distances
          between training points and testing points.
        Returns:
        - y: A numpy array of shape (num_test,)
             containing predicted labels for the
             test data, where y[i] is the predicted
             label for the test point X[i].
        """
        if num_loops == 0:
            dists = self.compute_distances_no_loops(X)
        elif num_loops == 1:
            dists = self.compute_distances_one_loop(X)
        elif num_loops == 2:
            dists = self.compute_distances_two_loops(X)
        else:
            raise ValueError('Invalid value %d for num_loops' % num_loops)
        return self.predict_labels(dists, k=k)
    def compute_distances_two_loops(self, X):
        """
        Compute the distance between each test point in X and each training
        point in self.X_train using a nested loop over both the training data
        and the test data.
        Inputs:
        - X: A numpy array of shape (num_test, D) containing test data.
        Returns:
        - dists: A numpy array of shape (num_test, num_train) where dists[i, j]
          is the Euclidean distance between the ith test point and the jth
          training point.
        """
        num_test = X.shape[0]
        num_train = self.X_train.shape[0]
        dists = np.zeros((num_test, num_train))
        for i in xrange(num_test):
            for j in xrange(num_train):
                ###############################################################
                # TODO:                                                             #
                # Compute the l2 distance between the ith test point and the jth    #
                # training point, and store the result in dists[i, j]. You should   #
                # not use a loop over dimension.                                    #
                ###############################################################
                dists[i, j] = np.sqrt(np.sum((X[i, :] - self.X_train[j, :]) ** 2))
            ###############################################################
            #                       END OF YOUR CODE                            #
            ###############################################################
        return dists
    def compute_distances_one_loop(self, X):
        """
        Compute the distance between each test point in X and each training point
        in self.X_train using a single loop over the test data.
        Input / Output: Same as compute_distances_two_loops
        """
        num_test = X.shape[0]
        num_train = self.X_train.shape[0]
        dists = np.zeros((num_test, num_train))
        for i in xrange(num_test):
            ###################################################################
            # TODO:                                                               #
            # Compute the l2 distance between the ith test point and all training #
            # points, and store the result in dists[i, :].                        #
            ###################################################################
            dists[i, :] = np.sqrt(np.sum(np.square(X[i, :] - self.X_train), axis=1))
            ###################################################################
            #                         END OF YOUR CODE                            #
            ###################################################################
        return dists
    def compute_distances_no_loops(self, X):
        """
        Compute the distance between each test point in X and each training point
        in self.X_train using no explicit loops.
        Input / Output: Same as compute_distances_two_loops
        """
        num_test = X.shape[0]
        num_train = self.X_train.shape[0]
        dists = np.zeros((num_test, num_train))
        #######################################################################
        # TODO:                                                                 #
        # Compute the l2 distance between all test points and all training      #
        # points without using any explicit loops, and store the result in      #
        # dists.                                                                #
        #                                                                       #
        # You should implement this function using only basic array operations; #
        # in particular you should not use functions from scipy.                #
        #                                                                       #
        # HINT: Try to formulate the l2 distance using matrix multiplication    #
        #       and two broadcast sums.                                         #
        #######################################################################
        dists = np.sqrt(np.multiply(np.dot(X, self.X_train.T), -2) +
                        np.sum(self.X_train ** 2, axis=1) +
                        np.sum(X ** 2, axis=1)[:, np.newaxis])
        #######################################################################
        #                         END OF YOUR CODE                              #
        #######################################################################
        return dists
    def predict_labels(self, dists, k=1):
        """
        Given a matrix of distances between test points and training points,
        predict a label for each test point.
        Inputs:
        - dists: A numpy array of shape (num_test, num_train) where dists[i, j]
          gives the distance betwen the ith test point and the jth training point.
        Returns:
        - y: A numpy array of shape (num_test,) containing predicted labels for the
          test data, where y[i] is the predicted label for the test point X[i].
        """
        num_test = dists.shape[0]
        y_pred = np.zeros(num_test)
        for i in xrange(num_test):
            # A list of length k storing the labels of the k nearest neighbors to
            # the ith test point.
            closest_y = []
            ###################################################################
            # TODO:                                                                 #
            # Use the distance matrix to find the k nearest neighbors of the ith    #
            # testing point, and use self.y_train to find the labels of these       #
            # neighbors. Store these labels in closest_y.                           #
            # Hint: Look up the function numpy.argsort.                             #
            ###################################################################
            k_nearest_index = np.argsort(dists[i, :])[:k]
            ###################################################################
            # TODO:                                                                 #
            # Now that you have found the labels of the k nearest neighbors, you    #
            # need to find the most common label in the list closest_y of labels.   #
            # Store this label in y_pred[i]. Break ties by choosing the smaller     #
            # label.                                                                #
            ###################################################################
            closest_y = self.y_train[k_nearest_index]
            labels_counts = {}
            for label in closest_y:
                if label in labels_counts.keys():
                    labels_counts[label] += 1
                else:
                    labels_counts[label] = 0
            sorted_labels_counts = sorted(
                labels_counts.items(), key=operator.itemgetter(1), reverse=True)
            y_pred[i] = sorted_labels_counts[0][0]
            ###################################################################
            #                           END OF YOUR CODE                            #
            ###################################################################
        return y_pred

Python data analysis Learning note Ch06

Posted on 2017-02-28 | | Visitors

Data loading, storage, and file formats

from __future__ import division
from numpy.random import randn
import numpy as np
import os
import sys
import matplotlib.pyplot as plt
np.random.seed(12345)
plt.rc('figure', figsize=(10, 6))
from pandas import Series, DataFrame
import pandas as pd
np.set_printoptions(precision=4)
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

%pwd

'C:\\Users\\Ewan\\Downloads\\pydata-book-master'

Reading and Writing Data in Text Format

1	!more ch06\ex1.csv

a,b,c,d,message
1,2,3,4,hello
5,6,7,8,world
9,10,11,12,foo

1 2	df = pd.read_csv('ch06/ex1.csv') df

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

1	pd.read_table('ch06/ex1.csv', sep=',')

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

1	!more ch06\ex2.csv

1,2,3,4,hello
5,6,7,8,world
9,10,11,12,foo

1 2	pd.read_csv('ch06/ex2.csv', header=None) pd.read_csv('ch06/ex2.csv', names=['a', 'b', 'c', 'd', 'message'])

	0	1	2	3	4
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

1 2	names = ['a', 'b', 'c', 'd', 'message'] pd.read_csv('ch06/ex2.csv', names=names, index_col='message')

	a	b	c	d
message
hello	1	2	3	4
world	5	6	7	8
foo	9	10	11	12

1
2
3

!more ch06\csv_mindex.csv
parsed = pd.read_csv('ch06/csv_mindex.csv', index_col=['key1', 'key2'])
parsed

key1,key2,value1,value2
one,a,1,2
one,b,3,4
one,c,5,6
one,d,7,8
two,a,9,10
two,b,11,12
two,c,13,14
two,d,15,16

		value1	value2
key1	key2
one	a	1	2
	b	3	4
	c	5	6
	d	7	8
two	a	9	10
	b	11	12
	c	13	14
	d	15	16

1	list(open('ch06/ex3.txt'))

['            A         B         C\n',
 'aaa -0.264438 -1.026059 -0.619500\n',
 'bbb  0.927272  0.302904 -0.032399\n',
 'ccc -0.264273 -0.386314 -0.217601\n',
 'ddd -0.871858 -0.348382  1.100491\n']

采用正则表达式作为分隔符

1 2	result = pd.read_table('ch06/ex3.txt', sep='\s+') result

	A	B	C
aaa	-0.264438	-1.026059	-0.619500
bbb	0.927272	0.302904	-0.032399
ccc	-0.264273	-0.386314	-0.217601
ddd	-0.871858	-0.348382	1.100491

1 2	!more ch06\ex4.csv pd.read_csv('ch06/ex4.csv', skiprows=[0, 2, 3])

# hey!
a,b,c,d,message
# just wanted to make things more difficult for you
# who reads CSV files with computers, anyway?
1,2,3,4,hello
5,6,7,8,world
9,10,11,12,foo

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

!more ch06\ex5.csv
result = pd.read_csv('ch06/ex5.csv')
result
pd.isnull(result)

something,a,b,c,d,message
one,1,2,3,4,NA
two,5,6,,8,world
three,9,10,11,12,foo

	something	a	b	c	d	message
0	one	1	2	3.0	4	NaN
1	two	5	6	NaN	8	world
2	three	9	10	11.0	12	foo

	something	a	b	c	d	message
0	False	False	False	False	False	True
1	False	False	False	True	False	False
2	False	False	False	False	False	False

1 2	result = pd.read_csv('ch06/ex5.csv', na_values=['NULL']) result

	something	a	b	c	d	message
0	one	1	2	3.0	4	NaN
1	two	5	6	NaN	8	world
2	three	9	10	11.0	12	foo

1 2	sentinels = {'message': ['foo', 'NA'], 'something': ['two']} pd.read_csv('ch06/ex5.csv', na_values=sentinels)

	something	a	b	c	d	message
0	one	1	2	3.0	4	NaN
1	NaN	5	6	NaN	8	world
2	three	9	10	11.0	12	NaN

逐块读取文本文件

1 2	result = pd.read_csv('ch06/ex6.csv') result

	one	two	three	four	key
0	0.467976	-0.038649	-0.295344	-1.824726	L
1	-0.358893	1.404453	0.704965	-0.200638	B
2	-0.501840	0.659254	-0.421691	-0.057688	G
3	0.204886	1.074134	1.388361	-0.982404	R
4	0.354628	-0.133116	0.283763	-0.837063	Q
5	1.817480	0.742273	0.419395	-2.251035	Q
6	-0.776764	0.935518	-0.332872	-1.875641	U
7	-0.913135	1.530624	-0.572657	0.477252	K
8	0.358480	-0.497572	-0.367016	0.507702	S
9	-1.740877	-1.160417	-1.637830	2.172201	G
10	0.240564	-0.328249	1.252155	1.072796	8
11	0.764018	1.165476	-0.639544	1.495258	R
12	0.571035	-0.310537	0.582437	-0.298765	1
13	2.317658	0.430710	-1.334216	0.199679	P
14	1.547771	-1.119753	-2.277634	0.329586	J
15	-1.310608	0.401719	-1.000987	1.156708	E
16	-0.088496	0.634712	0.153324	0.415335	B
17	-0.018663	-0.247487	-1.446522	0.750938	A
18	-0.070127	-1.579097	0.120892	0.671432	F
19	-0.194678	-0.492039	2.359605	0.319810	H
20	-0.248618	0.868707	-0.492226	-0.717959	W
21	-1.091549	-0.867110	-0.647760	-0.832562	C
22	0.641404	-0.138822	-0.621963	-0.284839	C
23	1.216408	0.992687	0.165162	-0.069619	V
24	-0.564474	0.792832	0.747053	0.571675	I
25	1.759879	-0.515666	-0.230481	1.362317	S
26	0.126266	0.309281	0.382820	-0.239199	L
27	1.334360	-0.100152	-0.840731	-0.643967	6
28	-0.737620	0.278087	-0.053235	-0.950972	J
29	-1.148486	-0.986292	-0.144963	0.124362	Y
…	…	…	…	…	…
9970	0.633495	-0.186524	0.927627	0.143164	4
9971	0.308636	-0.112857	0.762842	-1.072977	1
9972	-1.627051	-0.978151	0.154745	-1.229037	Z
9973	0.314847	0.097989	0.199608	0.955193	P
9974	1.666907	0.992005	0.496128	-0.686391	S
9975	0.010603	0.708540	-1.258711	0.226541	K
9976	0.118693	-0.714455	-0.501342	-0.254764	K
9977	0.302616	-2.011527	-0.628085	0.768827	H
9978	-0.098572	1.769086	-0.215027	-0.053076	A
9979	-0.019058	1.964994	0.738538	-0.883776	F
9980	-0.595349	0.001781	-1.423355	-1.458477	M
9981	1.392170	-1.396560	-1.425306	-0.847535	H
9982	-0.896029	-0.152287	1.924483	0.365184	6
9983	-2.274642	-0.901874	1.500352	0.996541	N
9984	-0.301898	1.019906	1.102160	2.624526	I
9985	-2.548389	-0.585374	1.496201	-0.718815	D
9986	-0.064588	0.759292	-1.568415	-0.420933	E
9987	-0.143365	-1.111760	-1.815581	0.435274	2
9988	-0.070412	-1.055921	0.338017	-0.440763	X
9989	0.649148	0.994273	-1.384227	0.485120	Q
9990	-0.370769	0.404356	-1.051628	-1.050899	8
9991	-0.409980	0.155627	-0.818990	1.277350	W
9992	0.301214	-1.111203	0.668258	0.671922	A
9993	1.821117	0.416445	0.173874	0.505118	X
9994	0.068804	1.322759	0.802346	0.223618	H
9995	2.311896	-0.417070	-1.409599	-0.515821	L
9996	-0.479893	-0.650419	0.745152	-0.646038	E
9997	0.523331	0.787112	0.486066	1.093156	K
9998	-0.362559	0.598894	-1.843201	0.887292	G
9999	-0.096376	-1.012999	-0.657431	-0.573315	0

10000 rows × 5 columns

1	pd.read_csv('ch06/ex6.csv', nrows=5)

	one	two	three	four	key
0	0.467976	-0.038649	-0.295344	-1.824726	L
1	-0.358893	1.404453	0.704965	-0.200638	B
2	-0.501840	0.659254	-0.421691	-0.057688	G
3	0.204886	1.074134	1.388361	-0.982404	R
4	0.354628	-0.133116	0.283763	-0.837063	Q

1 2	chunker = pd.read_csv('ch06/ex6.csv', chunksize=1000) chunker

<pandas.io.parsers.TextFileReader at 0x2035229de80>

chunker = pd.read_csv('ch06/ex6.csv', chunksize=1000)
tot = Series([])
for piece in chunker:
    tot = tot.add(piece['key'].value_counts(), fill_value=0)
tot = tot.sort_values(ascending=False)

tot[:10]

E    368.0
X    364.0
L    346.0
O    343.0
Q    340.0
M    338.0
J    337.0
F    335.0
K    334.0
H    330.0
dtype: float64

将数据写出到文本格式

1 2	data = pd.read_csv('ch06/ex5.csv') data

	something	a	b	c	d	message
0	one	1	2	3.0	4	NaN
1	two	5	6	NaN	8	world
2	three	9	10	11.0	12	foo

1 2	data.to_csv('ch06/out.csv') !more ch06\out.csv

,something,a,b,c,d,message
0,one,1,2,3.0,4,
1,two,5,6,,8,world
2,three,9,10,11.0,12,foo

1	data.to_csv(sys.stdout, sep='\|')

|something|a|b|c|d|message
0|one|1|2|3.0|4|
1|two|5|6||8|world
2|three|9|10|11.0|12|foo

1	data.to_csv(sys.stdout, na_rep='NULL')

,something,a,b,c,d,message
0,one,1,2,3.0,4,NULL
1,two,5,6,NULL,8,world
2,three,9,10,11.0,12,foo

1	data.to_csv(sys.stdout, index=False, header=False)

one,1,2,3.0,4,
two,5,6,,8,world
three,9,10,11.0,12,foo

1	data.to_csv(sys.stdout, index=False, columns=['a', 'b', 'c'])

a,b,c
1,2,3.0
5,6,
9,10,11.0

dates = pd.date_range('1/1/2000', periods=7)
dates
ts = Series(np.arange(7), index=dates)
ts
ts.to_csv('ch06/tseries.csv')
!more ch06\tseries.csv

DatetimeIndex(['2000-01-01', '2000-01-02', '2000-01-03', '2000-01-04',
               '2000-01-05', '2000-01-06', '2000-01-07'],
              dtype='datetime64[ns]', freq='D')






2000-01-01    0
2000-01-02    1
2000-01-03    2
2000-01-04    3
2000-01-05    4
2000-01-06    5
2000-01-07    6
Freq: D, dtype: int32



2000-01-01,0
2000-01-02,1
2000-01-03,2
2000-01-04,3
2000-01-05,4
2000-01-06,5
2000-01-07,6

1	Series.from_csv('ch06/tseries.csv', parse_dates=True)

2000-01-01    0
2000-01-02    1
2000-01-03    2
2000-01-04    3
2000-01-05    4
2000-01-06    5
2000-01-07    6
dtype: int64

手动处理分隔符格式

1	!more ch06\ex7.csv

"a","b","c"
"1","2","3"
"1","2","3","4"

import csv
f = open('ch06/ex7.csv')
reader = csv.reader(f)

1 2	for line in reader: print(line)

['a', 'b', 'c']
['1', '2', '3']
['1', '2', '3', '4']

lines = list(csv.reader(open('ch06/ex7.csv')))
header, values = lines[0], lines[1:]
for item in zip(*values):
    print(item)
    
for h, v in zip(header, zip(*values)):
    print(h, v)
data_dict = {h: v for h, v in zip(header, zip(*values))}
data_dict

('1', '1')
('2', '2')
('3', '3')
a ('1', '1')
b ('2', '2')
c ('3', '3')





{'a': ('1', '1'), 'b': ('2', '2'), 'c': ('3', '3')}

class my_dialect(csv.Dialect):
    lineterminator = '\n'
    delimiter = ';'
    quotechar = '"'
    quoting = csv.QUOTE_MINIMAL

with open('mydata.csv', 'w') as f:
    writer = csv.writer(f, dialect=my_dialect)
    writer.writerow(('one', 'two', 'three'))
    writer.writerow(('1', '2', '3'))
    writer.writerow(('4', '5', '6'))
    writer.writerow(('7', '8', '9'))

1	!more mydata.csv

one;two;three
1;2;3
4;5;6
7;8;9

JSON数据

obj = """
{"name": "Wes",
 "places_lived": ["United States", "Spain", "Germany"],
 "pet": null,
 "siblings": [{"name": "Scott", "age": 25, "pet": "Zuko"},
              {"name": "Katie", "age": 33, "pet": "Cisco"}]
}
"""

1
2
3

import json
result = json.loads(obj)
result

{'name': 'Wes',
 'pet': None,
 'places_lived': ['United States', 'Spain', 'Germany'],
 'siblings': [{'age': 25, 'name': 'Scott', 'pet': 'Zuko'},
  {'age': 33, 'name': 'Katie', 'pet': 'Cisco'}]}

1	asjson = json.dumps(result) # convert to json

1 2	siblings = DataFrame(result['siblings'], columns=['name', 'age']) siblings

	name	age
0	Scott	25
1	Katie	33

XML和HTML： Web信息收集

NB. The Yahoo! Finance API has changed and this example no longer works

from lxml.html import parse
from urllib.request import urlopen
parsed = parse(urlopen('http://finance.yahoo.com/q/op?s=AAPL+Options'))
doc = parsed.getroot()

1 2	links = doc.findall('.//a') links[15:20]

[<Element a at 0x20352cad598>,
 <Element a at 0x20352cad5e8>,
 <Element a at 0x20352cad638>,
 <Element a at 0x20352cad688>,
 <Element a at 0x20352cad6d8>]

lnk = links[28]
lnk
lnk.get('href')
lnk.text_content()

<Element a at 0x20352cad9a8>






'/quote/NFLX?p=NFLX'






'NFLX'

1 2	urls = [lnk.get('href') for lnk in doc.findall('.//a')] urls[-10:]

['//finance.yahoo.com/broker-comparison?bypass=true',
 'https://help.yahoo.com/kb/index?page=content&y=PROD_MAIL_ML&locale=en_US&id=SLN2310&actp=productlink',
 'http://help.yahoo.com/l/us/yahoo/finance/',
 'https://yahoo.uservoice.com/forums/382977',
 'http://info.yahoo.com/privacy/us/yahoo/',
 'http://info.yahoo.com/relevantads/',
 'http://info.yahoo.com/legal/us/yahoo/utos/utos-173.html',
 'http://twitter.com/YahooFinance',
 'http://facebook.com/yahoofinance',
 'http://yahoofinance.tumblr.com']

1
2
3

tables = doc.findall('.//table')
len(tables)
calls = tables[0]

1	rows = calls.findall('.//tr')

1
2
3

def _unpack(row, kind='td'):
    elts = row.findall('.//%s' % kind)
    return [val.text_content() for val in elts]

1 2	_unpack(rows[0], kind='th') _unpack(rows[1], kind='td')

[]




---------------------------------------------------------------------------

IndexError                                Traceback (most recent call last)

<ipython-input-87-7d371ed47023> in <module>()
      1 _unpack(rows[0], kind='th')
----> 2 _unpack(rows[1], kind='td')


IndexError: list index out of range

from pandas.io.parsers import TextParser
def parse_options_data(table):
    rows = table.findall('.//tr')
    header = _unpack(rows[0], kind='th')
    data = [_unpack(r) for r in rows[1:]]
    return TextParser(data, names=header).get_chunk()

1
2
3

call_data = parse_options_data(calls)
put_data = parse_options_data(puts)
call_data[:10]

Parsing XML with lxml.objectify

from lxml import objectify
path = '.\ch06\mta_perf\Performance_MNR.xml'
parsed = objectify.parse(open(path))
root = parsed.getroot()

data = []
skip_fields = ['PARENT_SEQ', 'INDICATOR_SEQ',
               'DESIRED_CHANGE', 'DECIMAL_PLACES']
for elt in root.INDICATOR:
    el_data = {}
    for child in elt.getchildren():
        if child.tag in skip_fields:
            continue
        el_data[child.tag] = child.pyval
    data.append(el_data)

1 2	perf = DataFrame(data) perf.ix[:10, 5:]

	INDICATOR_UNIT	MONTHLY_ACTUAL	MONTHLY_TARGET	PERIOD_MONTH	PERIOD_YEAR	YTD_ACTUAL	YTD_TARGET
0	%	96.9	95	1	2008	96.9	95
1	%	95	95	2	2008	96	95
2	%	96.9	95	3	2008	96.3	95
3	%	98.3	95	4	2008	96.8	95
4	%	95.8	95	5	2008	96.6	95
5	%	94.4	95	6	2008	96.2	95
6	%	96	95	7	2008	96.2	95
7	%	96.4	95	8	2008	96.2	95
8	%	93.7	95	9	2008	95.9	95
9	%	96.4	95	10	2008	96	95
10	%	96.9	95	11	2008	96.1	95

二进制数据格式

1
2
3

frame = pd.read_csv('ch06/ex1.csv')
frame
frame.to_pickle('ch06/frame_pickle')

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

1	pd.read_pickle('ch06/frame_pickle')

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

使用HDF5格式

store = pd.HDFStore('mydata.h5')
store['obj1'] = frame
store['obj1_col'] = frame['a']
store

<class 'pandas.io.pytables.HDFStore'>
File path: mydata.h5
/obj1                frame        (shape->[3,5])
/obj1_col            series       (shape->[3])

1	store['obj1']

	a	b	c	d	message
0	1	2	3	4	hello
1	5	6	7	8	world
2	9	10	11	12	foo

1 2	store.close() os.remove('mydata.h5')

使用数据库

import sqlite3
query = """
CREATE TABLE test
(a VARCHAR(20), b VARCHAR(20),
 c REAL,        d INTEGER
);"""
con = sqlite3.connect(':memory:')
con.execute(query)
con.commit()

<sqlite3.Cursor at 0x2035487c880>

data = [('Atlanta', 'Georgia', 1.25, 6),
        ('Tallahassee', 'Florida', 2.6, 3),
        ('Sacramento', 'California', 1.7, 5)]
stmt = "INSERT INTO test VALUES(?, ?, ?, ?)"
con.executemany(stmt, data)
con.commit()

<sqlite3.Cursor at 0x2035487c810>

1
2
3

cursor = con.execute('select * from test')
rows = cursor.fetchall()
rows

[('Atlanta', 'Georgia', 1.25, 6),
 ('Tallahassee', 'Florida', 2.6, 3),
 ('Sacramento', 'California', 1.7, 5)]

1 2	import pandas.io.sql as sql sql.read_sql('select * from test', con)

	a	b	c	d
0	Atlanta	Georgia	1.25	6
1	Tallahassee	Florida	2.60	3
2	Sacramento	California	1.70	5

Python data analysis Learning note ch05

Posted on 2017-02-28 | | Visitors

pandas入门

按照以下约定引用相关package

1 2	from pandas import Series, DataFrame import pandas as pd

from __future__ import division
from numpy.random import randn
import numpy as np
import os
import matplotlib.pyplot as plt
np.random.seed(12345)
plt.rc('figure', figsize=(10, 6))
from pandas import Series, DataFrame
import pandas as pd
np.set_printoptions(precision=4)
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

pandas数据结构介绍

Series

Series是一种类似于一维数组的对象，由一组数据以及一组与之相关的数据标签（类似于字典的键）组成，所以可以看成是一个有序的字典

1 2	obj = Series([4, 7, -5, 3]) obj

0    4
1    7
2   -5
3    3
dtype: int64

1 2	obj.values obj.index

array([ 4,  7, -5,  3], dtype=int64)






RangeIndex(start=0, stop=4, step=1)

1 2	obj2 = Series([4, 7, -5, 3], index=['d', 'b', 'a', 'c']) obj2

d    4
b    7
a   -5
c    3
dtype: int64

1	obj2.index

Index(['d', 'b', 'a', 'c'], dtype='object')

obj2['a']

-5

1 2	obj2['d'] = 6 obj2[['c', 'a', 'd']]

c    3
a   -5
d    6
dtype: int64

各种Numpy运算都是作用在数据上，同时索引与数据的链接会一直保持

1	obj2[obj2 > 0]

d    6
b    7
c    3
dtype: int64

obj2 * 2

d    12
b    14
a   -10
c     6
dtype: int64

1	np.exp(obj2)

d     403.428793
b    1096.633158
a       0.006738
c      20.085537
dtype: float64

1	'b' in obj2

True

1	'e' in obj2

False

因此可以直接根据Numpy的Dict来创建Series

1
2
3

sdata = {'Ohio': 35000, 'Texas': 71000, 'Oregon': 16000, 'Utah': 5000}
obj3 = Series(sdata)
obj3

Ohio      35000
Oregon    16000
Texas     71000
Utah       5000
dtype: int64

如果传入了index参数的话，那么就会与传入的Dict做键匹配，没有匹配上的就设为NaN

1
2
3

states = ['California', 'Ohio', 'Oregon', 'Texas']
obj4 = Series(sdata, index=states)
obj4

California        NaN
Ohio          35000.0
Oregon        16000.0
Texas         71000.0
dtype: float64

1	pd.isnull(obj4)

California     True
Ohio          False
Oregon        False
Texas         False
dtype: bool

1	pd.notnull(obj4)

California    False
Ohio           True
Oregon         True
Texas          True
dtype: bool

1	obj4.isnull()

California     True
Ohio          False
Oregon        False
Texas         False
dtype: bool

Series有一个非常重要的数据对齐的功能

obj3

Ohio      35000
Oregon    16000
Texas     71000
Utah       5000
dtype: int64

obj4

California        NaN
Ohio          35000.0
Oregon        16000.0
Texas         71000.0
dtype: float64

1	obj3 + obj4

California         NaN
Ohio           70000.0
Oregon         32000.0
Texas         142000.0
Utah               NaN
dtype: float64

Series本身以及其索引都有一个叫做name的属性，这个属性十分重要，以后很多高级功能都会用到

1
2
3

obj4.name = 'population'
obj4.index.name = 'state'
obj4

state
California        NaN
Ohio          35000.0
Oregon        16000.0
Texas         71000.0
Name: population, dtype: float64

可以通过直接赋值的方式修改index属性

1 2	obj.index = ['Bob', 'Steve', 'Jeff', 'Ryan'] obj

Bob      4
Steve    7
Jeff    -5
Ryan     3
dtype: int64

DataFrame

DataFrame是一个表格型的数据结构，可以看成由Series组成的字典，只不过这些Series共用一套索引

data = {'state': ['Ohio', 'Ohio', 'Ohio', 'Nevada', 'Nevada'],
        'year': [2000, 2001, 2002, 2001, 2002],
        'pop': [1.5, 1.7, 3.6, 2.4, 2.9]}
frame = DataFrame(data)

数据会被排序

frame

	pop	state	year
0	1.5	Ohio	2000
1	1.7	Ohio	2001
2	3.6	Ohio	2002
3	2.4	Nevada	2001
4	2.9	Nevada	2002

1	DataFrame(data, columns=['year', 'state', 'pop'])

	year	state	pop
0	2000	Ohio	1.5
1	2001	Ohio	1.7
2	2002	Ohio	3.6
3	2001	Nevada	2.4
4	2002	Nevada	2.9

1
2
3

frame2 = DataFrame(data, columns=['year', 'state', 'pop', 'debt'],
                   index=['one', 'two', 'three', 'four', 'five'])
frame2

	year	state	pop	debt
one	2000	Ohio	1.5	NaN
two	2001	Ohio	1.7	NaN
three	2002	Ohio	3.6	NaN
four	2001	Nevada	2.4	NaN
five	2002	Nevada	2.9	NaN

1	frame2.columns

Index(['year', 'state', 'pop', 'debt'], dtype='object')

DataFrame每一个Key对应的Value都是一个Series

1	frame2['state']

one        Ohio
two        Ohio
three      Ohio
four     Nevada
five     Nevada
Name: state, dtype: object

1	frame2.year

one      2000
two      2001
three    2002
four     2001
five     2002
Name: year, dtype: int64

注意到name属性也已经被设置好了

ix相当于一个行索引？

1	frame2.ix['three']

year     2002
state    Ohio
pop       3.6
debt      NaN
Name: three, dtype: object

可以利用Numpy的广播功能

1 2	frame2['debt'] = 16.5 frame2

	year	state	pop	debt
one	2000	Ohio	1.5	16.5
two	2001	Ohio	1.7	16.5
three	2002	Ohio	3.6	16.5
four	2001	Nevada	2.4	16.5
five	2002	Nevada	2.9	16.5

也可以赋值一个列表，但是长度必须匹配

1 2	frame2['debt'] = np.arange(5.) frame2

	year	state	pop	debt
one	2000	Ohio	1.5	0.0
two	2001	Ohio	1.7	1.0
three	2002	Ohio	3.6	2.0
four	2001	Nevada	2.4	3.0
five	2002	Nevada	2.9	4.0

如果是赋值一个Series，则会匹配上索引，没有匹配上的就是置为NaN

1
2
3

val = Series([-1.2, -1.5, -1.7], index=['two', 'four', 'five'])
frame2['debt'] = val
frame2

	year	state	pop	debt
one	2000	Ohio	1.5	NaN
two	2001	Ohio	1.7	-1.2
three	2002	Ohio	3.6	NaN
four	2001	Nevada	2.4	-1.5
five	2002	Nevada	2.9	-1.7

也可以进行逻辑操作

1 2	frame2['eastern'] = frame2.state == 'Ohio' frame2

	year	state	pop	debt	eastern
one	2000	Ohio	1.5	NaN	True
two	2001	Ohio	1.7	-1.2	True
three	2002	Ohio	3.6	NaN	True
four	2001	Nevada	2.4	-1.5	False
five	2002	Nevada	2.9	-1.7	False

del用于删除一列

1 2	del frame2['eastern'] frame2.columns

Index(['year', 'state', 'pop', 'debt'], dtype='object')

只要是通过索引方式进行的操作，都是直接在原数据上进行的操作，不是一个副本

嵌套的字典也可直接生成DataFrame，只不过内层的键被当作index，外层的键被当作colums

1 2	pop = {'Nevada': {2001: 2.4, 2002: 2.9}, 'Ohio': {2000: 1.5, 2001: 1.7, 2002: 3.6}}

1 2	frame3 = DataFrame(pop) frame3

	Nevada	Ohio
2000	NaN	1.5
2001	2.4	1.7
2002	2.9	3.6

同样可以进行转置，这样的话index和column就会互换

frame3.T

	2000	2001	2002
Nevada	NaN	2.4	2.9
Ohio	1.5	1.7	3.6

显式地指定索引，不匹配的会置为NaN

1	DataFrame(pop, index=[2001, 2002, 2003])

	Nevada	Ohio
2001	2.4	1.7
2002	2.9	3.6
2003	NaN	NaN

frame3

	Nevada	Ohio
2000	NaN	1.5
2001	2.4	1.7
2002	2.9	3.6

还可以这样构建

1
2
3

pdata = {'Ohio': frame3['Ohio'][:-1],
         'Nevada': frame3['Nevada'][:2]}
DataFrame(pdata)

	Nevada	Ohio
2000	NaN	1.5
2001	2.4	1.7

frame3

	Nevada	Ohio
2000	NaN	1.5
2001	2.4	1.7
2002	2.9	3.6

name属性也会在表格中显示出来

1 2	frame3.index.name = 'year'; frame3.columns.name = 'state' frame3

state	Nevada	Ohio
year
2000	NaN	1.5
2001	2.4	1.7
2002	2.9	3.6

values方法只返回数据，不返回index以及key

1	frame3.values

array([[ nan,  1.5],
       [ 2.4,  1.7],
       [ 2.9,  3.6]])

frame2

	year	state	pop	debt
one	2000	Ohio	1.5	NaN
two	2001	Ohio	1.7	-1.2
three	2002	Ohio	3.6	NaN
four	2001	Nevada	2.4	-1.5
five	2002	Nevada	2.9	-1.7

1	frame2.values

array([[2000, 'Ohio', 1.5, nan],
       [2001, 'Ohio', 1.7, -1.2],
       [2002, 'Ohio', 3.6, nan],
       [2001, 'Nevada', 2.4, -1.5],
       [2002, 'Nevada', 2.9, -1.7]], dtype=object)

索引对象

Index是一个可以单独提取出来的对象

1
2
3

obj = Series(range(3), index=['a', 'b', 'c'])
index = obj.index
index

Index(['a', 'b', 'c'], dtype='object')

index[1:]

Index(['b', 'c'], dtype='object')

不可修改~！

1	index[1] = 'd'

---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

<ipython-input-52-676fdeb26a68> in <module>()
----> 1 index[1] = 'd'


C:\Users\Ewan\Anaconda3\lib\site-packages\pandas\indexes\base.py in __setitem__(self, key, value)
   1243 
   1244     def __setitem__(self, key, value):
-> 1245         raise TypeError("Index does not support mutable operations")
   1246 
   1247     def __getitem__(self, key):


TypeError: Index does not support mutable operations

直接创建Index对象

1
2
3

index = pd.Index(np.arange(3))
obj2 = Series([1.5, -2.5, 0], index=index)
obj2.index is index

True

frame3

state	Nevada	Ohio
year
2000	NaN	1.5
2001	2.4	1.7
2002	2.9	3.6

1	'Ohio' in frame3.columns

True

1	2003 in frame3.index

False

基本功能

重新索引

1 2	obj = Series([4.5, 7.2, -5.3, 3.6], index=['d', 'b', 'a', 'c']) obj

d    4.5
b    7.2
a   -5.3
c    3.6
dtype: float64

重排索引形成新对象

1 2	obj2 = obj.reindex(['a', 'b', 'c', 'd', 'e']) obj2

a   -5.3
b    7.2
c    3.6
d    4.5
e    NaN
dtype: float64

1	obj.reindex(['a', 'b', 'c', 'd', 'e'], fill_value=0)

a   -5.3
b    7.2
c    3.6
d    4.5
e    0.0
dtype: float64

1 2	obj3 = Series(['blue', 'purple', 'yellow'], index=[0, 2, 4]) obj3.reindex(range(6), method='ffill') # 前向填充

0      blue
1      blue
2    purple
3    purple
4    yellow
5    yellow
dtype: object

1
2
3

frame = DataFrame(np.arange(9).reshape((3, 3)), index=['a', 'c', 'd'],
                  columns=['Ohio', 'Texas', 'California'])
frame

	Ohio	Texas	California
a	0	1	2
c	3	4	5
d	6	7	8

1 2	frame2 = frame.reindex(['a', 'b', 'c', 'd']) frame2

	Ohio	Texas	California
a	0.0	1.0	2.0
b	NaN	NaN	NaN
c	3.0	4.0	5.0
d	6.0	7.0	8.0

1 2	states = ['Texas', 'Utah', 'California'] frame.reindex(columns=states)

	Texas	Utah	California
a	1	NaN	2
c	4	NaN	5
d	7	NaN	8

插值只能按行

1 2	frame.reindex(index=['a', 'b', 'c', 'd'], method='ffill', columns=states)

	Texas	Utah	California
a	1	NaN	2
b	1	NaN	2
c	4	NaN	5
d	7	NaN	8

用ix方法进行重新索引操作会使得代码很简洁

1	frame.ix[['a', 'b', 'c', 'd'], states]

	Texas	Utah	California
a	1.0	NaN	2.0
b	NaN	NaN	NaN
c	4.0	NaN	5.0
d	7.0	NaN	8.0

丢弃指定轴上的项

1	frame.ix[['a', 'b', 'c', 'd'], states]

	Texas	Utah	California
a	1.0	NaN	2.0
b	NaN	NaN	NaN
c	4.0	NaN	5.0
d	7.0	NaN	8.0

1
2
3

obj = Series(np.arange(5.), index=['a', 'b', 'c', 'd', 'e'])
new_obj = obj.drop('c')
new_obj

a    0.0
b    1.0
d    3.0
e    4.0
dtype: float64

1	obj.drop(['d', 'c'])

a    0.0
b    1.0
e    4.0
dtype: float64

data = DataFrame(np.arange(16).reshape((4, 4)),
                 index=['Ohio', 'Colorado', 'Utah', 'New York'],
                 columns=['one', 'two', 'three', 'four'])
data

	one	two	three	four
Ohio	0	1	2	3
Colorado	4	5	6	7
Utah	8	9	10	11
New York	12	13	14	15

1	data.drop(['Colorado', 'Ohio'])

	one	two	three	four
Utah	8	9	10	11
New York	12	13	14	15

1	data.drop('two', axis=1)

	one	three	four
Ohio	0	2	3
Colorado	4	6	7
Utah	8	10	11
New York	12	14	15

1	data.drop(['two', 'four'], axis=1)

	one	three
Ohio	0	2
Colorado	4	6
Utah	8	10
New York	12	14

索引，选取和过滤

多种索引方式

1
2
3

obj = Series(np.arange(4.), index=['a', 'b', 'c', 'd'])
obj
obj['b']

a    0.0
b    1.0
c    2.0
d    3.0
dtype: float64






1.0

obj[1]

1.0

obj[2:4]

c    2.0
d    3.0
dtype: float64

1	obj[['b', 'a', 'd']]

b    1.0
a    0.0
d    3.0
dtype: float64

1	obj[[1, 3]]

b    1.0
d    3.0
dtype: float64

1	obj[obj < 2] # 直接对data进行操作

a    0.0
b    1.0
dtype: float64

这种切片方式…
末端包含

1	obj['b':'c']

b    1.0
c    2.0
dtype: float64

1 2	obj['b':'c'] = 5 obj

a    0.0
b    5.0
c    5.0
d    3.0
dtype: float64

data = DataFrame(np.arange(16).reshape((4, 4)),
                 index=['Ohio', 'Colorado', 'Utah', 'New York'],
                 columns=['one', 'two', 'three', 'four'])
data

	one	two	three	four
Ohio	0	1	2	3
Colorado	4	5	6	7
Utah	8	9	10	11
New York	12	13	14	15

1	data['two']

Ohio         1
Colorado     5
Utah         9
New York    13
Name: two, dtype: int32

1	data[['three', 'one']]

	three	one
Ohio	2	0
Colorado	6	4
Utah	10	8
New York	14	12

1	data[:2] # axis=0

	one	two	three	four
Ohio	0	1	2	3
Colorado	4	5	6	7

1	data[data['three'] > 5]

	one	two	three	four
Colorado	4	5	6	7
Utah	8	9	10	11
New York	12	13	14	15

data < 5

	one	two	three	four
Ohio	True	True	True	True
Colorado	True	False	False	False
Utah	False	False	False	False
New York	False	False	False	False

1	data[data < 5] = 0

data

	one	two	three	four
Ohio	0	0	0	0
Colorado	0	5	6	7
Utah	8	9	10	11
New York	12	13	14	15

索引的另外一种选择

1	data.ix['Colorado', ['two', 'three']]

two      5
three    6
Name: Colorado, dtype: int32

1	data.ix[['Colorado', 'Utah'], [3, 0, 1]]

	four	one	two
Colorado	7	0	5
Utah	11	8	9

1	data.ix[2]

one       8
two       9
three    10
four     11
Name: Utah, dtype: int32

1	data.ix[:'Utah', 'two']

Ohio        0
Colorado    5
Utah        9
Name: two, dtype: int32

1	data.ix[data.three > 5, :3]

	one	two	three
Colorado	0	5	6
Utah	8	9	10
New York	12	13	14

总结一下就是说，DataFrame是一个二维的数组，只不过每一维的索引方式除了序号之外，还可以用name属性来进行索引，且一切行为与序号无异

算术运算和数据对齐

1 2	s1 = Series([7.3, -2.5, 3.4, 1.5], index=['a', 'c', 'd', 'e']) s2 = Series([-2.1, 3.6, -1.5, 4, 3.1], index=['a', 'c', 'e', 'f', 'g'])

s1

a    7.3
c   -2.5
d    3.4
e    1.5
dtype: float64

s2

a   -2.1
c    3.6
e   -1.5
f    4.0
g    3.1
dtype: float64

数据对齐操作就是一种特殊的并集操作

s1 + s2

a    5.2
c    1.1
d    NaN
e    0.0
f    NaN
g    NaN
dtype: float64

df1 = DataFrame(np.arange(9.).reshape((3, 3)), columns=list('bcd'),
                index=['Ohio', 'Texas', 'Colorado'])
df2 = DataFrame(np.arange(12.).reshape((4, 3)), columns=list('bde'),
                index=['Utah', 'Ohio', 'Texas', 'Oregon'])
df1

	b	c	d
Ohio	0.0	1.0	2.0
Texas	3.0	4.0	5.0
Colorado	6.0	7.0	8.0

df2

	b	d	e
Utah	0.0	1.0	2.0
Ohio	3.0	4.0	5.0
Texas	6.0	7.0	8.0
Oregon	9.0	10.0	11.0

并且数据对齐操作是在所有维度上同时进行的

df1 + df2

	b	c	d	e
Colorado	NaN	NaN	NaN	NaN
Ohio	3.0	NaN	6.0	NaN
Oregon	NaN	NaN	NaN	NaN
Texas	9.0	NaN	12.0	NaN
Utah	NaN	NaN	NaN	NaN

在算术方法中填充词

下面这种定义column的方式值得注意

1
2
3

df1 = DataFrame(np.arange(12.).reshape((3, 4)), columns=list('abcd'))
df2 = DataFrame(np.arange(20.).reshape((4, 5)), columns=list('abcde'))
df1

	a	b	c	d
0	0.0	1.0	2.0	3.0
1	4.0	5.0	6.0	7.0
2	8.0	9.0	10.0	11.0

df2

	a	b	c	d	e
0	0.0	1.0	2.0	3.0	4.0
1	5.0	6.0	7.0	8.0	9.0
2	10.0	11.0	12.0	13.0	14.0
3	15.0	16.0	17.0	18.0	19.0

df1 + df2

	a	b	c	d	e
0	0.0	2.0	4.0	6.0	NaN
1	9.0	11.0	13.0	15.0	NaN
2	18.0	20.0	22.0	24.0	NaN
3	NaN	NaN	NaN	NaN	NaN

要想填充值必须使用add方法

1	df1.add(df2, fill_value=0)

	a	b	c	d	e
0	0.0	2.0	4.0	6.0	4.0
1	9.0	11.0	13.0	15.0	9.0
2	18.0	20.0	22.0	24.0	14.0
3	15.0	16.0	17.0	18.0	19.0

reindex方法与add方法还是存在差异的

1	df1.reindex(columns=df2.columns, fill_value=0)

	a	b	c	d
0	0.0	1.0	2.0	3.0
1	4.0	5.0	6.0	7.0
2	8.0	9.0	10.0	11.0

DataFrame和Series之间的运算

1 2	arr = np.arange(12.).reshape((3, 4)) arr

array([[  0.,   1.,   2.,   3.],
       [  4.,   5.,   6.,   7.],
       [  8.,   9.,  10.,  11.]])

arr[0]

array([ 0.,  1.,  2.,  3.])

广播操作

1	arr - arr[0]

array([[ 0.,  0.,  0.,  0.],
       [ 4.,  4.,  4.,  4.],
       [ 8.,  8.,  8.,  8.]])

frame = DataFrame(np.arange(12.).reshape((4, 3)), columns=list('bde'),
                  index=['Utah', 'Ohio', 'Texas', 'Oregon'])
series = frame.ix[0]
frame

	b	d	e
Utah	0.0	1.0	2.0
Ohio	3.0	4.0	5.0
Texas	6.0	7.0	8.0
Oregon	9.0	10.0	11.0

Series的name等于DataFrame的切片属性

series

b    0.0
d    1.0
e    2.0
Name: Utah, dtype: float64

默认情况下，DataFrame和Series之间的算术运算会将Series的index匹配到DataFrame的column, 然后沿着行向下广播

1	frame - series

	b	d	e
Utah	0.0	0.0	0.0
Ohio	3.0	3.0	3.0
Texas	6.0	6.0	6.0
Oregon	9.0	9.0	9.0

如果Series的index与DataFrame的column不匹配，则进行数据对齐

1 2	series2 = Series(range(3), index=['b', 'e', 'f']) frame + series2

	b	d	e	f
Utah	0.0	NaN	3.0	NaN
Ohio	3.0	NaN	6.0	NaN
Texas	6.0	NaN	9.0	NaN
Oregon	9.0	NaN	12.0	NaN

1 2	series3 = frame['d'] frame

	b	d	e
Utah	0.0	1.0	2.0
Ohio	3.0	4.0	5.0
Texas	6.0	7.0	8.0
Oregon	9.0	10.0	11.0

series3

Utah       1.0
Ohio       4.0
Texas      7.0
Oregon    10.0
Name: d, dtype: float64

匹配行并且在列上进行广播，就必须要指定axis

1 2	frame.sub(series3, axis=0) frame.sub(series3, axis=1)

	b	e
Utah	-1.0	1.0
Ohio	-1.0	1.0
Texas	-1.0	1.0
Oregon	-1.0	1.0

	Ohio	Oregon	Texas	Utah	b	d	e
Utah	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Ohio	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Texas	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Oregon	NaN	NaN	NaN	NaN	NaN	NaN	NaN

函数应用和映射

1 2	frame = DataFrame(np.random.randn(4, 3), columns=list('bde'), index=['Utah', 'Ohio', 'Texas', 'Oregon'])

frame

	b	d	e
Utah	-0.204708	0.478943	-0.519439
Ohio	-0.555730	1.965781	1.393406
Texas	0.092908	0.281746	0.769023
Oregon	1.246435	1.007189	-1.296221

Numpy的元素级方法也可以应用到DataFrame上，直接把DataFrame当作二维的Numpy.array即可

1	np.abs(frame)

	b	d	e
Utah	0.204708	0.478943	0.519439
Ohio	0.555730	1.965781	1.393406
Texas	0.092908	0.281746	0.769023
Oregon	1.246435	1.007189	1.296221

1	f = lambda x: x.max() - x.min()

apply方法将函数映射到由各行或者各列形成的一维数组上

1	frame.apply(f) # axis=0

b    1.802165
d    1.684034
e    2.689627
dtype: float64

1	frame.apply(f, axis=1)

Utah      0.998382
Ohio      2.521511
Texas     0.676115
Oregon    2.542656
dtype: float64

def f(x):
    return Series([x.min(), x.max()], index=['min', 'max'])
frame.apply(f)
frame.apply(f, axis=1)

	b	d	e
min	-0.555730	0.281746	-1.296221
max	1.246435	1.965781	1.393406

	min	max
Utah	-0.519439	0.478943
Ohio	-0.555730	1.965781
Texas	0.092908	0.769023
Oregon	-1.296221	1.246435

元素级的函数映射applymap，之所以叫这个名字是因为Series有一个元素级的映射函数map

1 2	format = lambda x: '%.2f' % x frame.applymap(format)

	b	d	e
Utah	-0.20	0.48	-0.52
Ohio	-0.56	1.97	1.39
Texas	0.09	0.28	0.77
Oregon	1.25	1.01	-1.30

1	frame['e'].map(format)

Utah      -0.52
Ohio       1.39
Texas      0.77
Oregon    -1.30
Name: e, dtype: object

排序和排名

对索引或者column进行（字典）排序

1 2	obj = Series(range(4), index=['d', 'a', 'b', 'c']) obj.sort_index()

a    1
b    2
c    3
d    0
dtype: int32

1
2
3

frame = DataFrame(np.arange(8).reshape((2, 4)), index=['three', 'one'],
                  columns=['d', 'a', 'b', 'c'])
frame.sort_index()

	d	a	b	c
one	4	5	6	7
three	0	1	2	3

1	frame.sort_index(axis=1)

	a	b	c	d
three	1	2	3	0
one	5	6	7	4

降序

1	frame.sort_index(axis=1, ascending=False)

	d	c	b	a
three	0	3	2	1
one	4	7	6	5

按照data进行排序

1 2	obj = Series([4, 7, -3, 2]) obj.sort_values()

2   -3
3    2
0    4
1    7
dtype: int64

在排序时，任何缺失值默认都会被放到Series的末尾

1 2	obj = Series([4, np.nan, 7, np.nan, -3, 2]) obj.sort_values()

4   -3.0
5    2.0
0    4.0
2    7.0
1    NaN
3    NaN
dtype: float64

1 2	frame = DataFrame({'b': [4, 7, -3, 2], 'a': [0, 1, 0, 1]}) frame

	a	b
0	0	4
1	1	7
2	0	-3
3	1	2

对指定index或者column进行排序

1	frame.sort_values(by='b')

	a	b
2	0	-3
3	1	2
0	0	4
1	1	7

或者根据multi-index亦或multi-column进行排序

1	frame.sort_values(by=['a', 'b'])

	a	b
2	0	-3
0	0	4
3	1	2
1	1	7

默认情况下，rank方法通过“为各组分配一个平均排名”的方式破坏平级关系。也就是说，如果有多个相同的值，则这些值的rank就是这些相同值rand的算术平均。

1 2	obj = Series([7, -5, 7, 4, 2, 0, 4, 4]) obj.rank()

0    7.5
1    1.0
2    7.5
3    5.0
4    3.0
5    2.0
6    5.0
7    5.0
dtype: float64

如果想按照一般方式排名

1	obj.rank(method='first')

0    7.0
1    1.0
2    8.0
3    4.0
4    3.0
5    2.0
6    5.0
7    6.0
dtype: float64

使用每个分组的最大排名

1	obj.rank(ascending=False, method='max')

0    2.0
1    8.0
2    2.0
3    5.0
4    6.0
5    7.0
6    5.0
7    5.0
dtype: float64

1
2
3

frame = DataFrame({'b': [4.3, 7, -3, 2], 'a': [0, 1, 0, 1],
                   'c': [-2, 5, 8, -2.5]})
frame

	a	b	c
0	0	4.3	-2.0
1	1	7.0	5.0
2	0	-3.0	8.0
3	1	2.0	-2.5

指定维度

1	frame.rank(axis=1)

	a	b	c
0	2.0	3.0	1.0
1	1.0	3.0	2.0
2	2.0	1.0	3.0
3	2.0	3.0	1.0

带有重复值的轴索引

1 2	obj = Series(range(5), index=['a', 'a', 'b', 'b', 'c']) obj

a    0
a    1
b    2
b    3
c    4
dtype: int32

1	obj.index.is_unique

False

返回一个Series

obj['a']

a    0
a    1
dtype: int32

obj['c']

1 2	df = DataFrame(np.random.randn(4, 3), index=['a', 'a', 'b', 'b']) df

	0	1	2
a	0.274992	0.228913	1.352917
a	0.886429	-2.001637	-0.371843
b	1.669025	-0.438570	-0.539741
b	0.476985	3.248944	-1.021228

1	df.ix['b']

	0	1	2
b	1.669025	-0.438570	-0.539741
b	0.476985	3.248944	-1.021228

汇总和计算描述统计

df = DataFrame([[1.4, np.nan], [7.1, -4.5],
                [np.nan, np.nan], [0.75, -1.3]],
               index=['a', 'b', 'c', 'd'],
               columns=['one', 'two'])
df

	one	two
a	1.40	NaN
b	7.10	-4.5
c	NaN	NaN
d	0.75	-1.3

1	df.sum() # axis=0 skipna=True

one    9.25
two   -5.80
dtype: float64

1	df.sum(axis=1) # skipna=True

a    1.40
b    2.60
c    0.00
d   -0.55
dtype: float64

1	df.mean(axis=1, skipna=False)

a      NaN
b    1.300
c      NaN
d   -0.275
dtype: float64

返回的是索引

1	df.idxmax()

one    b
two    d
dtype: object

1	df.cumsum()

	one	two
a	1.40	NaN
b	8.50	-4.5
c	NaN	NaN
d	9.25	-5.8

describe对于数值型和非数值型数据的行为不一样

1	df.describe()

C:\Users\Ewan\Anaconda3\lib\site-packages\numpy\lib\function_base.py:3834: RuntimeWarning: Invalid value encountered in percentile
  RuntimeWarning)

	one	two
count	3.000000	2.000000
mean	3.083333	-2.900000
std	3.493685	2.262742
min	0.750000	-4.500000
25%	NaN	NaN
50%	NaN	NaN
75%	NaN	NaN
max	7.100000	-1.300000

1
2
3

obj = Series(['a', 'a', 'b', 'c'] * 4)
obj
obj.describe()

0     a
1     a
2     b
3     c
4     a
5     a
6     b
7     c
8     a
9     a
10    b
11    c
12    a
13    a
14    b
15    c
dtype: object






count     16
unique     3
top        a
freq       8
dtype: object

唯一值，值计数以及成员资格

1	obj = Series(['c', 'a', 'd', 'a', 'a', 'b', 'b', 'c', 'c'])

1 2	uniques = obj.unique() uniques

array(['c', 'a', 'd', 'b'], dtype=object)

1	obj.value_counts()

c    3
a    3
b    2
d    1
dtype: int64

1	pd.value_counts(obj.values, sort=False)

a    3
d    1
b    2
c    3
dtype: int64

1 2	mask = obj.isin(['b', 'c']) mask

0     True
1    False
2    False
3    False
4    False
5     True
6     True
7     True
8     True
dtype: bool

obj[mask]

0    c
5    b
6    b
7    c
8    c
dtype: object

data = DataFrame({'Qu1': [1, 3, 4, 3, 4],
                  'Qu2': [2, 3, 1, 2, 3],
                  'Qu3': [1, 5, 2, 4, 4]})
data

	Qu1	Qu2	Qu3
0	1	2	1
1	3	3	5
2	4	1	2
3	3	2	4
4	4	3	4

1 2	result = data.apply(pd.value_counts).fillna(0) result

	Qu1	Qu2	Qu3
1	1.0	1.0	1.0
2	0.0	2.0	1.0
3	2.0	2.0	0.0
4	2.0	0.0	2.0
5	0.0	0.0	1.0

处理缺失数据

1 2	string_data = Series(['aardvark', 'artichoke', np.nan, 'avocado']) string_data

0     aardvark
1    artichoke
2          NaN
3      avocado
dtype: object

1	string_data.isnull()

0    False
1    False
2     True
3    False
dtype: bool

1 2	string_data[0] = None string_data.isnull()

0     True
1    False
2     True
3    False
dtype: bool

滤除缺失数据

1
2
3

from numpy import nan as NA
data = Series([1, NA, 3.5, NA, 7])
data.dropna()

0    1.0
2    3.5
4    7.0
dtype: float64

1	data[data.notnull()]

0    1.0
2    3.5
4    7.0
dtype: float64

data = DataFrame([[1., 6.5, 3.], [1., NA, NA],
                  [NA, NA, NA], [NA, 6.5, 3.]])
cleaned = data.dropna()
data

	0	1	2
0	1.0	6.5	3.0
1	1.0	NaN	NaN
2	NaN	NaN	NaN
3	NaN	6.5	3.0

cleaned

	0	1	2
0	1.0	6.5	3.0

1	data.dropna(how='all')

	0	1	2
0	1.0	6.5	3.0
1	1.0	NaN	NaN
3	NaN	6.5	3.0

1 2	data[4] = NA data

	0	1	2	4
0	1.0	6.5	3.0	NaN
1	1.0	NaN	NaN	NaN
2	NaN	NaN	NaN	NaN
3	NaN	6.5	3.0	NaN

1	data.dropna(axis=1, how='all')

	0	1	2
0	1.0	6.5	3.0
1	1.0	NaN	NaN
2	NaN	NaN	NaN
3	NaN	6.5	3.0

1
2
3

df = DataFrame(np.random.randn(7, 3))
df.ix[:4, 1] = NA; df.ix[:2, 2] = NA
df

	0	1	2
0	-0.204708	NaN	NaN
1	-0.555730	NaN	NaN
2	0.092908	NaN	NaN
3	1.246435	NaN	-1.296221
4	0.274992	NaN	1.352917
5	0.886429	-2.001637	-0.371843
6	1.669025	-0.438570	-0.539741

1	df.dropna(thresh=2)

	0	1	2
3	1.246435	NaN	-1.296221
4	0.274992	NaN	1.352917
5	0.886429	-2.001637	-0.371843
6	1.669025	-0.438570	-0.539741

填充缺失数据

1	df.fillna(0)

	0	1	2
0	-0.204708	0.000000	0.000000
1	-0.555730	0.000000	0.000000
2	0.092908	0.000000	0.000000
3	1.246435	0.000000	-1.296221
4	0.274992	0.000000	1.352917
5	0.886429	-2.001637	-0.371843
6	1.669025	-0.438570	-0.539741

1	df.fillna({1: 0.5, 3: -1})

	0	1	2
0	-0.204708	0.500000	NaN
1	-0.555730	0.500000	NaN
2	0.092908	0.500000	NaN
3	1.246435	0.500000	-1.296221
4	0.274992	0.500000	1.352917
5	0.886429	-2.001637	-0.371843
6	1.669025	-0.438570	-0.539741

fillna默认返回新对象，但也可以对现有对象就地修改

1
2
3

# always returns a reference to the filled object
_ = df.fillna(0, inplace=True)
df

	0	1	2
0	-0.204708	0.000000	0.000000
1	-0.555730	0.000000	0.000000
2	0.092908	0.000000	0.000000
3	1.246435	0.000000	-1.296221
4	0.274992	0.000000	1.352917
5	0.886429	-2.001637	-0.371843
6	1.669025	-0.438570	-0.539741

1
2
3

df = DataFrame(np.random.randn(6, 3))
df.ix[2:, 1] = NA; df.ix[4:, 2] = NA
df

	0	1	2
0	0.476985	3.248944	-1.021228
1	-0.577087	0.124121	0.302614
2	0.523772	NaN	1.343810
3	-0.713544	NaN	-2.370232
4	-1.860761	NaN	NaN
5	-1.265934	NaN	NaN

1	df.fillna(method='ffill')

	0	1	2
0	0.476985	3.248944	-1.021228
1	-0.577087	0.124121	0.302614
2	0.523772	0.124121	1.343810
3	-0.713544	0.124121	-2.370232
4	-1.860761	0.124121	-2.370232
5	-1.265934	0.124121	-2.370232

1	df.fillna(method='ffill', limit=2)

	0	1	2
0	0.476985	3.248944	-1.021228
1	-0.577087	0.124121	0.302614
2	0.523772	0.124121	1.343810
3	-0.713544	0.124121	-2.370232
4	-1.860761	NaN	-2.370232
5	-1.265934	NaN	-2.370232

1 2	data = Series([1., NA, 3.5, NA, 7]) data.fillna(data.mean())

0    1.000000
1    3.833333
2    3.500000
3    3.833333
4    7.000000
dtype: float64

层次索引

data = Series(np.random.randn(10),
              index=[['a', 'a', 'a', 'b', 'b', 'b', 'c', 'c', 'd', 'd'],
                     [1, 2, 3, 1, 2, 3, 1, 2, 2, 3]])
data

a  1    0.332883
   2   -2.359419
   3   -0.199543
b  1   -1.541996
   2   -0.970736
   3   -1.307030
c  1    0.286350
   2    0.377984
d  2   -0.753887
   3    0.331286
dtype: float64

1	data.index

MultiIndex(levels=[['a', 'b', 'c', 'd'], [1, 2, 3]],
           labels=[[0, 0, 0, 1, 1, 1, 2, 2, 3, 3], [0, 1, 2, 0, 1, 2, 0, 1, 1, 2]])

data['b']

1   -1.541996
2   -0.970736
3   -1.307030
dtype: float64

1	data['b':'c']

b  1   -1.541996
   2   -0.970736
   3   -1.307030
c  1    0.286350
   2    0.377984
dtype: float64

1	data.ix[['b', 'd']]

b  1   -1.541996
   2   -0.970736
   3   -1.307030
d  2   -0.753887
   3    0.331286
dtype: float64

1	data[:, 2]

a   -2.359419
b   -0.970736
c    0.377984
d   -0.753887
dtype: float64

1	data.unstack()

	1	2	3
a	0.332883	-2.359419	-0.199543
b	-1.541996	-0.970736	-1.307030
c	0.286350	0.377984	NaN
d	NaN	-0.753887	0.331286

1	data.unstack().stack()

a  1    0.332883
   2   -2.359419
   3   -0.199543
b  1   -1.541996
   2   -0.970736
   3   -1.307030
c  1    0.286350
   2    0.377984
d  2   -0.753887
   3    0.331286
dtype: float64

frame = DataFrame(np.arange(12).reshape((4, 3)),
                  index=[['a', 'a', 'b', 'b'], [1, 2, 1, 2]],
                  columns=[['Ohio', 'Ohio', 'Colorado'],
                           ['Green', 'Red', 'Green']])
frame

		Ohio		Colorado
		Green	Red	Green
a	1	0	1	2
a	2	3	4	5
b	1	6	7	8
b	2	9	10	11

1
2
3

frame.index.names = ['key1', 'key2']
frame.columns.names = ['state', 'color']
frame

	state	Ohio		Colorado
	color	Green	Red	Green
key1	key2
a	1	0	1	2
a	2	3	4	5
b	1	6	7	8
b	2	9	10	11

1	frame['Ohio']

	color	Green	Red
key1	key2
a	1	0	1
a	2	3	4
b	1	6	7
b	2	9	10

创建MultiIndex对象复用

1 2	MultiIndex.from_arrays([['Ohio', 'Ohio', 'Colorado'], ['Green', 'Red', 'Green']], names=['state', 'color'])

重排分级顺序

1	frame.swaplevel('key1', 'key2')

	state	Ohio		Colorado
	color	Green	Red	Green
key2	key1
1	a	0	1	2
2	a	3	4	5
1	b	6	7	8
2	b	9	10	11

1	frame.sortlevel(1)

	state	Ohio		Colorado
	color	Green	Red	Green
key1	key2
a	1	0	1	2
b	1	6	7	8
a	2	3	4	5
b	2	9	10	11

1	frame.swaplevel(0, 1).sortlevel(0)

	state	Ohio		Colorado
	color	Green	Red	Green
key2	key1
1	a	0	1	2
1	b	6	7	8
2	a	3	4	5
2	b	9	10	11

根据级别汇总统计

1	frame.sum(level='key2')

state	Ohio		Colorado
color	Green	Red	Green
key2
1	6	8	10
2	12	14	16

1	frame.sum(level='color', axis=1)

	color	Green	Red
key1	key2
a	1	2	1
a	2	8	4
b	1	14	7
b	2	20	10

使用DataFrame的列

frame = DataFrame({'a': range(7), 'b': range(7, 0, -1),
                   'c': ['one', 'one', 'one', 'two', 'two', 'two', 'two'],
                   'd': [0, 1, 2, 0, 1, 2, 3]})
frame

	a	b	c	d
0	0	7	one	0
1	1	6	one	1
2	2	5	one	2
3	3	4	two	0
4	4	3	two	1
5	5	2	two	2
6	6	1	two	3

1 2	frame2 = frame.set_index(['c', 'd']) frame2

		a	b
c	d
one	0	0	7
	1	1	6
	2	2	5
two	0	3	4
	1	4	3
	2	5	2
	3	6	1

1	frame.set_index(['c', 'd'], drop=False)

		a	b	c	d
c	d
one	0	0	7	one	0
	1	1	6	one	1
	2	2	5	one	2
two	0	3	4	two	0
	1	4	3	two	1
	2	5	2	two	2
	3	6	1	two	3

1	frame2.reset_index()

	c	d	a	b
0	one	0	0	7
1	one	1	1	6
2	one	2	2	5
3	two	0	3	4
4	two	1	4	3
5	two	2	5	2
6	two	3	6	1

拓展话题

整数索引

1 2	ser = Series(np.arange(3.)) ser.iloc[-1]

2.0

ser

0    0.0
1    1.0
2    2.0
dtype: float64

1 2	ser2 = Series(np.arange(3.), index=['a', 'b', 'c']) ser2[-1]

2.0

1	ser.ix[:1]

0    0.0
1    1.0
dtype: float64

1 2	ser3 = Series(range(3), index=[-5, 1, 3]) ser3.iloc[2]

1
2
3

frame = DataFrame(np.arange(6).reshape((3, 2)), index=[2, 0, 1])
frame
frame.iloc[0]

	0	1
2	0	1
0	2	3
1	4	5

0    0
1    1
Name: 2, dtype: int32

面板数据

import pandas_datareader.data as web
pdata = pd.Panel(dict((stk, web.get_data_yahoo(stk))
                       for stk in ['AAPL', 'GOOG', 'MSFT', 'DELL']))

pdata

<class 'pandas.core.panel.Panel'>
Dimensions: 4 (items) x 1820 (major_axis) x 6 (minor_axis)
Items axis: AAPL to MSFT
Major_axis axis: 2010-01-04 00:00:00 to 2017-02-27 00:00:00
Minor_axis axis: Open to Adj Close

1 2	pdata = pdata.swapaxes('items', 'minor') pdata['Adj Close'].iloc[:10]

	AAPL	DELL	GOOG	MSFT
Date
2010-01-04	27.727039	14.06528	313.062468	25.555485
2010-01-05	27.774976	14.38450	311.683844	25.563741
2010-01-06	27.333178	14.10397	303.826685	25.406859
2010-01-07	27.282650	14.23940	296.753749	25.142634
2010-01-08	27.464034	14.36516	300.709808	25.316031
2010-01-11	27.221758	14.37483	300.255255	24.994007
2010-01-12	26.912110	14.56830	294.945572	24.828866
2010-01-13	27.291720	14.57797	293.252243	25.060064
2010-01-14	27.133657	14.22005	294.630868	25.563741
2010-01-15	26.680198	13.92985	289.710772	25.481172

1	pdata.ix[:, '6/1/2012', :]

	Open	High	Low	Close	Volume	Adj Close
AAPL	569.159996	572.650009	560.520012	560.989983	130246900.0	72.681610
DELL	12.150000	12.300000	12.045000	12.070000	19397600.0	11.675920
GOOG	571.790972	572.650996	568.350996	570.981000	6138700.0	285.205295
MSFT	28.760000	28.959999	28.440001	28.450001	56634300.0	24.942239

1	pdata.ix['Adj Close', '5/22/2012':, :].iloc[:10]

	AAPL	DELL	GOOG	MSFT
Date
2012-05-22	72.160786	14.58765	300.100412	26.090721
2012-05-23	73.921494	12.08221	304.426106	25.520864
2012-05-24	73.242607	12.04351	301.528978	25.485795
2012-05-25	72.850038	12.05319	295.470050	25.477028
2012-05-28	NaN	12.05319	NaN	NaN
2012-05-29	74.143041	12.24666	296.873645	25.915380
2012-05-30	75.037005	12.14992	293.821674	25.722505
2012-05-31	74.850442	11.92743	290.140354	25.591000
2012-06-01	72.681610	11.67592	285.205295	24.942239
2012-06-04	73.109156	11.60821	289.006480	25.029908

1 2	stacked = pdata.ix[:, '5/30/2012':, :].to_frame() stacked

		Open	High	Low	Close	Volume	Adj Close
Date	minor
2012-05-30	AAPL	569.199997	579.989990	566.559990	579.169998	132357400.0	75.037005
	DELL	12.590000	12.700000	12.460000	12.560000	19787800.0	12.149920
	GOOG	588.161028	591.901014	583.530999	588.230992	3827600.0	293.821674
	MSFT	29.350000	29.480000	29.120001	29.340000	41585500.0	25.722505
2012-05-31	AAPL	580.740021	581.499985	571.460022	577.730019	122918600.0	74.850442
	DELL	12.530000	12.540000	12.330000	12.330000	19955600.0	11.927430
	GOOG	588.720982	590.001032	579.001013	580.860990	5958800.0	290.140354
	MSFT	29.299999	29.420000	28.940001	29.190001	39134000.0	25.591000
2012-06-01	AAPL	569.159996	572.650009	560.520012	560.989983	130246900.0	72.681610
	DELL	12.150000	12.300000	12.045000	12.070000	19397600.0	11.675920
	GOOG	571.790972	572.650996	568.350996	570.981000	6138700.0	285.205295
	MSFT	28.760000	28.959999	28.440001	28.450001	56634300.0	24.942239
2012-06-04	AAPL	561.500008	567.499985	548.499977	564.289978	139248900.0	73.109156
	DELL	12.110000	12.112500	11.800000	12.000000	17015700.0	11.608210
	GOOG	570.220958	580.491016	570.011006	578.590973	4883500.0	289.006480
	MSFT	28.620001	28.780001	28.320000	28.549999	47926300.0	25.029908
2012-06-05	AAPL	561.269989	566.470001	558.330002	562.830025	97053600.0	72.920005
	DELL	11.950000	12.240000	11.950000	12.160000	15620900.0	11.762980
	GOOG	575.451008	578.131003	566.470986	570.410999	4697200.0	284.920579
	MSFT	28.510000	28.750000	28.389999	28.510000	45715400.0	24.994841
2012-06-06	AAPL	567.770004	573.849983	565.499992	571.460022	100363900.0	74.038104
	DELL	12.210000	12.280000	12.090000	12.215000	20779900.0	11.816190
	GOOG	576.480979	581.970971	573.611004	580.570966	4207200.0	289.995487
	MSFT	28.879999	29.370001	28.809999	29.350000	46860500.0	25.731273
2012-06-07	AAPL	577.290009	577.320023	570.500000	571.720001	94941700.0	74.071787
	DELL	12.320000	12.410000	12.120000	12.130000	20074000.0	11.733960
	GOOG	587.601014	587.891038	577.251006	578.230986	3530100.0	288.826666
	MSFT	29.639999	29.700001	29.170000	29.230000	37792800.0	25.626067
2012-06-08	AAPL	571.599998	580.580017	568.999992	580.319984	86879100.0	75.185997
2012-06-08	DELL	12.130000	12.225000	12.020000	12.120000	18155600.0	11.724290
…	…	…	…	…	…	…	…
2017-02-13	AAPL	133.080002	133.820007	132.750000	133.289993	23035400.0	133.289993
	GOOG	816.000000	820.958984	815.489990	819.239990	1198100.0	819.239990
	MSFT	64.239998	64.860001	64.129997	64.720001	22920100.0	64.330000
2017-02-14	AAPL	133.470001	135.089996	133.250000	135.020004	32815500.0	135.020004
	GOOG	819.000000	823.000000	816.000000	820.450012	1053600.0	820.450012
	MSFT	64.410004	64.720001	64.019997	64.570000	23065900.0	64.570000
2017-02-15	AAPL	135.520004	136.270004	134.619995	135.509995	35501600.0	135.509995
	GOOG	819.359985	823.000000	818.469971	818.979980	1304000.0	818.979980
	MSFT	64.500000	64.570000	64.160004	64.529999	16917000.0	64.529999
2017-02-16	AAPL	135.669998	135.899994	134.839996	135.350006	22118000.0	135.350006
	GOOG	819.929993	824.400024	818.979980	824.159973	1281700.0	824.159973
	MSFT	64.739998	65.239998	64.440002	64.519997	20524700.0	64.519997
2017-02-17	AAPL	135.100006	135.830002	135.100006	135.720001	22084500.0	135.720001
	GOOG	823.020020	828.070007	821.655029	828.070007	1597800.0	828.070007
	MSFT	64.470001	64.690002	64.300003	64.620003	21234600.0	64.620003
2017-02-21	AAPL	136.229996	136.750000	135.979996	136.699997	24265100.0	136.699997
	GOOG	828.659973	833.450012	828.349976	831.659973	1247700.0	831.659973
	MSFT	64.610001	64.949997	64.449997	64.489998	19384900.0	64.489998
2017-02-22	AAPL	136.429993	137.119995	136.110001	137.110001	20745300.0	137.110001
	GOOG	828.659973	833.250000	828.640015	830.760010	982900.0	830.760010
	MSFT	64.330002	64.389999	64.050003	64.360001	19259700.0	64.360001
2017-02-23	AAPL	137.380005	137.479996	136.300003	136.529999	20704100.0	136.529999
	GOOG	830.119995	832.460022	822.880005	831.330017	1470100.0	831.330017
	MSFT	64.419998	64.730003	64.190002	64.620003	20235200.0	64.620003
2017-02-24	AAPL	135.910004	136.660004	135.279999	136.660004	21690900.0	136.660004
	GOOG	827.729980	829.000000	824.200012	828.640015	1386600.0	828.640015
	MSFT	64.529999	64.800003	64.139999	64.620003	21705200.0	64.620003
2017-02-27	AAPL	137.139999	137.440002	136.279999	136.929993	20196400.0	136.929993
	GOOG	824.549988	830.500000	824.000000	829.280029	1099500.0	829.280029
	MSFT	64.540001	64.540001	64.050003	64.230003	15850400.0	64.230003

3952 rows × 6 columns

1	stacked.to_panel()

<class 'pandas.core.panel.Panel'>
Dimensions: 6 (items) x 1207 (major_axis) x 4 (minor_axis)
Items axis: Open to Adj Close
Major_axis axis: 2012-05-30 00:00:00 to 2017-02-27 00:00:00
Minor_axis axis: AAPL to MSFT

Correlation Analysis

Posted on 2017-02-24 | | Visitors

问题描述

最近有一个关于关联分析的小作业，问题描述如下：

孩子的情感状况、品行症状、亲社会行为等等与很多因素相关，这个作业主要着眼点在于孩子每天看视频的时间是否对于以上这些指标有着重要的影响，因此需要对输入特征与输出标签之间的相关关系做一个分析。

属性描述如下：

由以上属性描述表可以看出，属性类型同时包括类别以及数值，因此进行相关分析时针对这两种情况需要采用不同的方法，这里采用的策略如下：

对于离散-离散的情况，采用卡方检验
对于连续-离散的情况，采用ANOVA方差检验

下面对于这两种检验方法进行介绍。

卡方检验

一维情况

假设一条河里有三种鱼，gumpies, sticklebarbs, 以及spotheads。如果这条河的生态环境没有遭到干扰，那么这三种鱼的数量是相等的（下表第三行）。现在从河里进行300次抽样，最后抽样得到的结果如下表第二行所示：

	gumpies	sticklebarbs	spothheads	Totals
观察到的频数	89	120	91	300
期望的频数	100	100	100	300

现在需要解决的问题是，是否这条河的生态环境收到了干扰？我们假设生态环境正常，这是我们的原假设。

（注意，期望的频数是根据实际情况自行定义的）

很容易想到的是，我们可以建立一种度量来衡量现实情况与原假设的偏离程度，例如：

$$\frac{观察值 - 期望值}{期望值}$$

把实际的数据带入，可以得到以下的结果：

gumpies: $\frac{89 - 100}{100} = -0.11$

sticklebarbs: $\frac{120 - 100}{100} = +0.20$

spothheads: $\frac{91 - 100}{100} = -0.09$

发现结果还不错，但是这只能用来衡量单个类别的偏差程度，而不能用来衡量整体的偏差程度，因为这三者的加和为零。既然这样，很容易想到可以对之前的度量进行简单的修改，变成这样：

$$\frac{(观察值 - 期望值)^2}{期望值}$$

再把数据带入看看：

gumpies: $\frac{(89 - 100)^2}{100} = 1.21$

sticklebarbs: $\frac{(120 - 100)^2}{100} = 4.0$

spothheads: $\frac{(91 - 100)}{100} = 0.81$

sum: $1.21 + 4.0 + 0.81 = 6.02$

这样一来问题便得到了解决，而这也正是卡方检验所采取的方式。而这个sum值就是卡方（chi-square），记作$\chi^2$，为了更加形式化地表示，我们把观察值记为$O$，期望值记为$E$，那么有如下等式成立：

$$\chi^2 = \sum \frac{(O - E) ^ 2}{E}$$

下面的问题是，我们达到了卡方值，但是这个值到底好还是不好呢？是更加支持原假设还是拒绝原假设呢？

可以设想这样一种情况，我们假设河里的鱼服从原假设的分布，也就是三种鱼出现的概率相等。然后我们把三百次采样看作一次实验，每一次实验完毕之后记录下采样出的300条鱼中每一种鱼的频数，然后计算卡方值。在进行很多次这样的实验之后，我们可以画一个柱状图，这个图记录下了卡方值的分布情况。然后再把实际的观察值（也就是上面表格的第二行）计算的卡方值（6.02）带入进去，看看大于或者等于这个值在柱状图中所有卡方值中占有多少比例。如果占有的比例很大，说明这个值是由跟原假设很近似的假设得出的，这就证明了原假设是对的；反之。如果这个比例很小，说明如果分布服从原假设，那么所计算出的卡方值基本不可能包含这个观测出的卡方值，表明原假设是不对的，我们就可以拒绝原假设。

其实统计检验的基本思想就是这样。但是存在一个问题，我们不可能进行真实的采样（从河里抓鱼），所以一般采用计算机模拟的方式，具体步骤如下所示：

等概率的产生a（代表gumpies）, b（代表sticklebarbs）, c（代表spothheads）的序列
计算一个大小为300的序列中a, b, c三者的频数，作为观察值，然后将a = b = c = 100作为期望值，计算并记录下算出的卡方值
重复1~2步10000次，画出柱状图如下：

可以看出只有5%左右的值大于6.02，说明我们可以以95%的置信度拒绝原假设。

二维情况

对于二维的情况，卡方检验又被称为卡方关联度检验，也就是检验两个变量之间的相关程度（独立程度），考虑下面这个数据表$O$：

		Alzhemer’s onset -	during 5-year period
		no	yes
recieved-	yes	147	9	156
estrogen	no	810	158	968
		957	167	1124

这是观察值，为了计算卡方值，很明显我们需要计算期望值。

为了方便表示，把上表变成如下形式：

		（A）Alzhemer’s onset -	during 5-year period
		no	yes
(R) recieved-	yes	[cell a]	[cell b]	156
estrogen	no	[cell c]	[cell d]	968
		957	167	1124

拿a做例子：

$$E_a = \frac{156}{1124} \times \frac{957}{1124} \times 1124$$

这个式子如何解释呢？如果这两个变量是独立的，那么$A$变量取$no$值与$R$变量取$yes$值这两个事件之间就是独立的，那么$[cell_a]$事件发生的概率就是两者相乘，也就是上述等式右边前两个数相乘，最后的期望值自然就是概率乘以实验总数。

上面的解释比较不正式，换一种较为正式的表达方式。假设我们要求$cell$的期望值，设$R$为$cell$所在列的边缘事件总数，$C$为$cell$所在行的边缘事件总数，$N$为实验总数目，这样就有：

$$E_{cell} =\frac{R}{N} \times \frac{C}{N} \times N$$

所以就有期望值$E$数据表如下：

		（A）Alzhemer’s onset -	during 5-year period
		no	yes
(R) recieved-	yes	$E_a = \frac{156 \times 957}{1124} = 132.82$	$E_b = \frac{156 \times 167}{1124} = 23.18$	156
estrogen	no	$E_c = \frac{968 \times 957}{1124} = 824.18$	$E_d = \frac{968 \times 167}{1124} = 143.82$	968
		957	167	1124

这样就可以调用公式：

$$\chi^2 = \sum \frac{(O - E) ^ 2}{E}$$

特别地，当行数以及列数都为2时，上述公式需要进行一下修改：

$$\chi^2 = \sum \frac{(|O - E| - 0.5) ^ 2}{E}$$

这样算出来的卡方值为$11.01$

最后涉及到自由度的问题，因为比较简单，所以只写出结论：

$$df = (r - 1)(c - 1)$$

r = number of rows

c = number of columns

ANOVA方差检验

为了方便比较，同样采用上述阿尔兹海默病的例子。研究表明，老年痴呆症患者患病之后会经常经历情绪非常不稳定地阶段，原因是因为患者患病之前的生活中经常有恐惧或者焦虑的体验，正是这些一直存在于脑海中的记忆出发了患病后的不稳定情绪的产生。

现在我们假设有一个实验团队发明了一种药物，可以缓解这种情绪问题，他们对小白鼠进行了实验。实验设计如下：

将小白鼠随机分为四组A， B， C， D
A组作为参照组，不给药；B， C， D三组分别注射一个单位，两个单位，三个单位的药剂
记录实验结果，数值越低表明实验效果越好

实验结果如下：

A	B	C	D	Total
27.0	22.8	21.9	23.5
26.2	23.1	23.4	19.6
28.8	27.7	20.1	23.7
33.5	27.6	27.8	20.8
28.8	24.0	19.3	23.9
$M_a = 28.86$	$M_b = 25.04$	$M_c = 22.50$	$M_d = 22.30$	$M_T = 24.68$

下面需要进行一下相关性分析，判断药物是否对症状的缓解产生作用。

和卡方检验一样，ANOVA检验最后也有一个衡量指标，记为$F$，定义如下：

$$F = \frac{MS_{bg}}{MS_{wg}} = \frac{组间相似度}{组内相似度}$$

具体的计算步骤如下所示（推导过程省略）

（1）首先计算出如下值：

A	B	C	D	Total
$N_A = 5$	$N_B = 5$	$N_C = 5$	$N_D = 5$	$N_T = 20$
$\sum X_{Ai} = 144.30$	$\sum X_{Bi} = 125.20$	$\sum X_{Ci} = 112.50$	$\sum X_{Di} = 111.50$	$\sum X_{Ti} = 493.50$
$\sum X^2_{Ai} = 4196.57$	$\sum X^2_{Bi} = 3158.50$	$\sum X^2_{Ci} = 2576.51$	$\sum X^2_{Di} = 2501.95$	$\sum X^2_{Ti} = 12433.53$
$SS_A = 32.07$	$SS_B = 23.49$	$SS_C = 45.26$	$SS_D = 15.50$	$SS_T = 256.42$

其中：

$$SS = \sum X^2_i - \frac{(\sum X_i)^2}{N}$$

（2）计算$SS_{wg}$以及$SS_{bg}$

$$SS_{wg} = SS_A + SS_B + SS_C + SS_D$$

$$SS_{bg} = SS_T - SS_{wg}$$

（4）计算相关自由度

$$df_{bg} = k - 1 = 4 - 1 = 3$$

$$df_{wg} = (N_A - 1) + (N_B - 1) + (N_C - 1) + (N_D - 1)$$

（5）计算 $MS_{bg}$以及$MS_{wg}$

$$MS_{bg} = \frac{SS_{bg}}{df_{bg}}$$

$$MS_{wg} = \frac{SS_{wg}}{df_{wg}}$$

（6）计算$F$

最后得出F = 6.42 （df = 3, 16）

代码实现

import os
import xlrd
import numpy as np
import scipy.stats as stats
def GetChildDetailLab(combine=True):
    """ Get the detail labs of each child data.
    Returns:
        child_detail_lab: the detail labs.
    """
    CHILDS_FEATURE_LAB = ['tid', 'age_m', 'female', 'onlychild', 'divorce',
                          'medu_newcat', 'income_newcat', 'scr_ave', 'edu_ave',
                          'edu_of_scr', 'scr_h_cat', 'mediacoview',
                          'mediacontact']
    CHILDS_CAT_LAB = ['emo_cat', 'con_cat', 'hyp_cat', 'pee_cat',
                      'difficulties_cat', 'pro_cat']
    if combine:
        CHILDS_DETAIL_LAB = CHILDS_FEATURE_LAB
        CHILDS_DETAIL_LAB.extend(CHILDS_CAT_LAB)
        return CHILDS_DETAIL_LAB
    else:
        return CHILDS_FEATURE_LAB, CHILDS_CAT_LAB
def ReadChildInfoFromExcel(
        file_name='屏幕暴露与SDQ.xlsx', sheet_name='data'):
    """ Read the screen-exposed vs.SDQ detail information
        of each child from the excel-type file.
    Arguments:
        file_name: the name of the excel-type file.
        sheet_name: the name of the sheet of the excel file.
    Returns:
        child_scr_exp_sdq: A list that contains the detail
                           information of each child.
        labs: The lab corresponds to each colume of the data.
    """
    CHILDS_FILE_NAME = 'child_scr_exp_sdq.npy'
    CHILDS_DETAIL_LAB = GetChildDetailLab()
    # print(CHILDS_DETAIL_LAB)
    NOT_INT_LAB_INDEIES = [CHILDS_DETAIL_LAB.index('age_m'),
                           CHILDS_DETAIL_LAB.index('scr_ave'),
                           CHILDS_DETAIL_LAB.index('edu_ave'),
                           CHILDS_DETAIL_LAB.index('edu_of_scr')]
    child_scr_exp_sdq = []
    if (os.path.isfile(CHILDS_FILE_NAME)):
        with open(CHILDS_FILE_NAME, 'rb') as f:
            child_scr_exp_sdq = np.load(f)
    else:
        workBook = xlrd.open_workbook(file_name)
        bookSheet = workBook.sheet_by_name(sheet_name)
        # read from second row because of the first row has tabs
        for row in range(1, bookSheet.nrows):
            child_row = []
            for col in range(bookSheet.ncols):
                cel = bookSheet.cell(row, col)
                try:
                    val = str(cel.value)
                    # tolerant the value error
                    if val == '':
                        val = '-1.0'
                except Exception as e:
                    print(e)
                # because of the type is different
                try:
                    if col in NOT_INT_LAB_INDEIES:
                        val = float(val)
                    else:
                        # in Excel, if cel.value is 1, then str(cel.value) is
                        # '1.0'
                        val = val.split('.')[0]
                        val = int(val)
                except Exception as e:
                    print(e)
                    val = -1
                child_row.append(val)
            child_scr_exp_sdq.append(child_row)
        child_scr_exp_sdq = np.array(child_scr_exp_sdq)
        with open(CHILDS_FILE_NAME, 'wb') as f:
            np.save(f, child_scr_exp_sdq)
    return child_scr_exp_sdq, CHILDS_DETAIL_LAB
def SplitDataSet(data_set, feature, cat):
    """ Split the data set to two column, that is feature and cat
    Arguments:
        data_set: the source data set.
        feature: the input vector
        cat: the correspond category.
    Returns:
        splited_data_set: self-explation.
    """
    CHILDS_DETAIL_LAB = GetChildDetailLab()
    feature_index = CHILDS_DETAIL_LAB.index(feature)
    cat_index = CHILDS_DETAIL_LAB.index(cat)
    # print(feature_index, cat_index)
    return data_set[:, (feature_index, cat_index)]
def CalChi2(data_set):
    """ Calculate the chi-square value and p-value corresponds to the data set.
    Arguments:
        data_set: the object data set.
    Returns:
        chi2: the chi-square value.
        p: the p-value.
    """
    rows_number = len(set(data_set[:, -1]))
    columns_number = len(set(data_set[:, 0]))
    # print(rows_number, columns_number)
    counts = np.zeros((rows_number, columns_number))
    for row in data_set:
        if row[-1] != -1:
            if row[0] != -1:
                try:
                    counts[int(row[-1])][int(row[0])] += 1
                except:
                    pass
    # drop the row that all item is 0
    del_row_index = []
    for index, count in enumerate(counts):
        if not count.any():
            del_row_index.append(index)
    counts = np.delete(counts, tuple(del_row_index), axis=0)
    # drop the column that all item is 0
    del_col_index = []
    for index, count in enumerate(counts.T):
        if not count.any():
            del_col_index.append(index)
    counts = np.delete(counts, tuple(del_col_index), axis=1)
    # print(counts)
    # calculate the chi-square value and correspond p-value
    chi2, p, dof, excepted = stats.chi2_contingency(counts)
    return chi2, p
def ANOVATest(data_set):
    """ Implement a ANOVA test on the data set.
    Arguments:
        data_set: the object data set.
    Return:
        f: The computed F-value of the test.
        p: The associated p-value from the F-distribution.
    """
    # Initial the three categories
    normal = []
    critical = []
    abnormal = []
    for data in data_set:
        if data[0] != -1:
            if data[-1] == 0:
                normal.append(data[0])
            elif data[-1] == 1:
                critical.append(data[0])
            elif data[-1] == 2:
                abnormal.append(data[0])
    f, p = stats.f_oneway(normal, critical, abnormal)
    return f, p
def GenerateCoffMatrix(data_set):
    """ Calculate the chi-square and p-value of each feature-category pair.
    Arguments:
        data_set: the source data set.
    Returns:
        coff_matrix: the final cofficient matrix.
    """
    coff_matrix = {}
    CHILDS_FEATURE_LAB, CHILDS_CAT_LAB = GetChildDetailLab(combine=False)
    NOT_KEEP_FEATURE_LAB = ['tid', 'age_m', 'scr_ave', 'edu_ave', 'edu_of_scr']
    for feature in CHILDS_FEATURE_LAB:
        if feature in NOT_KEEP_FEATURE_LAB:
            if feature != NOT_KEEP_FEATURE_LAB[0]:
                for cat in CHILDS_CAT_LAB:
                    splited_data_set = SplitDataSet(data_set, feature, cat)
                    # print(feature, cat)
                    f, p = ANOVATest(splited_data_set)
                    key = feature + '-' + cat
                    coff_matrix[key] = (f, p)
        else:
            for cat in CHILDS_CAT_LAB:
                splited_data_set = SplitDataSet(data_set, feature, cat)
                # print(feature, cat)
                chi2, p = CalChi2(splited_data_set)
                key = feature + '-' + cat
                coff_matrix[key] = (chi2, p)
    return coff_matrix
def SiftRelativeFeature(coff_matrix, conf=1e-5):
    """ Sift the feature that satisfy the caonfident condition.
    Arguemnts:
        coff_matrix: the calculated cofficient matrix for
                     all features and categories.
        conf: the confident.
    Returns:
        relative_feature_matrix: the satisfied feature and correspond
                                 category and chi2 and p-value.
    """
    relative_feature_matrix = {}
    for key in coff_matrix.keys():
        if coff_matrix[key][-1] <= conf:
            relative_feature_matrix[key] = coff_matrix[key]
    return relative_feature_matrix
def WriteResult(coff_matrix, file_name='result.txt'):
    """ Write the result to file.
    Arguments:
        relative_feature_matrix: the satisfied feature and correspond
                                 category and chi2 and p-value.
        file_name: the result file name.
    """
    # Sort
    sorted_coff_matrix = sorted(coff_matrix.items(),
                                key=lambda item: item[1][-1], reverse=False)
    # print(sorted_coff_matrix)
    with open(file_name, 'w') as f:
        for item in sorted_coff_matrix:
            f.write(str(item))
            f.write('\n')
if __name__ == '__main__':
    data_set, labels = ReadChildInfoFromExcel()
    # splited_data_set = SplitDataSet(data_set, 'female', 'difficulties_cat')
    # chi2, p = CalChi2(splited_data_set)
    coff_matrix = GenerateCoffMatrix(data_set)
    # relative_feature = SiftRelativeFeature(coff_matrix, 1)
    WriteResult(coff_matrix)

因子分析

先说说因子分析与主成分分析的区别，下面的话引自知乎(https://www.zhihu.com/question/24524693)

具体实现使用SPSS软件进行实现：

操作步骤如下(引自 SPSS数据分析从入门到精通-陈胜可)：