SAS Programming for Data Mining

Visualize your own image filter

2018-11-23T10:39:00.000-08:00

First load necessary packages


import pandas as pd
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Dropout, Conv2D
import keras.backend as K
import scipy, imageio
import matplotlib.pyplot as plt
from PIL import Image
%matplotlib inline

Then show original picture of my Jeep


# 首先将图片读入为矩阵
# 我们可以用pyplot的imshow()方法来展示图片
# 这是我曾经拥有的牧马人JK Rubicon Unlimited
#
img_data = imageio.imread('./pics/wranglerJK.jpg')
print(img_data.data.shape)

img = Image.fromarray(img_data, 'RGB')
plt.imshow(img)

添加图片说明

Now, build our 2-D convolutional function that takes a custom filter matrix and comput the filtered output image matrix.

def my_init(shape, dtype=None):
    new_mat = np.zeros((shape[0], shape[1], 3, 3))
    for i in range(shape[0]):
        for j in range(shape[1]):
            new_mat[:, :, i, j] = filter_mat
    return np.array(new_mat, dtype=dtype)


def MyFilter(filter_mat):
    print(len(filter_mat.shape))
    if len(filter_mat.shape)!=2:
        print('Invalid filter matrix. It must be 2-D')
        return []
    else:
        kernel_size=filter_mat.shape
        row, col, depth = img_data.shape
        input_shape=img_data.shape
        filter_size = row*col*depth
        print(filter_size)


        model = Sequential()
        model.add(Conv2D(depth, 
                         kernel_size=kernel_size, 
                         input_shape=input_shape, 
                         padding='same', 
                         activation='linear', 
                         data_format='channels_last',
                         kernel_initializer=my_init,
                         name='Conv')                         
                  )
        model.add(Dense(1, activation='linear'))
        model.compile(optimizer='sgd', loss='mse')
        model.summary()


        inX = model.input                                          
        outputs = [layer.output for layer in model.layers if layer.name=='Conv']     
        functions = [K.function([inX], [out]) for out in outputs]   

        layer_outs = [func([img_data.reshape(1, row, col, depth)]) for func in functions]
        activationLayer = layer_outs[0][0]
        temp = (activationLayer-np.min(activationLayer))
        normalized_activationLayer = temp/np.max( np.max(temp))
        return(normalized_activationLayer.reshape(row, col, depth))

Now, insert our own fixed filter matrix and get the output, using pyplot.imshow() to display the filtered picture. This time we throw in an edge detector.

filter_mat = np.array([-1, -2, -3, 0, 0, 0, 1, 2, 3]).reshape(3, 3)

outLayer = MyFilter(filter_mat)
plt.imshow(outLayer)

Below is the filtered picture.

Visualize activations functions using keras

2017-12-22T00:44:00.002-08:00

In keras, we can visualize activation functions' geometric properties using backend functions over layers of a model.

We all know the exact function of popular activation functions such as 'sigmoid', 'tanh', 'relu', etc, and we can feed data to these functions to directly obtain their output. But how to do that via keras without explicitly specifying their functional forms?

This can be done following the four steps below:

1. define a simple MLP model with a one dimension input data, a one neuron dense network as the hidden layer, and the output layer will have a 'linear' activation function for one neuron.
2. Extract layers' output of the model (fitted or not) via iterating through model.layers
3. Using backend function K.function() to obtain calculated output for a given input data
4. Feed desired data to the above functions to obtain the output from appropriate activation function.

The code below is a demo:




from keras.layers import Dense, Activation
from keras.models import Sequential
import keras.backend as K
import numpy as np
import matplotlib.pyplot as plt


# 以下设置显示中文文方法根据 http://blog.csdn.net/rumswell/article/details/6544377 
plt.rcParams['font.sans-serif'] = ['SimHei'] #指定默认字体  
plt.rcParams['axes.unicode_minus'] = False #解决图像中中文符号显示为方块的问题  

def NNmodel(activationFunc='linear'):
    '''
       定义一个神经网络模型。如果要定义不同的模型，可以直接修改该函数
    '''
    if (activationFunc=='softplus') | (activationFunc=='sigmoid'):
        winit='lecun_uniform'
    elif activationFunc=='hard_sigmoid':
        winit='lecun_normal'
    else:
        winit='he_uniform'
    model = Sequential()
    model.add(Dense(1, input_shape=(1,), activation=activationFunc, 
                    kernel_initializer=winit, 
                    name='Hidden'))
        
    model.add(Dense(1, activation='linear', name='Output'))   
    model.compile(loss='mse', optimizer='sgd')    
    return model

def VisualActivation(activationFunc='relu', plot=True):    
    x = (np.arange(100)-50)/10
    y = np.log(x+x.max()+1)
    
    model = NNmodel(activationFunc = activationFunc)
    
    inX = model.input                                          
    outputs = [layer.output for layer in model.layers if layer.name=='Hidden']     
    functions = [K.function([inX], [out]) for out in outputs]   
    
    layer_outs = [func([x.reshape(-1, 1)]) for func in functions]
    activationLayer = layer_outs[0][0]
        
    activationDf = pd.DataFrame(activationLayer)
    result=pd.concat([pd.DataFrame(x), activationDf], axis=1)
    result.columns=['X', 'Activated']
    result.set_index('X', inplace=True)
    if plot:
        result.plot(title=f)
    
    return result


# Now we can visualize them (assuming default settings) :
actFuncs = ['linear', 'softmax', 'sigmoid', 'tanh', 'softsign', 'hard_sigmoid', 'softplus', 'selu', 'elu']

from keras.layers import LeakyReLU
figure = plt.figure()
for i, f in enumerate(actFuncs):
    # 依次画图
    figure.add_subplot(3, 3, i+1)
    out=VisualActivation(activationFunc=f, plot=False)
    plt.plot(out.index, out.Activated)
    plt.title(u'激活函数：'+f)

This figure is the output from above code. As we can see, the geometric property of each activation function is well captured.

A VGG-like CNN in keras for Fashion-MNIST with 94% accuracy

2017-09-06T17:48:00.002-07:00

In many introductory to image recognition tasks, the famous MNIST data set is typically used. However, there are some issues with this data:

1. It is too easy. For example, a simple MLP model can achieve 99% accuracy, and a 2-layer CNN can achieve 99% accuracy.

2. It is over used. Literally every machine learning introductory article or image recognition task will use this data set as benchmark. But because it is so easy to get nearly perfect classification result, its usefulness is discounted and is not really useful for modern machine learning/AI tasks.

Therefore, there appears Fashion-MNIST dataset. This dataset is developed as a direct replacement for MNIST data in the sense that:

1. It is the same size and style: 28x28 grayscale image
2. Each image is associated with 1 out of 10 classes, which are:
       0:T-shirt/top,
       1:Trouser,
       2:Pullover,
       3:Dress,
       4:Coat,
       5:Sandal,
       6:Shirt,
       7:Sneaker,
       8:Bag,
       9:Ankle boot
3. 60000 training sample and 10000 testing sample Here is a snapshot of some samples:

Since its appearance, there have been multiple submissions to benchmark this data, and some of them are able to achieve 95%+ accuracy, most noticeably Residual network or separable CNN.

I am also trying to benchmark against this data, using keras. keras is a high level framework for building deep learning models, with selection of TensorFlow, Theano and CNTK for backend. It is easy to install and use. For my application, I used CNTK backend. You can refer to this article on its installation.

Here, I will benchmark two models. One is a MLP with layer structure of 256-512-100-10, and the other one is a VGG-like CNN. Code is available at my github: https://github.com/xieliaing/keras-practice/tree/master/fashion-mnist

The first model achieved accuracy of [0.89, 0.90] on testing data after 100 epochs, while the latter achieved accuracy of >0.94 on testing data after 45 epochs. First, read in the Fashion-MNIST data:


import numpy as np
import io, gzip, requests
train_image_url = "http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/train-images-idx3-ubyte.gz"
train_label_url = "http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/train-labels-idx1-ubyte.gz"
test_image_url = "http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/t10k-images-idx3-ubyte.gz"
test_label_url = "http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/t10k-labels-idx1-ubyte.gz"

def readRemoteGZipFile(url, isLabel=True):
    response=requests.get(url, stream=True)
    gzip_content = response.content
    fObj = io.BytesIO(gzip_content)
    content = gzip.GzipFile(fileobj=fObj).read()
    if isLabel:
        offset=8
    else:
        offset=16
    result = np.frombuffer(content, dtype=np.uint8, offset=offset)    
    return(result)

train_labels = readRemoteGZipFile(train_label_url, isLabel=True)
train_images_raw = readRemoteGZipFile(train_image_url, isLabel=False)

test_labels = readRemoteGZipFile(test_label_url, isLabel=True)
test_images_raw = readRemoteGZipFile(test_image_url, isLabel=False)

train_images = train_images_raw.reshape(len(train_labels), 784)
test_images = test_images_raw.reshape(len(test_labels), 784)

Let's first visual it using tSNE. tSNE is said to be the most effective dimension reduction tool.This plot function is borrowed from sklearn example.


from sklearn import manifold
from time import time
import matplotlib.pyplot as plt
from matplotlib import offsetbox
plt.rcParams['figure.figsize']=(20, 10)
# Scale and visualize the embedding vectors
def plot_embedding(X, Image, Y, title=None):
    x_min, x_max = np.min(X, 0), np.max(X, 0)
    X = (X - x_min) / (x_max - x_min)

    plt.figure()
    ax = plt.subplot(111)
    for i in range(X.shape[0]):
        plt.text(X[i, 0], X[i, 1], str(Y[i]),
                 color=plt.cm.Set1(Y[i] / 10.),
                 fontdict={'weight': 'bold', 'size': 9})

    if hasattr(offsetbox, 'AnnotationBbox'):
        # only print thumbnails with matplotlib > 1.0
        shown_images = np.array([[1., 1.]])  # just something big
        for i in range(X.shape[0]):
            dist = np.sum((X[i] - shown_images) ** 2, 1)
            if np.min(dist) < 4e-3:
                # don't show points that are too close
                continue
            shown_images = np.r_[shown_images, [X[i]]]
            imagebox = offsetbox.AnnotationBbox(
                offsetbox.OffsetImage(Image[i], cmap=plt.cm.gray_r),
                X[i])
            ax.add_artist(imagebox)
    plt.xticks([]), plt.yticks([])
    if title is not None:
        plt.title(title)

tSNE is very computationally expensive, so for impatient people like me, I used 1000 samples for a quick run. If your PC is fast enough and have time, you can run tSNE against the full dataset.


sampleSize=1000
samples=np.random.choice(range(len(Y_train)), size=sampleSize)
tsne = manifold.TSNE(n_components=2, init='pca', random_state=0)
t0 = time()
sample_images = train_images[samples]
sample_targets = train_labels[samples]
X_tsne = tsne.fit_transform(sample_images)
t1 = time()
plot_embedding(X_tsne, sample_images.reshape(sample_targets.shape[0], 28, 28), sample_targets,
               "t-SNE embedding of the digits (time %.2fs)" %
               (t1 - t0))
plt.show()

We see that several features, including mass size, split on bottom and semetricity, etc, separate the categories. Deep learning excels here because you don't have to manually engineering the features but let the algorithm extracts those.

In order to build your own networks, we first import some libraries


from keras.models import Sequential
from keras.layers.convolutional import Conv2D, MaxPooling2D, AveragePooling2D
from keras.layers.advanced_activations import LeakyReLU
from keras.layers import Activation

We also do standard data preprocessing:


X_train = train_images.reshape(train_images.shape[0], 28, 28, 1).astype('float32')
X_test = test_images.reshape(test_images.shape[0], 28, 28, 1).astype('float32')

X_train /= 255
X_test /= 255

X_train -= 0.5
X_test -= 0.5

X_train *= 2.
X_test *= 2.

Y_train = train_labels
Y_test = test_labels
Y_train2 = keras.utils.to_categorical(Y_train).astype('float32')
Y_test2 = keras.utils.to_categorical(Y_test).astype('float32')

Here is the simple MLP implemented in keras:


mlp = Sequential()
mlp.add(Dense(256, input_shape=(784,)))
mlp.add(LeakyReLU())
mlp.add(Dropout(0.4))
mlp.add(Dense(512))
mlp.add(LeakyReLU())
mlp.add(Dropout(0.4))
mlp.add(Dense(100))
mlp.add(LeakyReLU())           
mlp.add(Dropout(0.5))
mlp.add(Dense(10, activation='softmax'))
mlp.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
mlp.summary()

This model achieved almost 90% accuracy on test dataset at about 100 epochs. Now, let's build a VGG-like CNN model. We use an architecture that is similar to VGG but still very different. Because the figure data is small, if we use original VGG architecture, it is very likely to overfit and won't perform very well in testing data which is observed in publically submitted benchmarks listed above. To build such a model in keras is very natural and easy:


num_classes = len(set(Y_train))
model3=Sequential()
model3.add(Conv2D(filters=32, kernel_size=(3, 3), padding="same", 
          input_shape=X_train.shape[1:], activation='relu'))
model3.add(Conv2D(filters=64, kernel_size=(3, 3), padding="same", activation='relu'))
model3.add(MaxPooling2D(pool_size=(2, 2)))
model3.add(Dropout(0.5))
model3.add(Conv2D(filters=128, kernel_size=(3, 3), padding="same", activation='relu'))
model3.add(Conv2D(filters=256, kernel_size=(3, 3), padding="valid", activation='relu'))
model3.add(MaxPooling2D(pool_size=(3, 3)))
model3.add(Dropout(0.5))
model3.add(Flatten())
model3.add(Dense(256))
model3.add(LeakyReLU())
model3.add(Dropout(0.5))
model3.add(Dense(256))
model3.add(LeakyReLU())
#model2.add(Dropout(0.5))
model3.add(Dense(num_classes, activation='softmax'))
model3.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
model3.summary()

This model has 1.5million parameters. We can call 'fit' method to train the model:


model3_fit=model3.fit(X_train, Y_train2, validation_data = (X_test, Y_test2), epochs=50, verbose=1, batch_size=500)

After 40 epochs, this model archieves accuracy of 0.94 on testing data.Obviously, there is also overfitting problem for this model. We will address this issue later.

如何零基础使用KERAS搭建实用深度模型

2017-09-03T00:32:00.001-07:00

在这篇小文章中，我们将简要讨论如何使用KERAS这个现在最新的深度学习框架来构造实用的深度学习模型。

深度学习是目前最热门的高级分析技术之一，在很多方面表现出了超越传统机器学习方法的有效性。但是在常用的TensorFlow，CNTK，Theano等计算环境中实现不同的深度学习模型仍然需要耗费很多时间来编写程序。KERAS的出现提供了一个高度抽象的环境来搭建深度学习模型，特别是其简单易用，跟网络结构一一对应的特点使得其迅速在数据科学家这个使用人群中流行起来。

## 什么是KERAS

KEARS是Google工程师François Chollet为主创人员，基于Python开发和维护的一个抽象的神经网络建模环境，提供了一系列的API供用户调用构造自己的深度学习网络。KERAS的出发点就是为用户提供一个能够快速实现模型的手段，从而缩短建模迭代的时间，加快模型试验的频率。用KERAS开发者的话说，就是要做好的科研必须尽可能地缩短从想法到实现结果的时间。在业界工作中这也是成功的关键要素之一。

相比较于常见的深度学习环境，比如TensorFlow，CNTK，Theano，Caffe等，KERAS有以下几个不同：

1. 设计初衷就是方便以模块化地方式快速构造深度学习模型的原型；

2. 可以很方便地在CPU和GPU之间切换

3. KERAS本身只是描述模型的环境，其计算平台目前依赖于TensorFlow，CNTK和Theano这三种，以后会拓展到其他流行的计算平台上，比如mxNet等；

4. KERAS的拓展性既可以通过自定义KERAS里的激活函数或者损失函数等能自定义的部分进行，也可以通过引用对应的计算平台的自定义部分进行，具有一定的灵活性；

跟这些流行的计算平台一样，KERAS也支持常见的深度学习模型，比如卷积神经网络，循环神经网络以及二者的组合等。

使用KERAS构造深度神经网络有一系列相对固定的步骤：

1. 首先要将原始数据处理成KERAS的API能够接受的格式，一般是一个张量的形式，通常在维度上表示为（批量数，[单一样本对应张量的维度]）。这里[单一样本对应张量的维度] 是一个通用的说法，对应于不同类型的模型，数据有不同的要求。通常，如果是一个简单的全链接模型，则单一样本对应张量的维度就是特征个数；如果是一维的时间序列数据，并要用循环神经网络模型训练的话，则单一样本对应张量的维度是时间步和每个时间步对应的回看序列长度；如果输入数据是图像，并使用卷积神经网络模型进行训练，则单一样本张量对应图像的高，宽和色彩频道三个维度。但是如果是使用全连接模型训练图像数据，则单一样本对应张量是该图像扁化（Flatten）以后的向量长度，其为高，宽和色彩频道各个维度数量的乘积。一般卷积神经网络最靠近输出层的那层都设置一个全连接层，因此也需要扁化输入张量。

2. 其次要构造需要的深度学习模型。这一步又分为模型的选择和模型的细化两个步骤：
   - 选择模型的类型。KERAS里定义了两大类模型
             1）序列模型(Sequential)；
             2）通用模型(Model)。

            序列模型指的是深度模型每一层之间都是前后序列关系，如下图所示：

Figure 1。MLP是一个典型的序列模型，[图片来源](http://article.sapub.org/10.5923.j.ajis.20120204.01.html) 可以看到从左到右，输入层到隐含层到输出层每一层之间都是前后依次相连的简单关系。这个简单的网络结构可以用三句KERAS命令实现：


    model=Sequential()
    model.add(Dense(5, input_shape=(4,), activation=’sigmoid’))
    model.add(Dense(1, activation=’sigmoid’))

而通用模型则是对应更广义的模型，具备更大的灵活性。上面提到的序列模型也可以用通用模型来表达，这个我们在后一节详解。当然通用模型更能用来描述层与层之间有较复杂关系的情况，比如非相邻的层之间进行连接，或者多个神经网络的合并等。比如我们可以使用通用模型进行矩阵分解：


    user_in = Input(shape=(1,), dtype='int64', name='user_in')
    u = Embedding(n_users, n_factors, input_length=1)(user_in)
    movie_in = Input(shape=(1,), dtype='int64', name='movie_in')
    v = Embedding(n_movies, n_factors, input_length=1)(movie_in)
    x = merge([u, v], mode='dot')
    x = Flatten()(x)
    model = Model([user_in, movie_in], x)
    model.compile(Adam(0.001), loss='mse')

这里构造了一个基于矩阵分解的推荐系统的一个深度学习模型，其对应的网络结构如下图所示：

Figure 2。矩阵分解的深度学习模型

- 细化模型的结构。其实上面已经展示了模型结构细化之后的情况。一般来说，确定了模型类型以后，其结构不外乎每层的类型是什么，是全连接层还是卷积层还是放弃（Dropout）层；每层的其他参数是什么，比如如果需要指定激活函数，那么使用什么样的激活函数，如果是卷积层，那么需要多少过滤器，每个过滤器的大小是怎样的？等等这些都可以通过设定不同的参数进行细化。

3. 然后对模型进行编译，编译完成以后可以查看模型的基本信息，特别是参数的数量；

4. 最后带入数据对模型进行拟合。一般来讲，如果数据是静态的张量数据，通过使用fit方法。如果数据特别大，可是使用可迭代的data generator对象，并使用fit_generator方法来拟合。

## KERAS和深度学习模型的对应关系

KERAS既然是开发出来快速构造深度学习模型的工具，那么它的API和深度学习模型的要素都有很强的对应关系。正如上面所说，目前的深度学习模型都可以纳入序列模型或者通用模型的，那么我们用图示的方式来表示这个对应关系，方便读者理解。这里网络图为了方便与按行排列的代码对应，对每一层都进行了标注。下图展示的是一个典型的全连接序列模型：

Figure 3。全连接序列模型，修改自[博客](http://www.samyzaf.com/ML/pima/pima.html)

这个序列模型可以使用如下的KERAS命令快速搭建：


Model = Sequential()
Model.add(Dense(10, activation=’sigmoid’, 
                   input_shape=(8, ))             【隐含层1+输入层】
Model.add(Dense(8, activation=’relu’))          【隐含层2】
Model.add(Dense(10, activation=’relu’))         【隐含层3】
Model.add(Dense(5, activation=’softmax’))       【输出层】

上面的序列模型也可以用通用模型的API描述的结果，其与图中的网络结构有更强的对应关系：


x = Input(shape=(8,))                          【输入层】
b = Dense(10, activation=’sigmoid’)(x)       【隐含层1】
c = Dense(8, activation=’relu’)(b)           【隐含层2】
d = Dense(10, activation=’relu’)(c )         【隐含层3】
out = Dense(5, activation=’softmax’)(d)      【输出层】
model = Model(inputs=x, outputs=out)

上面也举了另外的比较复杂的例子。在后面的具体案例中，我们也会强调网络结构和对应的KERAS命令，使读者能建立起较强的联系。 ##使用KERAS构造深度推荐系统推荐系统是机器学习最广泛的应用领域之一，大家熟悉的亚马逊、迪士尼、谷歌、Netflix 等公司都在网页上有其推荐系统的界面，帮助用户更快、更方便地从海量信息中找到有价值的信息。比如亚马逊（www.amazon.com）会给你推荐书、音乐等，迪士尼（video.disney.com）给你推荐最喜欢的卡通人物和迪士尼电影，谷歌搜索更不用说了， Google Play、 Youtube 等也有自己的推荐引擎、推荐视频和应用等。下面是我登陆亚马逊之后的一个推荐页面，可见我之前应该是购买了咖啡机，所以会有相关的产品推荐出来。

Figure 4。亚马逊的推荐页面局部

推荐系统的最终目的是从百万甚至上亿内容或者商品中把有用的东西高效地显示给用户，这样可以为用户节省很多自行查询的时间，也可以提示用户可能忽略的内容或商品，使用户更有黏性，更愿意花时间待在网站上，从而使商家可以从内容或者商品中赚取更多的利润，即使流量本身也会使商家从广告中受益。

传统上，推荐系统是基于矩阵分解的协同过滤算法，前面也展示了这样的一个简单模型。下面我们着重介绍深度学习推荐系统。这个模型除了能将用户和可选产品联系起来意外，还能将其他辅助数据，比如用户年龄，地区，上网设备以及各种产品属性，联系起来。

这里通过嵌入（Embedding）这种技术将不同的信息串在一起作为输入层，再继续搭建不同的神经网络模型，最后一层用预测评分作为输出层。虽然这里的数据只有用户编码和电影产品编码，但是这样的结构可以拓展到包含其他相关数据。下图展示了这样的一个深度模型的结构示意图：

Figure 5。深度模型

有了这个示意图，我们就可以很方便地用KERAS依次构造。这里我们假设已经将用户和电影产品做了按照One Hot编码形式组织好了。

首先用嵌入层对用户和电影进行嵌入映射：


k = 128
model1 = Sequential()
model1.add(Embedding(n_users + 1, k, input_length = 1))
model1.add(Reshape((k,)))
model2 = Sequential()
model2.add(Embedding(n_movies + 1, k, input_length = 1))
model2.add(Reshape((k,)))

这里的k是映射到的空间的维度。在一般的业务系统中我们可能有上百万的用户和产品，经过嵌入映射到128维的实数域上以后显著地降低了整个系统的维度和大小。以上几句命令实现了上图从最低下到“用户嵌入”和“电影嵌入”这一阶段的编程。

其次，我们需要用第三个神经网络把前面的两个嵌入网络映射所得到的向量叠加在一起：


model = Sequential()
model.add(Merge([model1, model2], mode = 'concat'))

至此完成了到第一个粗箭头的网络构造。两个网络已经合并为一个网络。下面的命令依次完成“隐含层128”和“隐含层32”的构造：


model.add(Dropout(0.2))
model.add(Dense(k, activation = 'relu'))
model.add(Dropout(0.5))
model.add(Dense(int(k/4), activation = 'relu'))
model.add(Dropout(0.5))

下面继续构造“隐含层8”：


model.add(Dense(int(k/16), activation = 'relu'))
model.add(Dropout(0.5))

隐含层构造完毕之后，需要构造输出层。因为是预测连续变量评分，最后一层直接上线性变化：


model.add(Dense(1, activation = 'linear'))

至此，模型构造完毕，可以编译了：


model.compile(loss = 'mse', optimizer = "adam")

这里使用了均方差（MSE）作为损失函数，并使用了ADAM优化算法。下面，为了能训练模型，需要将数据构造为[users, movies]的形式：


users = ratings['user_id'].values
movies = ratings['movie_id'].values
X_train = [users, movies]

最后训练模型：


model.fit(X_train, y_train, batch_size = 100, epochs = 50)

使用movielens的用户观看电影评分数据进行训练和验证，我们发现这个模型的误差在0.8226左右，大约一个评分等级不到。

即使这样一个简单的模型，效果还是比较好的。如果进一步优化结构，或者引入其他信息，误差还可以进一步降低。

##使用KERAS构造图像识别系统

图像识别是深度学习最典型的应用之一。关于深度学习的图像识别可以追溯很长的历史，其中最具有代表性的例子是手写字体识别和图片识别。手写字体识别主要是用机器正确区别手写体数字 0~9。银行支票上的手写体识别技术就是基于这个技术。图片识别的代表作就是 ImageNet。这个比赛需要团队识别图片中的动物或者物体，把它们正确地分到一千个类别中的其中一个。图像识别有很多种技术可以实现，目前最主流的技术是深度神经网络，其中尤以卷积神经网络（CNN）最为出名。卷积神经网络（见图1）是一种自动化特征提取的机器学习模型。从数学的角度看，任何一张图片都可以对应到 224 × 224 × 3 或者 32 × 32 × 3 等三维向量，这取决于像素。我们的目标是把这个三维向量（又被称为张量）映射到 N个类别中的一类。神经网络就是建立了这样一个映射关系，或者称为函数。它通过建立网状结构，辅以矩阵的加、乘等运算，最后输出每个图像属于每个类别的概率，并且取概率最高的作为我们的决策依据。下面是一个典型的序列卷积神经网络模型的结构：

Figure 6。卷积神经网络结构。来源于CNTK教程

上面这个网络依次展示了卷积网络模型的主要要素：
- 输入层的图像；
- 卷积操作；
- 激活函数的应用；
- 池化操作；
- 将数据展平（Flatten），为输出到全连接层做准备；
- 全连接层准备输出；
- softmax应用于分类问题的全连接层作为输出层。

下面详细介绍一下在KERAS中如何对应地进行编程处理。

- 首先，这是一个序列模型，因此先要声明一个序列模型的对象：
model.=sequential()

- 卷积是应用一个局部的过滤器到原始数据的过程，比如下图就展示了一个3x3的过滤器应用在一个7x7的图像上过程。假设在当前步，这个过滤器的权重经过学习得到如图所示的权重结果，在当前步继续进行卷积操作就是将这个3x3的过滤器从左上角每次要么向右要么向下移动一格，和对应的图像局部的3x3的区域进行依元素点乘求和得到卷积的结果。因为依次移动，到最边上的时候过滤器会超出图像的边界，一般会将这些对应的卷积结果删除，从而卷积后的张量维度会少于原始图像。比如这个例子中原图为7x7，使用一个3x3的过滤器卷积之后最后两列和两行进行卷积的时候会使过滤器超过边界，因此最后的卷积结果是一个5x5的图像。这里可以使用多个过滤器，每个过滤器应用一次，每次应用产生的卷积结果构成隐含层的一层。比如采用16个过滤器，如果不删除边界的过滤结果，则得到新的[7x7x16]的张量。

Figure 7．卷积演示，来源于CNTK教程

在KERAS里，对于图像这种二维数据，一般使用Conv2D这个二维卷积层。Conv2D有几个必备的参数：
1. 首先指定过滤器数量 filter，是一个整数；
2. 第二是要指定二维过滤器的大小，比如(3,3)；
3. 第三要指定步长，即延某轴移动的时候是依次移动一个像素还是多个像素，默认为1；
4. 第四要指定补齐策略padding，即是否要将在边界的卷积结果去掉。如果值为”same”，则不去掉，卷积结果和输入图像有同样的高和宽；如果值为”valid”，则不会处理过滤器会超出边界的像素。
5. 最后，如果卷积层是第一层，那么还需要指明输入数据的维度input\_shape。因为一般用TensorFlow或者CNTK做后台，输入数据要求是channel_last，因此输入的原始维度为[样本量，高，宽，频道]，那么这里的维度只需要去掉样本量即可，即为[高，宽，频道数]，一般用X.shape[1:]即可得到。

对于上面的例子，KERAS里的典型写法是：


model.add(Conv2D(filters=16, kernel_size=(3, 3),
                  strides=1, padding=”valid”, 
                  input_shape=xtrain.shape[1:]))

- 再次要添加激活层引入激活函数，通常是一个非线性的函数。激活函数既可以通过在Conv2D里面指定activation=参数引入，也可以通过单独添加Activation层引入。

卷积神经网络常用的激活函数是Rectified Linear Unit，简称relu。该函数其实就是max(0, x)，在层次较深的网络中比以前常用的取值区间在（0，1）或者（-1， 1）之间的sigmoid类激活函数效果好，因为不存在梯度消失的问题。

如果是通过参数，则上面的代码改写为：


model.add(Conv2D(filters=16, kernel\_size=(3, 3),
                 strides=1, padding=”valid”, 
                 activation=’relu’)

如果通过添加激活层引入，则在上面的代码后添加：


model.add(Activation(‘relu’))

- 然后进行的池化操作是在卷积神经网络中对图像特征的一种处理，通常在卷积操作和激活函数之后进行。池化操作是将原有输入按照一定大小切分成互不交叉的局部区域，目的是为了计算特征在局部的充分统计量，从而降低总体的特征数量，防止过度拟合和减少计算量。下图展示了最大池化方法的应用。在一个6x6的图像上应用3x3的池化操作，将原输入矩阵切割为不相交叉的2x2区域，每个区域的取值是对应原输入局部的最大值。

Figure 8。最大池化操作对应于图像的最大池化层通过MaxPooling2D，KERAS也支持平均池化层，区别在于取对应局部的平均值作为池化后结果，方法为AveragePooling2D。对应上面的例子，KERAS的命令如下：


model.add(MaxPooling2D(pool_size=(3, 3))

- 为了输出到全连接层，先要对数据进行展平（Flatten）。这是因为全连接层只处理包含样本数在内一共二维的数据，要求第一维是样本数，第二维是所有特征的个数。因此对于一个包含2000个样本，每个样本是28x28x3的小图像的数据，展平之后是一个2000x2352的矩阵，其中2352是28,28,3的乘积。在KERAS里进行展平非常简单，在上面的MaxPooling2D层之后添**model.add(Flatten())** 即可，KERAS能自己分析出输入和输出的维度。

- 在前面这些处理之后，但是在输出之前，通常会添加一个或者多个全连接层进一步处理数据。全连接层可以通过Dense指出，指明输出神经元个数和激活函数即可：
model.add(Dense(1000, activation=’relu’))

- 最后使用一个全连接层作为输出层，同样要求使用softmax激活函数，并使用跟输出类别同样多的神经元个数。比如识别0—9十个数字，那么就应该写作：
model.add(Dense(10, activation=’relu’))

把所有步骤组合到一起，我们就可以将图6显示的一个卷积神经网络模型相应地写为KERAS代码了：


model=Sequential()
model.add(Conv2D(filters=32, kernel_size=(3, 3), 
          padding="same", 
          input_shape=X_train.shape[1:], 
          activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Conv2D(filters=64, kernel_size=(3, 3), padding="valid"))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dense(num_classes, activation='softmax'))

是不是很简单？要训练这个模型非常简单。我们先编译这个模型并显示其关键信息：


model.compile(loss='categorical_crossentropy', optimizer='adagrad', metrics=['accuracy'])
model.summary()

Figure 9。模型信息我们看到这个模型一共有421642个参数，大多来自于倒数第二层的全连接层。拟合这个模型也很简单：


model.fit(X_train, y_train, 
          epochs=20, verbose=1, 
          batch_size=10, 
          validation_data = (X_test, y_test))

这里使用最标准的fit方法。其中指定几个核心参数：

1. 训练用特征数据X\_train，
2. 训练用结果数据y\_train，
3. 迭代次数epochs，
4. 批量大小用batch\_size指定，
5. verbose表示显示训练过程的信息，如果值为0不显示任何中间信息，如果为1显示按批量拟合的进度，如果为2则显示拟合结果信息。
6. 另外读者还可以指定验证数据集，用validation_data这个参数表示，其包含一个tuple，第一个元素是验证用特征数据，第二个是验证用结果数据。

下面我们使用这个模型训练识别0-9这十个数字，使用著名的MNIST数据。不过在训练之前还需要提及对数据的处理：

* 首先将数据重塑为[样本数，高，宽，色彩通道数]的格式。这个可以通过numpy.reshape方法来实现。因为keras自带的MNIST数据已经是numpy的多维矩阵，并且是单色的，因此色彩通道数为1，因此数据重塑可以用下面的命令实现。读者可自行重塑验证用数据。


X_train = X_train.reshape(X_train.shape[0], 
                          X_train.shape[1], 
                          X_train.shape[2], 1).astype(float)

* 其次，需要将数据的取值压缩到[0， 1]之间。这有利于拟合时用的随机梯度递降算法的稳定和收敛。这可以使用X_train /= 255.0 来实现。。

* 最后要将结果数据变为One Hot Code形式。KERAS提供了一个非常方便的方法to_categorical来实现这个功能：


y_train = keras.utils.to_categorical(y_train, len(set(y_train)))

下图的结果显示即使是这个非常简单的模型，其在验证数据上的预测准确率都是非常高的，达到了99.14%。

Figure 10。简单卷积模型在MNIST数据上的拟合结果。

使用KERAS可以非常方便的构造自己的卷积神经网络，对于比较复杂的情况，也可以使用已经训练好的一些常见的高效模型，比如VGG16，Xception 等做迁移训练来拟合自己的数据。

Figure 11． VGG16结构，来源于：https://www.cs.toronto.edu/~frossard/post/vgg16

上图是著名的VGG16模型的结构。根据刚才的学习结果，读者可以很快地模仿这个结构搭建自己的类似模型，但是KERAS在application库里已经提供了现成训练好的VGG16模型供读者读入使用。读者可以引用这个模型，将顶层去掉用自己的数据重新训练，但是底层的参数借用VGG16已经训练好的权重。这就是迁移学习的思路，可以大大降低需要训练的参数数量，加快新模型的开发。这里使用了通用模型以便在现有的VGG16模型上进行修改：


model_vgg = VGG16(include_top = False, 
                  weights = 'imagenet', 
                  input_shape =(224,224,3))
model = Flatten(name = 'flatten')(model_vgg.output)
model = Dense(10, activation = 'softmax')(model)
model_vgg_mnist = Model(model_vgg.input, model, 
                        name = 'vgg16')

这里首先引用VGG16模型，但是通过参数include_top=False指定迁移除顶层以外的其余网络结构到自己的模型中。Weights=’imagenet’表示借用的权重是用ImageNet数据训练出来的额。其次，通过函数方法在修改过的VGG16模型上构造一个新的扁平层用来连接新构造的全连接层，这个全连接层跟前面的模型没有区别。最后把修改过的VGG16模型和新的顶层叠加起来并赋予新的名字vgg16。这样就得到了一个基于VGG16的新模型。

## 使用KERAS构造时间序列预测模型

时间序列是在商业数据或者工程数据中经常出现的一种数据形式，是以时间为次序排列，用来描述和计量一系列过程或者行为的数据的统称。比如每天商店的收入流水或者某个工厂每小时的产品产出都是时间序列数据。一般研究的时间序列数据有两种类型。最常见的是跟踪单一的计量数据随时间变化的情况，即每个时间点上收集的数据是一个一维变量，这种是最常见的，通常的时间序列默认就是这种数据，也是本章研究的对象。另外一种时间序列数据是多个对象或者多个维度的计量数据随时间变化的情况，即每个时间点上收集的数据是一个多维变量，这种一般也被称为纵向数据（Longitudinal Data），但是不属于这里介绍的对象。

在这里我们介绍如何搭建一个LSTM深度学习模型来对在汉口测量的长江每月流量数据进行预测建模。该数据来源于DataMarket 的[时间序列数据库](https://datamarket.com/data/list/?q=provider:tsdl)，由澳大利亚莫纳什大学的统计学教授Rob Hyndman 创建，收集了数十个公开的时间序列数据集。

汉口长江月流量数据包含从 1865 年 1 月到 1978 年 12 月在汉口记录的长江每月的流量，总计 1368 个数据点。计量单位未知。

Figure 12。长江月度流量时间序列在一般的时间序列建模中，都需要检验数据的平稳性，因为传统时间序列建模都是建立在平稳数据的假设之上。这个数据具备非常强的年度周期性。使用传统的统计技术建模的时候都需要侦测周期性，并消除之，对消除周期性之后的数据运用ARIMA模型建模。

Figure 13。长江月度流量局部和移动平滑结果

我们可以通过周期图谱法（Periodogram）来得到主要的周期幅度。在Python中可以使用scipy.signal.periodogram来得到周期图谱。在这里我们不是使用原始数据，而是使用原始数据的自相关函数的周期图谱来计算主要周期，这样可以抵消噪音的影响。对读入pandas DataFrame的原始数据ts运行下面的程序我们可以得到如下的周期图谱和计算得到的主要周期长度。


import statsmodels.api as sm
from statsmodels.tsa.stattools import acf
from scipy import signal
import peakutils as peak
acf_x, acf_ci = acf(ts, alpha=0.05, nlags=36)
fs=1
f, Pxx_den = signal.periodogram(acf_x, fs)
index = peak.indexes(Pxx_den)
cycle=(1/f[index[0]]).astype(int)
fig = plt.figure()
ax0 = fig.add_subplot(111)
plt.vlines(f, 0, Pxx_den)
plt.plot(f, Pxx_den, marker='o', linestyle='none', color='red')
plt.title("Identified Cycle of %i" % (cycle))
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.show()
print( index, f, Pxx_den)

Figure 14。周期图谱很明显有一个周期为 12 个月的季节性。

虽然考虑到这个数据的本质是长江水文资料， 12 个月的周期是非常自然的预期，但是这个方法展示了对 ACF 序列运用周期图法（periodogram）找季节性周期的可靠性。在传统方法里，这里需要通过取间隔为12 的差分来消除周期性，得到一个尽可能平稳的时间序列，进而采用ARIMA模型建模。在Python里，单周期的时间序列数据，知道周期的长度以后可以直接使用季节性ARIMA模型（SARIMA）来训练。但是在使用循环神经网络模型的时候我们不用考虑这些情况，可以直接使用长短记忆模型。此外，在使用LSTM这种序列模型的时候在使用LSTM对这种单一时间序列进行建模的时候，一般通过一下步骤：

1. 将数据标准化为[0，1]区间。
2. 按照LSTM的要求，将输入数据组织为[样本数，时间步，特征变量数]的三位格式来组织。
3. 定义一个LSTM深度学习模型，通常为一个序列模型对象，逐层添加LSTM层或者其他层，最后通过一个全连接层输出到输出层。
4. 最后对需要的时间段进行预测。首先对数据进行标准化，我们使用sklearn包里的MinMaxScaler函数：


scaler = MinMaxScaler(feature_range=(0, 1))
trainstd = scaler.fit_transform(train.values.astype(float).reshape(-1, 1))
teststd = scaler.transform(test.values.astype(float).reshap

其次，我们将训练数据和测试数据组织成需要的格式，这个格式与我们将要建立的LSTM模型有关。这里我们对每个输入构造一个LSTM神经元，一个60个输入单元，每一个对应一个时间步。这60个单元的输出会作为一个全连接层的输入，这个全连接层直接产生下K个连续时间步的输出预测。作为防止过度拟合的正则化手段，我们在LSTM层和全连接层之间加了一个Dropout层。这个Dropout层在训练的时候会随机放弃一部分权重的更新，但是在进行预测的时候所有权重都会被用到。

Figure 15。LSTM网络结构（修改自CNTK Tutorial）

对于这样的网络结构，我们需要如下的一个函数来定义我们的数据，即将数据组织成为[批量数，时间步数，滞后特征数]的形式。这个可以通过如下的函数来实现：


def create_dataset(dataset, timestep=1, look_back=1, look_ahead=1):
    from statsmodels.tsa.tsatools import lagmat
    import numpy as np
    ds = dataset.reshape(-1, 1)    
    dataX = lagmat(dataset, 
                   maxlag=timestep*look_back, 
                   trim="both", original='ex')
    dataY = lagmat(dataset[(timestep*look_back):], 
                  maxlag=look_ahead, 
                  trim="backward", original='ex')
    dataX = dataX.reshape(dataX.shape[0], 
                          timestep, look_back)[:-(look_ahead-1)]
    return np.array(dataX), np.array(dataY[:-(look_ahead-1)])

执行下面的命令就可以生成所需数据：


lookback=1
lookahead=24
timestep=60
trainX, trainY = create_dataset(trainstd, 
                        timestep=timestep, 
                        look_back=lookback,  look_ahead=lookahead)
trainX, trainY = trainX.astype('float32'), trainY.astype('float32')
truthX, truthY = create_dataset(truthstd, 
                        timestep=timestep, 
                        look_back=lookback, look_ahead=lookahead)

有了如图11的网络图以后，就可以开始定义我们的LSTM深度学习模型。


batch_size=100
model = Sequential()
model.add(LSTM(48, batch_size=batch_size, \
input\_shape=(timestep, lookback), kernel_initializer='he_uniform'))
model.add(Dropout(0.15))
model.add(Dense(lookahead))
model.compile(loss='mean_squared_error', optimizer='adam')

调用fit方法就可以快速的训练这个模型。我们指定迭代20次，小批量数为100）：


model.fit(trainX, trainY, epochs=20, batch_size=batch_size, verbose=1)

下图展示了拟合过程的信息：

Figure 16。LSTM拟合过程信息

那么这个模型的拟合效果如何呢？

Figure 17。LSTM拟合结果

我们看到拟合效果还不错。平均绝对误差百分比（MAPE）只有25%不到，比用传统的SARIMA模型效果要好点。其次，LSTM模型一次输出未来24个时间点的预测值，使用起来比用SARIMA迭代预测方便很多。另外需要指出的是我们也可以直接在模型中指定损失函数为MAPE，这样更好优化衡量指标。

## 小结

在这篇短文中，我们介绍了一个目前正在流行起来的深度学习建模环境KERAS。这个建模环境相对于传统的计算环境，比如CNTK，TensorFlow，Theano等具有抽象性高，易用性好的特点，同时又依托于这几种计算环境，具有一定的可拓展性，非常适合于从事深度学习的实践者使用。我们看到使用KERAS可以非常直观地描述神经网络结构，几乎可以达到所见即所得的情况。我们在文中还分别介绍了三种流行的应用领域，分别是：

- 深度推荐模型，运用嵌入技术可以将不同类型的信息有机结合在一起构造一个深度神经网络推荐系统；
- 图像识别模型，运用多层卷积神经网络对图像进行切割分析，得到一个精度很好的手写数字分类器。同样的技术和模型可以直接移植到其他物体识别数据上，比如CIFAR10等。我们也介绍了运用已经训练好的现成模型进行迁移学习的手段，减少拟合的参数量，在保持一定精度的情况下提高训练速度；
- 简单时间序列预测模型，运用长短记忆（LSTM）神经网络模型来有效预测具备一定周期性的时间序列模型。一个非常简单的单层LSTM模型既可以达到定制的SARIMA模型的预测精度。

如果读者对如何安装KERAS这个软件以及KERAS应用于其他领域感兴趣，可以阅读由电子工业出版社出版的[**《KERAS快速上手：基于Python的深度学习》**]

CNTK integrated into Keras

2017-06-12T18:44:00.003-07:00

Keras is a very convenient tool to build your deep learning model from scratch, it is so easy to use that it almost becomes the de-facto deep learning modeling framework in Kaggle competition.

Keras used to support only TensorFlow and Theano. CNTK, now rebranding as Cognitive Toolkit in full name, from Microsoft is joining as a new backend choice, but you need to download Microsoft private fork of Keras at this moment while Microsoft is working with Keras author to fully integrate CNTK as its backend.

CNTK is a very good deep learning tool: 1. it is super fast especially for RNN type of model; 2. it scales well across multiple GPUs Here is our own speed comparsion (number in sec) on Windows platform using NVidia Titan Xp, YMMV though:

If you want to give CNTK a try, you can follow the installation instruction @ here. To install Keras that supports CNTK, follow the instructions @here.

Google's switch to G Suite

2017-02-20T15:24:00.001-08:00

My paymenet information for this custom URL is expiring and Google told me to update the payment information. That's OK, I will do it. However, by switching to G Suite, Google makes paying this $10.00 so difficult that I decided it is not worthy to continue the business with Google. Its email contains link that won't work for me, taking me to webpages that do not answer the question. All the steps become very bumpy and it is really a bad experience for customers. With this, now I understand why Google is "good at inventing technologies but bad at making products". I will switch to my WordPress sites eventually. This dying out of sas-programming is also chronical with the fading out of SAS from new emerging markets (Surely SAS will still thrive in traditional analytics market for a while).

Install Theano under Anaconda3 Python 3.5

2016-08-11T14:25:00.003-07:00

Update on 2017/01/05:

With the release of Anaconda3-4.2.0 in September 2016, users are able to install mingw and libpython under conda, which makes using Theano and keras in Python3.5 much easier.
Here are the simple steps:
1. Install Anaconda3-4.2.0. I used their Anaconda3-4.2.0-Windows-x86_64.exe installer, MD5=0ca5ef4dcfe84376aad073bbb3f8db00
2. In your Anaconda Prompt, execute : >>conda install -c anaconda mingw libpython
3. Install Theano and keras. I used their github repository
4. Remove the following statements from __init__.py file under Theano's installation folder:

if sys.platform == 'win32' and sys.version\_info[0:2] == (3, 5):

raise RuntimeError( "Theano do not support Python 3.5 on Windows. Use Python 2.7 or 3.4.")

5. Input the following statement in Environment Variable under system settings:

THEANO_FLAG=floatX=float32,device=cpu,blas.ldflags=-LC:\\openblas -lopenblas

Adding GPU capability follows the general GPU installation rule as well.

Now you are good to go, enjoying Theano backed deep learning and Python35

=======================

As of writing, Deep Learning package Theano can't be installed under Python 3.5, which causes some problems with my project. ruslanagit provided a viable solution that works for me under Windows 10.

For convenience, I copied his solution below:

Well, the main problem with Python 3.5 on windows, I guess, is that mingw and libpython are not available (not complied with Python 3.5), so you cannot run $ conda install mingw libpython step.

So you either need to downgrade to Python 3.4 (was not an option for me) and then follow standard instructions for installing Theano on Windows, or make a few tricks to make theano work with Python 3.5. For me the following steps worked on Windows 10 with Anaconda3 and Python 3.5:

Install mingw from https://sourceforge.net/projects/mingw-w64/
Add the bin directory of mingw to PATH, and make sure there is no other gcc compiler in PATH (i.e. TDM-GCC is not there).
In .theanorc file add

[gcc]
cxxflags = -shared -I"[TDM-GCC path]\include" -I"[TDM-GCC path]\x86_64-w64-mingw32\include"

You should update paths to TDM-GCC according to your system. Note, that TDM include directory is required, since the compilation will fail with mingw include directories for python 3.5 (I think they would work for Python 2, but I am not 100% sure)

Create libpython35.a manually and copy it to appropriate directory. For example:

Create temp directory
Copy python35.dll (I took it from Anaconda3 folder) into created directory
Navigate into created directory
Run: gendef python35.dll
Run: dlltool --dllname python35.dll --def python35.def --output-lib libpython35.a
Copy libpython35.a into Anaconda3\libs
All other installations/configurations were done as described in the guide for installaing Theano on Windows.

Add Table of Contents to your Jupyter Notebook

2016-08-10T16:13:00.002-07:00

Jupyter notebook does not support Table Of Contents (ToC) out of box, unlike R Markdown. There are nbextensions that ppl made to support this. It is easy to add this capability, following these steps to add Table of Contents to your notebook:

1. Find your Jupyter Notebook data directory by using this command in your DOS console:
          a. Jupyter --data-dir

2. Copy the following two files into your Jupyter data directory:
   https://rawgithub.com/minrk/ipython_extensions/master/nbextensions/toc.js
           https://rawgithub.com/minrk/ipython_extensions/master/nbextensions/toc.css

3. Run this command:
           a. Jupyter nbextension enable toc --user

4. In your Jupyter Notebook, add this command to a MarkDown cell on top of your notebook:
           a. < div id="toc" > < /div >

One caveat is that all headings will be listed in ToC, including your Title (if there is one). One easy work around is to write your title in Word, save a screenshot into an image, save the image file to the same folder as your notebook is, then use the following command in the top Markdown cell:

< div style="float: left;" > < img src="Title.png" height="200" width="800"> < /div >

Fast SQL moving average calculation without windowing functions

2015-05-10T18:09:00.001-07:00

In this post, I show a trick to do moving average calculation (can be extended to other operations requiring windowing functions) that is super fast.

Often, SAS analysts need to conduct moving average calculation and there are several options by the order of preference:

1. PROC EXPAND
2. DATA STEP
3. PROC SQL

But many sites may not licensed SAS/ETS to use PROC EXPAND and doing moving average in DATA STEP requires some coding and is error prone. PROC SQL is a natural choice for junior programmers and in many business cases the only solution, but SAS's PROC SQL lacks windowing functions that are available in many DBs to facilitate moving average calculation. One technique people usually use is CROSS JOIN, which is very expensive and not a viable solution for even medium sized data set. In this post, I show a trick to do moving average calculation (can be extended to other operations requiring windowing functions) that is super fast.

Consider the simplest moving average calculation where the trailing K observations are included in the calculation, namely MA(K), here we set K=5. We first generate a 20 obs sample data, where variable ID is to be used for windowing and the variable X is to be used in MA calculation, and then we apply the standard CROSS JOIN to first examine the resulting data, Non-Grouped, just to understand how to leverage the data structure.


%let nobs=20;
%let ndiff=-5;
data list;
     do id=1 to &nobs;
     x=id*2;
  output;
  end;
run;

options notes;
options fullstimer;
proc sql;
      create table ma as
   select a.id as aid, b.id as bid, a.id-b.id as diff, a.x as ax, b.x as bx
   from   list as a, list as b
   where a.id>=b.id and (a.id-b.id)<= abs(&ndiff)-1   
      having aid-bid>=(&ndiff+1)
   order by aid, bid
   ;
quit;

From the resulting data set, it is hard to find a clue, now let's sort by "bid" column in this data set:

From this sorted data, it is clear that we actually don't have to CROSS JOIN the whole original data set, but instead, we can generate an "operation" data set that contains the difference value, and let the original data set CROSS JOIN with this much smaller "operation" data set, and all the data we need to use for MA calculation will be there. Now let's do it: CROSS JOIN original data with "operation" data, sort by (a.id+ops), which is actually "bid' in sorted data set;



%let ndiff=5;
data operation;
     do ops = 0 to &ndiff by 1;
     weight=1;
     output;
  end;
run;
proc sql;
      create table ma2 as
   select a.id as aid, b.ops, a.id+b.ops as id2, a.x*b.weight as ax
      from   list as a, operation as b    
   order by id2, aid
      ;
quit;

Note that in above code, it is necessary to have a.x multiply by b.weight so that the data can be inter-leaved, otherwise the same X value from original table will be output and MA calculation will be failed. The explicit weight variable actually adds in more flexibility to the whole MA calculation. While setting it to be 1 for all obs result in a simple MA calculation, assign different weights will help to resolve more complex MA computing, such as giving further observations less weight for a decayed MA. If different K parameter in MA(K) calculations are required, only the operation data set need to be updated which is trivial job. Now the actual code template for MA(K) calculation will be:


%let ndiff=5;
data operation;
     do ops = 0 to &ndiff by 1;
     weight=1;
     output;
  end;
run;
proc sql noprint;
      select max(id) into :maxid
   from   list;
quit;
proc sql;
      create table ma2 as
   select a.id+b.ops as id2, avg(a.x*b.weight) as MA
      from   list as a, operation as b 
      group by id2 
   having id2>=&ndiff
     and  id2<=&maxid
   order by id2  
      ;
quit;

With this new method, it is interesting to compare it to the expensive self CROSS JOIN as well as to PROC EXPAND. On my workstation (Intel i5 3.8Ghz, 32GB memory, 1TB 72K HDD), self CROSS JOIN is prohibitively long in running time (if data is large) while the new method uses only 2X as much time as PROC EXPAND, both time consumptions are trivial comparing to self CROSS JOIN. Time consumption shown below is in "second".

Below is the code readers can run and compare yourselves.



%macro test(nobs, ndiff);
options nonotes;
data list;
     do id=1 to &nobs;
     x=id*2;
  output;
  end;
run;



%let t0 = %sysfunc(time());
options fullstimer;
proc sql;
      create table ma as
   select a.id, avg(b.x) as ma
   from   list as a, list as b
   where a.id>=b.id and (a.id-b.id)<= &ndiff-1
   group by a.id
   having id>=abs(&ndiff)
   ;
quit;
%let t1 = %sysfunc(time());


proc expand data=list out=ma2 method=none;
     convert x=ma / transformout=(movave 5);
run;


%let t2 = %sysfunc(time());

%let ndiff=5;
data operation;
     do ops = 0 to &ndiff-1 by 1;
     weight=1;
     output;
  end;
run;
proc sql noprint;
      select max(id) into :maxid
   from   list;
quit;
proc sql;
      create table ma3 as
   select a.id+b.ops as id2, avg(a.x*b.weight) as ma
   from   list as a, operation as b
   group by id2   
      having id2>=abs(&ndiff)
      ;
quit;

%let t3 = %sysfunc(time());

%let d1 = %sysfunc(round( %sysevalf(&t1 - &t0), 0.001));
%let d2 = %sysfunc(round( %sysevalf(&t2 - &t1), 0.001));
%let d3 = %sysfunc(round( %sysevalf(&t3 - &t2), 0.001));
%put SelfXJoin : &d1, EXPAND: &d2, Operation: &d3;
options notes;
%mend;



dm log 'clear';
%test(5000, -40);

%test(10000, -40);

%test(15000, -40);

%test(20000, -40);

%test(25000, -40);

%test(30000, -40);

sklearn DecisionTree plot example needs pydotplus

2015-04-25T17:33:00.002-07:00

In Python, sklearn (scikit-learn)'s DecisionTree example uses pydot for plotting the generated tree: @here.

But for Python 3, pydot has some issues with the string from dot_data.getvalue(), for example it will report "TypeError: startswith first arg must be str or a tuple of str, not bytes", @here .

While I was trying to modify related functions in pydot package, I figured that "pydotplus" is a better alternative, which you can download from @here.


>>> 
>>> from sklearn import tree
>>> from sklearn.datasets import load_iris
>>> from sklearn.externals.six import StringIO
>>> 
>>> iris = load_iris()
>>> 
>>> classifier = tree.DecisionTreeClassifier(random_state=99999)
>>> 
>>> fitresult = classifier.fit(iris.data, iris.target)
>>> 
>>> import pydotplus
>>> 
>>> dot_data=StringIO()
>>> tree.export_graphviz(fitresult, out_file=dot_data)
>>> graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
>>> 
>>> 
>>> 
>>> graph.write_pdf(r"c:\\data\\iris.pdf")
True
>>>

Migrating code pieces to GitHub

2015-02-05T11:08:00.002-08:00

One of the original reasons for this blog was to keep track of my SAS code as well as its relevant context.

That was the mindset when I was a SAS analyst, but now working in professional software company, using the right tool for versioning, collaborating and code change tracking is the professional way to go.

So I decided to set up my own GitHub repository and will migrate my code pieces to it gradually. You can find my repository @ : https://github.com/xieliang/SAS

%SVD macro with BY-Processing

2014-12-18T01:30:00.000-08:00

For the Regularized Discriminant Analysis Cross Validation, we need to compute SVD for each pair of $(\lambda, \gamma)$, and the factorization result will be feed to the downdating algorithm to obtain leave one out variance-covariance matrix $\hat{\Sigma}_{k\v}(\lambda, \gamma)$. The code is a direct modification of my original SVD macro @here. It is much more than simply put "BY-" statement in PROC PRINCOMP, but still easy to do. I used the within class CSSCP from the Iris data shipped with SAS and cross verified with R. %CallR macro is used here.

The screenshot below shows the SAS macro generates the same result as from R. Here is the test code.


dm log 'clear';
proc discrim data=sashelp.iris noclassify noprint
               outstat=iris_stat(WHERE=(_TYPE_ = "CSSCP" & Species ^= " "))
      ;
   class Species;
   var   SepalLength SepalWidth  PetalLength  PetalWidth;
run;
            
%svd_bystmt(sigma, output_v, output_s, output_u, &covars, Species, , nfac=0);

data _null_;
     file "d:\data\sigma.csv"  dlm=",";
  set iris_stat;
  if _n_=1 then do;
     put 'Species'      ',' 
            '_NAME_'       ',' 
            'SepalLength'  ',' 
            'SepalWidth'   ','   
            'PetalLength'  ','   
            'PetalWidth';
  end;
  put Species  _NAME_ SepalLength SepalWidth  PetalLength  PetalWidth;
run;

%RScript(d:\data\rscript.r)
cards4;
dat=read.csv("d:\\data\\sigma.csv", header=1, as.is=1)
head(dat)
sigma=dat[, c(-1, -2)]
sigmasvd=by(sigma, dat$Species, svd)
sigmasvd$Setosa$d
sigmasvd$Virginica$d
sigmasvd$Setosa$v
;;;;
run;

%CallR(d:/data/rscript.r, d:/data/rlog1.txt);

data _null_;
     infile "d:\data\rlog1.txt";
     input;
     put _infile_;
run;
data _null_;
     do until (eof1);
        set output_s  end=eof1;
     put @1 Species= @24 Number=  @40 Singuval=  8.2;
  end;
  put "             ";
  do until (eof2);
     set output_v  end=eof2;
  where Species="Setosa";
  put  Prin1  8.5 Prin2 8.5 Prin3 8.5 Prin4 8.5;
  end;
  stop;
run;


%macro SVD_bystmt(
           input_dsn,
           output_V,
           output_S,
           output_U,     
           input_vars,
    by_var,
           ID_var,
           nfac=0
           );

%local blank   para  EV  USCORE  n  pos  dsid nobs nstmt
       shownote  showsource  nbylevel  bystmt;

%let shownote=%sysfunc(getoption(NOTES));
%let showsource=%sysfunc(getoption(SOURCE));
options nonotes  nosource;

%let blank=%str( );
%let EV=EIGENVAL;
%let USCORE=USCORE;

%if &by_var eq &blank %then %do;
    %let bystmt = ␣
 %let nbylevel = 1;
%end;
%else %do;
    %let bystmt = by &by_var;
%end;

%let n=%sysfunc(countW(&input_vars));

%let dsid=%sysfunc(open(&input_dsn));
%let nobs=%sysfunc(attrn(&dsid, NOBS));
%let dsid=%sysfunc(close(&dsid));
%if  &nfac eq 0 %then %do;
     %let nstmt=␣ %let nfac=&n;
%end;     
%else %do;
     %let x=%sysfunc(notdigit(&nfac, 1)); 
  %if  &x eq 0 %then %do;
          %let nfac=%sysfunc(min(&nfac, &n));
          %let nstmt=%str(n=&nfac);
  %end;
  %else %do;
          %put ERROR: Only accept non-negative integer.;
          %goto exit;
  %end;
%end;

/* calculate U=XV/S */
%if &output_U ne %str() %then %do;
    %let outstmt=  out=&output_U.(keep=&by_var &ID_var  Prin:);
%end;
%else %do;
    %let outstmt=␣
%end;

/* build index for input data set */
proc datasets library=work nolist;
     modify &input_dsn;
  index create &by_var;
run;quit;

%let options=noint cov noprint  &nstmt;

proc princomp data=&input_dsn  
             /* out=&input_dsn._score */
              &outstmt
              outstat=&input_dsn._stat(where=(_type_ in ("&USCORE", "&EV")))  &options;
  &bystmt;
     var &input_vars;
run;

/* 
   need to check if the by_var is Char or Numeric, then set the
   format accordingly
*/
data _null_;
     set &input_dsn;
     type=vtype(&by_var);
  if type="C" then 
     call symput("bylevelfmtvar", "$bylevelfmt");
  else
     call symput("bylevelfmtvar", "bylevelfmt");
run;


data __bylevelfmt;
     set &input_dsn._stat  end=eof;
  &bystmt;
  retain _nbylevel 0;
  retain  fmtname "&bylevelfmtvar";
  if first.&by_var then do;
        _nbylevel+1;
  start=&by_var;  
  label=_nbylevel;
     output;
  keep label start fmtname  ;
  end;
  if eof then call symput("_nbylevel", _nbylevel);
run;
proc format cntlin=__bylevelfmt;
run;
%put &_nbylevel;


%if &output_V ne &blank %then %do;
 proc transpose data=&input_dsn._stat(where=(_TYPE_="&USCORE")) 
                out=&output_V.(rename=(_NAME_=variable))
                name=_NAME_;
   &bystmt;
      var &input_vars;
      id _NAME_;
      format &input_vars 8.6;
 run;
%end;

/* recompute Proportion */
%if &output_S ne %str() %then %do;
 data &output_S;
       set &input_dsn._stat ;
    &bystmt;
       where _TYPE_="EIGENVAL";
       array _s{*} &input_vars;
       array _x{&nfac, 3} _temporary_; 
    Total=sum(of &input_vars, 0);
    _t=0;
    do _i=1 to &nfac;
       _x[_i, 1]=_s[_i]; 
          _x[_i, 2]=_s[_i]/Total; 
       if _i=1 then 
             _x[_i, 3]=_x[_i, 2]; 
       else 
             _x[_i, 3]=_x[_i-1, 3]+_x[_i, 2];
       _t+sqrt(_x[_i, 2]);
    end;
    do _i=1 to &nfac;
       Number=_i;  
       EigenValue=_x[_i, 1]; Proportion=_x[_i, 2]; Cumulative=_x[_i, 3];
       S=sqrt(_x[_i, 2])/_t;  SinguVal=sqrt(_x[_i, 1] * &nobs/&_nbylevel);
       keep &by_var Number EigenValue  Proportion Cumulative  S SinguVal;
       output;
    end;
 run;
%end;

%if &output_U ne %str() %then %do; 
 data &output_U;
      array _S{&_nbylevel, &nfac}  _temporary_;  
      if _n_=1 then do;
         do until (eof);
            set  &output_S(keep=&by_var Number SinguVal)  end=eof;
      where Number<=&nfac;
      %if &by_var eq &blank %then %do;
          _row=1;
      %end;
      %else %do;
                   _row = input(put(&by_var, &bylevelfmtvar..), best.);
      %end;
            _S[_row, number]=SinguVal; 
            if abs(_S[_row, number]) < CONSTANT('MACEPS') then _S[_row, j]=CONSTANT('BIG');
         end;
     end;
     set &output_U;
     array _A{*}  Prin1-Prin&nfac;
     %if &by_var eq &blank %then %do;
        _row=1;
  %end;
        %else %do;
           _row = input(put(&by_var, &bylevelfmtvar..), best.);
  %end;
     do _j=1 to dim(_A);
         _A[_j]=_A[_j]/_S[_row, _j];
     end;
     keep &by_var &ID_var Prin1-Prin&nfac ;
 run;
%end;

proc datasets library=work nolist;
     modify &input_dsn;
  index delete &by_var;
run;quit;


%exit: 
options &shownote  &showsource;
%mend;

Experient downdating algorithm for Leave-One-Out CV in RDA

2014-12-14T22:42:00.001-08:00

In this post, I want to demonstrate a piece of experiment code for downdating algorithm for Leave-One-Out (LOO) Cross Validation in Regularized Discriminant Analysis [1].

In LOO CV, the program needs to calculate the inverse of $\hat{\Sigma}_{k\v}(\lambda, \gamma)$ for each observation v to obtain its new scoring function at a given pair of regularizing pair parameter $ (\lambda, \gamma) $, where $k\v$ indicates class $k$ excluding observation $v$.

As mentioned in [1], to obtain the inverse of this variance covariance matrix for class k, it is not sufficient to just remove observation $v$ and recalculate it, but instead, we try to calculate $$W_{k\v}(\lambda, \gamma)\hat{\Sigma}^{-1}_{k\\v}(\lambda, \gamma)$$, which is equivalent to compute the inverse of :

\[
W_k(\lambda)\hat{\Sigma}^{-1}_k (\lambda) - \gamma \|Z_v\|^2 \cdot I - (1-\gamma)Z_v Z_v^{t} \]

In this formula, the first element $W_k(\lambda)\hat{\Sigma}^{-1}_k (\lambda) $ does not involve quantities varying along with observation $v$ and can be pre-computed. The latter two can be easily updated on the fly as we scan through each individual observations.

With this in mind, we use the iris data for demonstration. Note that the SAS code here is for illustrating the concept of downdating algorithm above for a given pair of [ (\lambda, \gamma) ]. The computing is divided into several steps. (1) We first use PROC DISCRIM to calculate CSSCP matrix, total weight and variable mean by class. These quantities correspond to $ \Sigma_k $, $ \bar{X_k}$ and $ W_k$ in the original paper. (2) We then use the %SVD Macro @here to obtain key elements for later updating.

Note that in order to obtain CSSCP, MEAN and total weights, in addition to PROC DISCRIM, we can also use PROC CORR. The advantage is that we can keep all computing within the framework of SAS/Base so to benefit budget-sensitive business and maximize compatibility. The down side is that the data set need to sorted or indexed to obtain CSSCP by class and the data has to be passed twice, one for pooled CSSCP and the other for class specific CSSCP.

(to be completed...)

 

/* study downdating algorithm */
%let lambda = 0.3;
%let gamma = 0.6;
%let target = Species;
%let covars = SepalLength SepalWidth  PetalLength  PetalWidth;

/* step 0. Obtain key statistics from DISCRIM */
proc discrim data=iris  outstat=stat noprint  noclassify;
      class ⌖
   var   &covars;
run;

/* step 1. Obtain \Sigma_k(\lambda) & \Sigma_k(\lambda, \gamma)
           \S_k = (_TYPE_="CSSCP' & &target = "something")
           \S  = (_TYPE_="CSSCP" & &target = " ")
           \W_k = (_TYPE_="N" & &target = "something") -1
           \W  = (_TYPE_="N" & &target = " ") - &nClass + 1
           \Sigma_k(\lambda) = \S_k(\lambda) / \W_k(\lambda)
              where:
                    \S_k(\lambda) = (1-\lambda)*S_k + \lambda*S
                    \W_k(\lambda) = (1-\lambda)*W_k + \lambda*W

          \Sigma_k(\lambda, \gamma) = (1-\lambda)*\Sigma_k(\lambda) + \gamma* trace[\Sigma_k(\lambda)]/p*I         
*/ 
proc sql noprint;
     select distinct &target into :targetvalues separated by " "
  from   stat
  where _TYPE_="CSSCP"
  ;
  select distinct _NAME_ into :covars separated by " "
  from   stat
  where _TYPE_="CSSCP"
  ;
quit;
%put &targetvalues;
%put &covars;

%let nClass=%sysfunc(countw("&targetvalues"));
%put &nClass;
%let nVars=%sysfunc(countw("&covars"));
%put &nVars;

data weights;
     set stat;
  where _TYPE_="N";
  array _nv{*} _numeric_;
  if &target="" then _total=_nv[1]-&nClass+1;
  retain lambda λ
  retain _total;
  weight=(1-&lambda) * (_nv[1]-1) + &lambda*_total;  
  keep &target   weight  lambda;
  if &target ^= "";
run;

proc sort data=stat  out=sigma;
      by _TYPE_  ⌖
   WHERE _TYPE_ = "CSSCP" & &target ^= " ";
run;
proc sort data=weights;
     by  ⌖
run;

/*-  step 1.1 Update \Sigma_k(\lambda) = ((1-\lambda)*S_k + \lambda*S) / ( (1-\lambda)*W_k + \lambda*W )  -*/
                                                                                          
data sigma;
     array _sigma{&nVars, &nVars} _temporary_;
  array _x{&nVars} &covars;
  array _y{&nClass} _temporary_;  
  if _n_=1 then do;
     _iter=1;
     do until (eof1);
        set stat(where=(_TYPE_="CSSCP" & &target = " "))  end=eof1;
     do _j=1 to &nVars;
        _sigma[_iter, _j]=_x[_j];
     end;
     _iter+1;
     end; 
        
     _iter=1; 
     do until (eof2);
     set weights end=eof2;
        _y[_iter]=weight;
     _iter+1;
  end;
  put _y[3]=;
  
  _target=⌖
  _iter=0;
  end;
  modify sigma ;  
  retain  weight;
  retain  _target;
  _TYPE_="COV";
 
  if &target^=_target then do;
        _iter+1;
        _i=0; 
  weight=_y[_iter];
  _target=⌖
  end;
  _i+1;
  put "&covars";
  put &covars;
  do _j=1 to &nVars;
     _x[_j]= ((1-&lambda)*_x[_j] + &lambda*_sigma[_i, _j])  ;
  end;
  put &covars;
  put "-----------------------------------------------------------------";
  replace;  
run;

/*- trace is sum of every elemnt of a matrix -*/
data _trace;
     set sigma;
  by ⌖
  array _x{*} &covars;
  retain trace _iter;
  if first.&target then do;
     _iter = 1;
  trace=0;
  end;
  trace+_x[_iter];
  _iter+1;
  if last.&target then output;
  keep &target  trace;
run;


/*- step 1.2 Update \Sigma_k (\lambda, \gamma) = (1-\gamma)*\Sigma_k(\lambda) + \gamma *c*I -*/
/*- where  c = trace(\Sigma_k(\lambda)) / p, where p=&nVars                             -*/
data sigma;
      if _n_=1 then do;
      declare hash h(dataset:'_trace');
   h.defineKey("&target");
   h.defineData("trace");
   h.defineDone(); 
   call missing(trace);
   end;
   set sigma; by ⌖
   array _x{*} &covars;
   retain _adj  trace _iter;
   if first.&target then do;
      _iter=1;
   end;
   if _iter=1 then do;
         _rc=h.find();
         _adj= &gamma/&nVars * trace;
   end;
   put _iter=  &target=;
   do _j=1 to &nVars;
      if _j=_iter then do;
      _x[_iter] = (1-&gamma)*_x[_iter] + _adj;
   end;
   else do;
            _x[_iter] = (1-&gamma)*_x[_iter];
   end;
   end;
   _iter=_iter+1;
   drop _type_   _adj _iter _j  _rc  trace;
run;



/*-step 1.3 Obtain the inverse of W_k \Sigma_k(\lambda, \gamma) - c*I  using SVD                                     -*/  
/*-This inverse is the same as first inverse \Sigma_k(\lambda, \gamma), then with S matrix deflated by W_k(\lambda)  -*/
/*-and differenced by c, where c=\gamma* |Z_v|^2/p, and Z_v = sqrt(bk(v))*(X_v- mean(X_k\v),                         -*/
/*-bk(v)=sk(v)*Wc(v) wv/(Wc(v)-wv)                                                                                   -*/
/*-v is the current observation                                                                                      -*/

/*-       1. apply SVD to obtain eigenvalues   
%SVD(sigma, out_u, out_S, out_V)                                                                                     -*/


/*-       2. The whole CV will be done in one Data Step
             Because all related quantitaties of the rest, such as such as |Z_v|, sk(v), Wc(v) 
             will require at least one pass through of the original data, therefore it is
             best to incorporate all computing in this step.
-*/
data _iris;
     array _m{*} &covars;
     array _mt{&nVars} _temporary_;
     if _n_=1 then do;
     _class=1;
        do until eof;
     set stat(where=(_TYPE_='MEAN'& &target^=""))  end=eof;
     do _j=1 to dim(_m);
        _mt[_class, _j]=_m[_j];
     end;
     _class+1;
  end;
  end;
  array _rank1{&nVars, &nVars} _temporary_;
  set iris point=1;  
  do _j=1 to dim(_m);
     _m[_j]=_m[_j]-_mt[1, _j];
  end;
  do _j=1 to &nVars;
     do _k=1 to &nVars;
           _rank[_j, _k]=_m[_j] * _[_k];

           /*  MORE CODE TO COME*/

     end;
  end;
run;

[1] J. H. Friedman. Regularized discriminant analysis. Journal of the American Statistical Association, 84(405):165–175, 1989.

Control Excel via SAS DDE & Python win32com

2014-12-14T22:41:00.006-08:00

Excel is probably the most used interface between human and data. Whenever you are dealing with business people, Excel is the de facto means for all things about data processing.

I used to only use SAS and Python for number crunching but in one of my recent projects, I need to work out a way to send summarized data to Excel for charting, specifically charting in the Excel flavor where the style has particularity and clients need data to be associated with the chart, just like a copied Excel chart. It is extremely difficult to fine tune your SAS or Python figures to mimic the exact required Excel style, no to mention interactive data display on chart that comes handy with Excel, so I decided to study how to control Excel via SAS and Python.

Python is for standalone PC implementation while SAS runs on the analytic server, each serves different customer groups. In my situation, we have MS Office 2013 Professional installed on both client PC and the server. In Python, we use win32com to control Excel while in SAS, I use the old school Dynamic Data Exchange (DDE) . DDE is obsolete technology but I don't have much experience in using COM within SAS, nor did I found good examples in this aspect.

In this blog, I show how to generate the following Excel worksheets and copy the Charts with data to PowerPoint, in Python and in SAS.

For python scripting, I largely follow the posts from Dan's Python Excel blogs, and the Office IPython blog @here. You are welcome to check out his scripts there and develop your own. In order for win32com to work in your Python environment, it is required to install pywin32 package, which can be found on its project page @here . You may also want to create static dispatch use makepy.py from the project code and installs all generated constants from a type library in the object win32com.client.constants.

For SAS scripting, I largely follow SAS manual outlined above and some SAS user papers on DDE, especially using DDECMD. But since in PowerPoint 2013 macro recording is not longer available, the work will stop at Excel.

In the following, I will demonstrate, via comparing both languages, how to:
1. Invoke Excel Application
2. Reading data and generate contingency tables
3. Send the table data to Excel
4. Generate charts in Excel driving by SAS or Python
5. Copy Excel charts to PowerPoint with data embedded.

(to be completed .. )

%HPGLIMMIX SAS macro is available online at JSS website

2014-07-01T10:52:00.000-07:00

My paper "%HPGLIMMIX: A High-Performance SAS Macro for GLMM Estimation" is now available at Journal of Statistical Software website @here.

SAS macro and code can also be found there. If you use it, please kindly send me an email so that I know my work is appreciated and is helpful to more people.
===========

Update:

The paper and the macro are well accepted: 400+ downloads for the paper and 90+ downloads for the macro for the first month it is online and I've got many good feedbacks. Several key points to make:

1. For binary dependent variable, you have to code your data into numeric {0, 1} and the macro will model the probability dependent variable having value 1. So if you have other requirements, you have to process the data before calling the macro.

2. This macro relies on PROC HPMIXED, so that it will generate two residual parameters by default, therefore for conditional distributions without scale parameters, such as Poisson, users need to hold the last residual variance to be 1 by using PARMS /HOLD= statement. Please see the explanation of second difference between this macro and PROC GLIMMIX on Page 9, as well as Example 2.

Market trend in advanced analytics for SAS, R and Python

2013-12-06T00:16:00.004-08:00

Disclaimer: This study is a view on the market trend on demand of advanced analytics software and their adoptions from the job market perspective, and should not be read as a conclusive statement on what is all happening there. The findings should be combined with other third party research results such IDC reports http://www.idc.com/getdoc.jsp?containerId=241689 to reach a balanced and comprehensive idea.

/* --------------------------------------------------------------------------------------------------------------------*/
The debate on the competition between SAS, R and Python for advanced analytics, Statistics or Data Science in general never ends. The question boils down to how far SAS, as a software and programming language for such purpose, can still be used by business in mass, given its current pricing structure.

Some recent discussions can be found at:

http://www.listserv.uga.edu/cgi-bin/wa?A2=ind1312a&L=sas-l&F=&S=&P=1171

http://r4stats.com/2013/03/19/r-2012-growth-exceeds-sas-all-time-total/

http://r4stats.com/articles/popularity/

and more...

indeed.com has been used to address such question before by Robert Muenchen @ http://r4stats.com/articles/popularity/. His analysis on popularity of SAS and R in the job market using indeed.com was simple. The terms he used were just "R SAS" and "SAS !SATA !storage !firmware".

In my analysis using indeed trend, I directly search combination of languages such as R, SAS and Python with analytics related job terms, such as "R AND Analytics" or "R AND Regression", etc. The goal is to understand the market dynamics in adopting each of the three languages in advanced business analytics. I try to be as fair as possible in my analysis.

The job description related terms that I am going to combine with the language names fall into three categories:

1. Techniques typically used, here I use "Regression", which almost all advanced business analytics job will look for;
2. Industry, here I use "pharmaceutical", and I used "JPMorgan Chase" to represent financial service (even though it is not representative"
3. General term such as "data mining", "machine learning" and "statistics". In general, data mining has different implication from machine learning. The former is a more general term, or a modern term for "analytics" while the latter shows up momentum only recently and largely focus on computer science related field and more hardcore on algorithm development, etc, so favoring Python over the other two.

The graph below shows indeed.com job trends using search term "R and regression, SAS and regression, Python and regression". It has three immediate pieces of information to tell:
a.) Python has been picking up very fast since 09;
b.) The market share gap between SAS and R/Python is consistently dropping since 2006 and in foreseeable future the lead of SAS is going to disappear.
c.) SAS's market share reached high peak in 2011 but kept dropping thereafter while Python maintains steady trend and R is picking up in the same time period.

The below graph searches the trend using term "R and pharmaceutical", etc. I wanted to see how the trend of adopting those three languages in pharmaceutical industry. What we found out are three points:
a.) R and SAS almost show up at the same time in pharmaceutical jobs
b.) The relative job demand for analytics in pharmaceutical industry is declining.
c.) Python is on the rise even though still trailing far behind SAS and R.

When using "JPmorgan Chase" combined with languages for trend analysis, I specifically excluded "storage" from the SAS search term because the bank may hire some IT folks and result will mix with SAS storage system. The term is ' R and "jpmorgan chase", SAS not Storage and "jpmorgan chase", Python and "jpmorgan chase" '. There are several interesting observations:
a.) Python is picking up some job shares in JPmorgan Chase from 2010. R is losing;
b.) There is noticeable seasonality in hiring, mostly in summer time;
c.) Before 2008, there was strong demand for analysts using SAS, then the crisis came and hiring stalled for almost 2 years. Beginning in 2011, the business is getting better and the demand for analysts using SAS is picking up all year around and attained peak in this year.

In our last part of trend analysis, we want to combine the languages with more general terms such as "analytics", "data mining" to see how the whole business perceives these three analytical languages.

We first combine with "data mining". Three points need to be noted:
1. The job market for analytics in general has been increasing since 2009 and all languages shows good momentum from 2009 to 2010;
2. SAS was the dominant player in this business, but the landscape is changing rapidly. Both R and Python see almost identical adoption by the business. There are some minimal gap between these two but I would consider it not significant;
3. SAS is not as favored by the market as before and its share not only stalled but also showed some declining trend this year. This is very alarming.

Now, let's be more focused. First, we study "Statistics". The message is not good for SAS here. SAS used to be THE de facto software for industrial statistical analysis. This fact is reflected by the trend before 2010. SAS is so dominant that its market share is more than 3X comparing to R and even more comparing to Python. But what's interesting is that while SAS saw some good days from 2010 to 2011, which also reflected in their annual revenue", its market share is on a steady down road from 2012. On the other hand, R and Python are still picking up without hesitation.

Turning to "machine learning", Python is the leader, followed by R. SAS sees some distance here. The seasonality of Python is a little different from R and SAS while the latter two have very close seasonal dynamics. The hypothesis here is that Python is dominantly used in different industries while R and SAS share similar industries.

I don't always do regression, but when I do, I do it in SAS 9.4

2013-07-18T21:39:00.000-07:00

There are several exciting add-ins from SAS Analytics products running on v9.4, especially the SAS/STAT high performance procedures, where "high performance" refers to either in single-machine multi-threading mode or full distributed mode.

HPGENSELECT: high performance model fitting and variable selection for GLM within standard exponential family;

HPLOGISTIC: high performance logistic regression

HPLMIXED: high performance Linear Mixed Model, unfortunately it still strictly a LMM, so does not support the DIST= option in MODEL statement

HPNLMOD: high performance nonlinear LS or MLE

HPREG: high performance for OLS, luckily HPREG supports CLASS statement and includes LAR/LASSO variable selection methods. On the other hand, HPREG does not supports best-subsets variable selection methods.

HPSPLIT: high-performance utility procedure that creates a decision tree model and saves results in output data sets and files for use in SAS Enterprise Miner.

Finding the closest pair in datat using PROC MODECLUS

2013-05-08T20:58:00.003-07:00

UPDATE: Rick Wicklin kindly shared his visualization efforts on the output to put a more straightforward sense on the results. Thanks. Here is the code, run after my code below. Note that this is designed for K=2.

proc iml;

use out;      read all var {ID x y}; close out;
use neighbor; read all var {Nbor}; close neighbor;
NBor = num(Nbor);         /* convert to numeric */
xEnd = x[Nbor];
yEnd = y[Nbor];
create Viz var {ID x y xEnd yEnd}; append; close Viz;
quit;
proc sgplot data=Viz noautolegend;
        scatter y=y x=x / datalabel=id;
        vector x=xEnd y=yEnd / xorigin=x yOrigin=y;
run;

***************************************************************************

More often than not, the analyst wants to find the pair of observations that are the closest in terms of certain metrics, such as Euclidean distance. For example, on the popular SAS-L archive group, Randall Powers posted such a question at here, but he was using PROC DISTANCE.

His idea is straightforward: First calculate the pairwise distance among all observations, and Second, find the one that has the shortest distance. For the first step, he was trying to use PROC DISTANCE, then then he planned to use DATA STEP to search through. This is not workable in even the modest data size, say 100K observations, because PROC DISTANCE will generate a dense double precision 100K-by-100K matrix and will blow up your workstation. That is the reason Randall turns to SAS-L for help.

Well, SAS actually has an old PROC that can readily solve his problem. Of course, before a solution can be crafted, there are several issues that have to be clarified first.

First, what if some observations share the same closest observations, do you allow shared case or you need a mutually exclusive solution that each observation only has its unique non-shared counterpart as the closest point, then how would you determine which one should take the shared observation? First come first serve or any other rule?

Second, when we talk about pairs, we know N should be even number, but what if the data we have has an odd number for N? How to deal with the left over?

That being said, let's consider the simplest and most tolerate case, that is each observation is allowed to serve as the closest point to multiple other observations and N is an even number.

The solution uses PROC MODECLUS. This is an old procedure for density based clustering analysis, but its build-in nearest neighbor search algorithm and output capability make it the perfect candidate for this type of job. We use the example data in the "Getting Started" section.

The ODS OUTPUT statement directly outputs the list of nearest neighbors (closest points) for each observation and your have to specify either ALL or NEIGHBOR option in the PROC MODECLUS statement in order to use this functionality. In the same statement, we also specify K=2. K=2 means 2 nearest neighbor but since this applies for all observations, for each one, the nearest neighbor is K-1=1. So if you specify K=3, you actually ask this procedure find 2 closest point for each 1 observation at hand.

The good thing about the nearest neighbor list data set is that it also contains the calculated distance, therefore, in cases you need to deal with more complex situation as listed above, such as non-shared nearest neighbor, etc, you need some room for post processing. Here you can specify K=3 and PROC MODECLUS will output 2 neighbors for each observation, like below.

You can see that observation 3 can be matched to either observation 1 or 2, but with 2, it yields the shortest distance, and if this is the rule you are going to apply, your post process can work on this data set to implement the rule.

Please feel free to verify the result. If you find error, let me know or post in the comment.

 
data example;
      input x y @@;
   ID=_n_;
      datalines;
   
   18 18  20 22  21 20  12 23  17 12  23 25  25 20  16 27
   20 13  28 22  80 20  75 19  77 23  81 26  55 21  64 24
   72 26  70 35  75 30  78 42  18 52  27 57  41 61  48 64
   59 72  69 72  80 80  31 53  51 69  72 81
   ;
run;

ods select none;
ods output neighbor=neighbor; /* ODS output dataset */
proc modeclus data=example method=0  k=2    /*Find a pair of nearest neighbor*/
                 out=out  
                 all  /* required option in order to output Nearest Neighbors*/;
     var x y;
  id  id;
run;
ods select all;

Large Scale Linear Mixed Model

2013-03-25T21:44:00.000-07:00

Update at the end:
****************************;
Bob at r4stats.com claimed that a linear mixed model with over 5 million observations and 2 million levels of random effects was fit using lme4 package in R:

I am always interested in large scale mixed model like this and would appreciate anyone who can point me to the example Bob quoted (update: Bob responded in comment, please check it out). ~~I left a reply to Bob on his blog but no response yet.~~ I also tried to Google this example but all examples returned are on a scale with hundreds up to 5 thousand of levels of random effects, which is very common and reasonable. While I don't have the data used in the quoted example, I decided to see how long it takes for a mixed model with similar scale to be fitted in SAS's PROC HPMIXED. Therefore I borrowed Tianlin Wang's simulated data example in his paper: All the Cows in Canada: Massive Mixed Modeling with the HPMIXED procedure in SAS 9.2. But I modified the code a little bit so that there will be 400K cows to be included in the random effect and I used SUBJECT= option to accelerate the computing with by-subject processing.

%let NFarm = 1000;
%let NAnimal = %eval(&NFarm*400);
%let Va = 4.0;
%let Ve = 8.0;
data Sim;
keep Species Farm Animal Yield;
array BV{&NAnimal};
array AnimalSpecies{&NAnimal};
array AnimalFarm{&NAnimal};
do i = 1 to &NAnimal;
BV {i} = sqrt(&Va)*rannor(12345);
AnimalSpecies{i} = 1 + int( 5 *ranuni(12345));
AnimalFarm {i} = 1 + int(&NFarm*ranuni(12345));
end;
do i = 1 to 40*&NAnimal;
Animal = 1 + int(&NAnimal*ranuni(12345));
Species = AnimalSpecies{Animal};
Farm = AnimalFarm {Animal};
Yield = 1 + Species
+ Farm
+ BV{Animal}
+ sqrt(&Ve)*rannor(12345);
output;
end;
run;

proc sort data=sim; by animal species farm; run;
options fullstimer;
ods listing close;
ods select none;
proc hpmixed data=Sim;
class Species Farm Animal;
model Yield = Species Farm*Species /ddfm=residual;
random animal /cl ;
ods output SolutionR=EBV;
run;
ods listing;

For this particular model, it took more than 8 hours to be fitted on a desktop with i5-3570K CPU running at 3.8Ghz. I certainly dare not to test against a case with 2 million levels of random effects and would assume it will be a LONG LONG TIME. Note that with the original 10K levels of random effects, it took only 20 seconds to be fitted on the same machine, and I do know the time consumption is more than super-linear in the number of random effects' levels with all the computing optimization in HPMIXED.

lme4 package supports Sparse matrix techniques only on the random effect design matrix, but HPMIXED procedure supports Sparse matrix techniques to the full Mixed Model Equation, and in many applications, will perform much better:

Take the above example, but using only 150 Farms and assuming 100 Animals within a Farm, we fit the same mixed model using lme4 package and the HPMIXED procedure. The lmer function is called within SAS using the %RScript and %CallR macros (see here).

Left panel shows the status of the machine running lmer function on the mixed model and the right panel shows the status when the machine ran HPMIXED. Because lme4 does not support sparse matrix techniques for fixed effect design matrix, it takes significant penalty on memory usage (12+GB vs 25MB) and computing time (244sec real time vs. 17sec real time). In a later test using only random effects, HPMIXED used 5.7 seconds while lmer used 11.3 second.

Left: Running R; Right: Running SAS

Poor man's HPQLIM?

2013-02-26T19:47:00.001-08:00

Tobit model is a type of censored regression and is one of the most important regression models you will encounter in business. Amemiya 1984 classified Tobit models into 5 categories and interested reader can refer to SAS online doc for details. In SAS, PROC QLIM is used to estimate various types of Tobit models as well as related class of models, such as Heckman's selection model, etc.

Jason (Blog @ Analytics-in-writing) in his recent post showed that PROC QLIM is not designed for large data, meaning for data set with millions or more records, and SAS developed HPQLIM, the High Performance version of QLIM, to deal with large scale modeling within Tobit model framework. While I celebrate the roll out of HPQLIM, I think the example Jason used is poorly chosen. He used standard Tobit model (type 1), while he really should use a more complicated one to demonstrate the full power of HPQLIM because there already is a solution in SAS to fit standard Tobit model on very large data set very fast.

Before we dive into the example Jason mentioned, it is useful to review the concept of censoring and truncation in regression analysis. Suppose there is a variable Y following a continuous distribution D and a constant K. After a value y is sampled from Y:

if y < K (y > K) and discard the whole record, then it is left (right) truncation;

if y < K (y > K) and report y=K, then it is left (right) censoring.

Jason showed that, when facing millions of records, PROC QLIM is very slow in estimating even the Type 1 Tobit model, the simplest left censoring at 0, and suggested to use HPQLIM which shows significant improvement in computing time. Unfortunately, not every business can afford High Performance (HPx) solution from SAS but Tobit model is so prevailing in observational study, people may wonder is there a way to get around this within the SAS environment but with much improved speed?

Yes, at least for Type 1 Tobit model. SAS Manual provide an exampling using PROC LIFEREG for Type 1 Tobit model since it is left censoring and the log likelihood function is no difference from a parametric survival model with corresponding distribution:

$$ l = \sum_{i \in (y=0)} \left [ 1- \Phi(X_i\beta/\sigma)\right ] + \sum_{i \in (y>0)} \left [ \Phi(X_i\beta/\sigma)/\sigma \right ] $$

PROC LIFEREG is significantly faster comparing to PROC QLIM as it is optimized for this particular type of model. In order to demonstrate the speed gain, I used a similar data as Jason used, the training sample from KDD CUP 1998. This is the original data where the sample Jason used is from, with 95K observations and 400+ variables. In consistent with Jason's test, we will use 10 randomly selected numerical variables in the model fitting, and stack the data 100 times so that the final data feed to PROC QLIM & PROC LIFEREG has about 9.5 million observations.

The SAS log below shows that, on a laptop with i5 2520M CPU clocked at 2.5Ghz and with 4G memory, PROC QLIM used almost 39 minutes while PROC LIFEREG only used less than 37 sec, more than 60 times difference. In comparison, HPQLIM in Jason's example uses 1/3 of the time of QLIM. Based on the real time data, Jason's OC (old client) machine runs about the same speed as my laptop using QLIM, so I expect the same results from HPQLIM if I could have run one, namely, 12 minutes real time, 27 minutes CPU time. Still not comparable to LIFEREG's less than 1 minute time usage.

6092

6093 %let keepvars=target_: age wealth1 hit CARDPROM numprom cardpm12 numprm12

6094 RAMNTALL NGIFTALL CARDGIFT MINRAMNT LASTGIFT

6095 ;

6096 data pva;

6097 set train(keep=&keepvars);

6098 do i=1 to 100;

6099 output;

6100 end;

6101 run;

NOTE: There were 95412 observations read from the data set WORK.TRAIN.

NOTE: The data set WORK.PVA has 9541200 observations and 15 variables.

NOTE: DATA statement used (Total process time):

real time 1:11.55

cpu time 5.26 seconds
6102

6103 proc qlim data=pva;

6104 model target_d = age wealth1 hit CARDPROM numprom cardpm12 numprm12

6105 RAMNTALL NGIFTALL CARDGIFT MINRAMNT LASTGIFT ;

6106 endogenous target_d~censored(lb=0);

6107 run;
NOTE: Convergence criterion (FCONV=2.220446E-16) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: PROCEDURE QLIM used (Total process time):

real time 39:51.24

cpu time 39:33.01

6108

6109 /* make it suitable for parametric survival model */

6110 data pva;

6111 set pva;

6112 if target_d=. then y=0; else y=target_d;

6113 run;

NOTE: There were 9541200 observations read from the data set WORK.PVA.

NOTE: The data set WORK.PVA has 9541200 observations and 16 variables.

NOTE: DATA statement used (Total process time):

real time 1:57.75

cpu time 11.82 seconds

6114

6115 proc lifereg data=pva;

6116 model (target_d, y) = age wealth1 hit CARDPROM numprom cardpm12 numprm12

6117 RAMNTALL NGIFTALL CARDGIFT MINRAMNT LASTGIFT ;

6118 run;

NOTE: Algorithm converged.

NOTE: PROCEDURE LIFEREG used (Total process time):

real time 36.53 seconds

cpu time 13.06 seconds

One additional advantage of PROC LIFEREG is that it also can very easily handle observation-by-observation censoring, which is interval censoring in a survival analysis.

The disadvantage of PROC LIFEREG comparing to PROC QLIM in the Tobit model framework is that PROC LIFEREG only computes predicted value for uncensored observations, but in Tobit model, the expected value of censored observations is of particular interests. To do this, analyst will have to use DATA STEP to manually compute this value for each of censored observations. But this is easy.

For example, standard Tobit model usually assumes normal distribution for latent variable, and in this case, the expected observed value can be calculated as: $$ E(Y_i) = \Phi(X_i\beta/\sigma)(X_i\beta + \sigma \lambda)$$

where $$ \sigma=SCALE $$ is from estimation and $$\lambda=\phi(X_i\beta/\sigma)/\Phi(X_i\beta/\sigma)$$, where $$\phi(), \Phi()$$ indicate density function and cumulative probability function from standard normal distribution, respectively.

Alternatively, if a logistic regression is assumed, then $$ E(Y_i) = \sigma \log(1 + exp(X_i\beta/\sigma))$$

Even if with this extra step to calculate necessary quantities, the total time it would take is still significantly smaller than using PROC QLIM. As seen from the log, go through the whole data set computing these quantities with 9.5 million records took less than 2 minutes, so the total time using PROC LIFEREG is less than 3 minutes, that is still only 8% of what it takes by PROC QLIM. Recall that in Jason's poset, HPQLIM runs on a Greenplum Application with 35 working nods still takes 1 minutes 24 seconds to finish the fitting and scoring process on a data set with 7.7 million records.

Lastly but not the least, for this particular example, a Heckman's selection model is more appropriate because the intention to donate (the latent variable) is correlated with the amount of donation if the subject decides to donate, while a standard Tobit model literally assumes independence. And applying a Full Information Maximum Likelihood (FIML) estimation of Heckman's selection model on such a large data set is where HPQLIM really shine and there is no immediate replacement available in current SAS solution.

Reference:
SAS Institute, Inc., SAS/ETS User's Guide for v9.2, Cary, NC
SAS Institute, Inc., SAS/STAT User's Guide for v9.2, Cary, NC.

Amemiya, T. (1984), "Tobit Models: A Survey," Journal of Econometrics, 24, 3–61.

Kaggle Digit Recoginizer: SAS k-Nearest Neighbor solution

2012-12-10T18:58:00.003-08:00

Kaggle is hosting an educational data mining competition: Kaggle Digit Recognizer, using MNIST data. Handwritten digit recognition is one of the few applications that kNN classifier performs well. Of course, the benchmark kNN classifier provided by the competition organizer is not even close to being sophisticated. So I decided to try the kNN Classifier in SAS using PROC DISCRIM. For details about how to do kNN classifier in SAS, see here and here .

The data is pre-processed from raw images using NIST standardization program, but it noteworthy some extra efforts to conduct more exploratory data analysis (EDA). Because this is a pre-processed data, therefore no missing data issues and EDA will focus on descriptive statistics of input variables. The code below read in the data which is stored in folder "D:\DATA\MNIST";


%let folder=D:\Data\MNIST;

libname digits "&folder";
data digits.MNISTtrain;
infile "&folder.\train.csv" lrecl=32567 dlm=',' firstobs=2;
input label x1-x784;
run;
data digits.MNISTtest;
infile "&folder.\test.csv" lrecl=32567 dlm=',' firstobs=2;
input x1-x784;
run;

Figure 1 below shows the distribution of variances of input variables in the Training set. Apparently there are many input variables of little information and they should be removed before further data processing and actual modeling. Next we examine Coefficient of Variation (CoV) of input variables. CoV gives a more consistent way to compare across variables. The result is demonstrated in figure 2. Note that variables with 0 variance will have their CoV be set to 0.


ods select none;

proc means data=digits.MNISTtrain ;
var x1-x784;
output out=stat(where=(_STAT_ in ('STD', 'MEAN')));
run;
ods select all;
proc transpose data=stat out=stat_t;
var _STAT_ x:;
run;
data stat_t;
set stat_t;
ID=_n_-1;
if col2=0 then CoV=0; else CoV=col1/col2;
run;
proc sgplot data=stat_t;
scatter y=col2 x=ID;
label ID="X #" col2="StdErr";
run;
proc sgplot data=stat_t;
scatter y=CoV x=ID;
label ID="X #";
run;

Figure 1. StdErr by input features

Figure 2. Coefficient of Variation by input features

In this demonstration, we select features that have a variance larger than 10% of quantiles. Even so, we still end up with more than 700 features. Note that in SAS, the computing time of kNN classification increases in superlinear and kNN itself suffers curse of dimensionality, therefore it is desirable to conduct dimension reduction first before applying kNN classifiers. While there are so many different techniques to choose from, we decided to go with SVD, which is both widely used and simple to implement in SAS, see here for sample code and examples. The code below applies SVD to the data that combining both training and testing. This is important because SVD is sensitive to input data and include both data sets in the computing won't violate any best practice rules in predictive model because no information from dependent variable is used. The left eigenvectors are also standardized to have 0 mean and unit standard deviation.

data all;
set digits.MNISTtrain(in=_1)
digits.MNISTtest(in=_2);
if _1 then group='T';
else group='S';
run;
proc princomp data=all out=allpca(drop=x:) outstat=stat noint cov noprint;
var x1-x784;
run;

proc standard data=allpca out=allpcastd mean=0 std=1;
var Prin:;
run;

data allpcastd;
set allpcastd;
array _p{*} p1-p10 (10*0.1);
cv=rantbl(97897, of _p[*]);
drop p1-p10;
run;

After the data is compressed, we split them into 10 folds of subsets, flag each of them and search for best tuning parameter. In building a simple kNN classifier, there are two tuning parameters to decide, 1.what features to use and 2. the number of neighbors to use. In this exercise, we simplify the first problem into searching for the number of features, always starting from feature #1. We then calculate the error rate for each combination. Before we start a brutal force search, we have to understand that search through the whole space of combination of 700 features and 1 to 25 neighbors is not practical, it just takes too much time, may run for over 2 days on a regular PC. But we can narrow down to a range of number of features by testing number of features at every 10th interval, say test 1 feature, 11 features, 21 features,..., 201 features. In classification, kNN is more sensitive to features used instead of number of neighbors. The code below executes such experiment and will take about 18 hours to finish on my PC equipped with i5-3570 CPU, clocked at 3.8Ghz. In actual execution, I used my proprietary %foreach MACRO that distribute the job to 3 cores of the CPU, and shaved the total running time to only 6 hours.

/*@@@@@@@@@@ USE 10-FOLD CV @@@@@@@@@@*/

%macro loopy;
%do k=2 %to 8;
%do i=1 %to 10;
%let t0=%sysfunc(time());
options nosource nonotes;
proc discrim data=allpcastd(where=(group='T' & cv^=&i))
method=npar k=&k noprint
testdata=allpcastd(where=(group='T' & cv=&i))
testout=out(drop=Prin:)
;
class label;
var Prin1-Prin30;
run;

proc freq data=out noprint;
table label*_INTO_/out=table;
run;

proc sql;
create table _tmp as
select &i as cv, &k as K, 31 as dim, sum(PERCENT) as correct
from table
where label=_INTO_
;
quit;
%let t1=%sysfunc(time());
data _null_;
dt=round(&t1-&t0, 0.1);
call symput('dt', dt);
run;

%put k=&k i=&i time=&dt;
%if %eval(&k+&i)=3 %then %do;
data result; set _tmp; run;
%end;
%else %do;
data result; set result _tmp; run;
%end;
%end;
%end;
options source notes;
%mend;
options nosource nonotes;
%loopy;
options source notes;

proc means data=result noprint;
by k;
var correct;
output out=result_summary mean(correct)=mean std(correct)=std;
run;

Figure 3. Correct% vs. Number of Features by Number of Neighbors

Now we see that the best possible number of features may exist somewhere between 21 to 41 while the number of neighbors has a much smaller influence except for using only 2 NN. Now we are in a good situation to conduct more greedy search over the feature space. We are going to search through all combinations of number of features from 21 to 41 at single interval and neighbors 3 to 7, which results in a total number of 120 iterations each runs over the 10 folds, and the code ran about 3 hours on my PC. The code below executes this exercise and the result is summarized in table 1. We can see that using up to 31 features and 3 nearest neighbors will be the best choice for the data at hand.

%macro loopy;
%do k=2 %to 7;
%do d=21 %to 41;
%do i=1 %to 5;
%let t0=%sysfunc(time());
options nosource nonotes;
proc discrim data=allpcastd(where=(group='T' & cv^=&i ))
method=npar k=&k noprint
testdata=allpcastd(where=(group='T'))
testout=out(drop=Prin:)
;
class label;
var Prin1-Prin&d;
run;

proc freq data=out noprint;
table label*_INTO_/out=table;
run;

proc sql;
create table _tmp as
select &i as cv, &k as K, &d as dim, sum(PERCENT) as correct
from table
where label=_INTO_
;
quit;

%let t1=%sysfunc(time());
data _null_;
dt=round(&t1-&t0, 0.1);
call symput('dt', dt);
run;

%put k=&k d=&d i=&i time=&dt;
%if %eval(&k=2) & %eval(&d=21) %then %do;
data result; set _tmp; run;
%end;
%else %do;
data result; set result _tmp; run;
%end;
%end;
%end;
%end;
options source notes;
%mend;
options nosource nonotes;
%loopy;
options source notes;

proc means data=result noprint nway;
class k dim;
var correct;
output out=result_summary2 mean(correct)=mean std(correct)=std;
run;

proc means data=result_summary2 nway noprint;
class k;
var mean;
output out=result_max idgroup(max(mean) out[1](mean std dim)=mean std dim);
run;

proc sgplot data=result ;
where i=.;
dot dim /response=correct stat=mean limitstat=stddev numstd=1 group=k;
label dim="# of Features";
run;

Figure 4. Correct% and Confidence Interval by # of features and # of neighbors.

Table 1. Tuning Parameters Selection

Now, re-score the data using chosen tuning parameters: up to 33 features and 3-NN, we end up an accuracy of 96.80% on the leader-board, which certainly is not the best result you can get from kNN method, but easily beat both the 1NN and ramdomForest benchmark results.

With a little bit more tweak, the accuracy from a kNN with similar parameters can be boosted to 97.7%.

KNN Classification and Regression in SAS

2012-11-25T13:47:00.004-08:00

PDF available at here. Related post on KNN classification using SAS is here.

In data mining and predictive modeling, it refers to a memory-based (or instance-based) algorithm for classification and regression problems. It is a widely used algorithm with many successfully applications in medical research, business applications, etc. In fact, according to Google Analytics, it is the second most viewed article on my SAS programming blog of all time, with more than 2400 views in past year.

In classification problems, the label of potential objects is determined by the labels of closest training data points in the feature space. The determination process is either through "majority voting" or "averaging", potentially with weighting. In "majority voting", the label of object is assigned to be the label which most frequent among the k closest training examples. In "averaging", the object is not assigned a label, but instead, the ratio of each class among the k closest training data points. In Regression problems, the property of object is obtained via a similar "averaging" process, where the value of the object is the average value of the k closest training points.

In practice, both the "majority voting" and "averaging" process can be refined by adding weights to the k closest training points, where the weights are proportional to the distance between object and the training point. In this way, the closest points will have the biggest influence on the final results.

In fact, the SAS implementation of kNN classification has the averaging process be weighted by volume size. Interested readers are encouraged to read the manual for details. kNN algorithm has several characteristics that worth addressing. Understanding its advantages and disadvantages in theory and in practice will help practicioners better leverage this tool. Generally, KNN algorithm has 4 advantages worth noting. It is simple to implement.

Theoretically, kNN algorithm is very simple to implement. The naive version of the algorithm is easy to implement. For every data point in the test sample, directly computing the desired distances to all stored vectors, and choose those shortest k examples among stored vectors. It is, however, computationally intensive, especially when the size of the training set grows. Over the years, many nearest neighbor search algorithms have been proposed, seeking to reduce the number of distance evaluations actually performed. Using an appropriate nearest neighbor search algorithm makes k-NN computationally feasible even for large data sets.

In SAS, the k-d Tree data structure and associated search algorithm of Friedman et. al. is implemented. It is analytically tractable. The error rate for 1-NN classification is twice as much as Bayesian Error in asymptotic sense. It is highly adaptive to local information. kNN algorithm uses the closest data points for estimation, therefore it is able to take full advantage of local information and form highly nonlinear, highly adaptive decision boundaries for each data point. It is easily implemented in parallel. Because it is instance based, for each data point to be scored, the algorithm check against the training table for the k nearest neighbor. Since each data point is independent of the others, the execution of search and score can be conducted in parallel.

On the other hand, kNN algorithm has 3 disadvantages, which will also be addressed in details in section PRACTICAL ISSUES ON IMPLEMENTATION. It is computationally very intensive, because for each records to be classified and/or scored, the program has to scan through available data points to find the nearest neighbors and then determine which ones to use for further calculation. It is demanding in storage, because all records in training sample have to be stored and used whenever scoring is necessary. It is highly susceptible to curse of dimensionality, because as dimension increases, the distance calculation and finding the appropriate space to contain at least one relevant record is becoming difficult.

In order to conduct either KNN Classification or KNN regression in SAS, there are two basic approaches. The first one utilize special options in certain SAS Procedures and conduct the KNN analysis directly. I call it "Direct Approach"; while the second method will first find the nearest data points explicitly leveraging search capability in some SAS procedures and manually calculate the final response values, which is simple anyway. I call it "Manual Approach". To take the manual approach, there are three PROCs that are capable to find nearest neighbor points in efficient way, namely PROC FASTCLUS, PROC KRIGE2D and PROC MODECLUS. Each has its pro and con and all need data preparation. The discussion of this issue will beyond the scope of this paper. To take the direct approach, PROC DISCRIM can be directly used for KNN Classification while PROC KRIGE2D can be directly used for KNN regression, even though the functionality SAS provides is only basic. For kNN-based classification, PROC DISCRIM offers intuitive programming interface and rich options. For kNN-based regression, even though there is not dedicated procedure available, PROC KRIGE2D can be used to fulfill this task, with some tweaks. In this section, we will demonstrate how to conduct these two data mining tasks in SAS and address some closely related issues. For kNN-based regression, we will discuss the details of how to tweak PROC KRIGE2D. SAS offers kNN classification capability via PROC DISCRIM, by invoking the following procedure options:
METHOD=NPAR K=

The METHOD=NPAR option asks SAS to use non-parametric discrimination function, together with K=' option PROC DISCRIM will use kNN classification, where K= tells SAS how many neighbors to use in determining the . A typical SAS program for kNN-based classification looks like the follows:

proc discrim data    = train ................................[1]
              method = npar k=5 .............................[2]
        test    = toscore ..............................[3]
              testout = toscore_out...........................[4]
;
      class y; ...............................................[5]
var x1-x10; ............................................[6]
run;

We explain the function of each statement below.

Statement [1] tells SAS that \verb|PROC DISCRIM| should be called to process the data named \emph{train}.

Statement [2] tells SAS to apply kNN Classification method using 5 nearest neighbors.

Statement [3] tells SAS that the classification rule be applied to a test data called 'toscore'.

Statement [4] tells SAS to output classification result for the test data and to name the output data as \emph{toscoreout}.

Statement [5] defines Y as the dependent variable to use in \verb|PROC DISCRIM|.

Statement [6] tells SAS to use x1 to x10 as input features in calculating the distance to neighbors.

kNN classification is computationally intensive. Based on SAS manual [4], for the specific tree search algorithm implemented in \verb|PROC DISCRIM|, the time usage, excluding I/O time, is roughly proportional to $$\log(N)*(N*P)$$, where N is the number of observations and P is the number of variables used. kNN is a memory-based method, when an analyst wants to score the test data or new data in production, the whole raw table will be loaded to memory in the scoring process. In computing, \verb|PROC DISCRIM| uses memory approximately proportional to the second order of number of variables, i.e. proportional to $$P^2$$. Therefore, it would be interesting to see how the cost of computing, in terms of memory usage and time consumption, increases as the number of observations in training data and the number of features are increasing. \subsection{kNN Regression using SAS} kNN technique can be applied to regression problems, too, but the coding in SAS is not as straightforward as in a classification problem. kNN regression uses the average value of dependent variable over the selected nearest neighbors to generate predicted value for scoring data point. In SAS, we can use \verb|PROC KRIGE2D| to conduct this analysis, but it needs more cautions and tweaks. The following example code instructs SAS to do a kNN regression on variable Y using 5 nearest neighbors and 2 features: X1 and X2.

 
proc krige2d data=train outest=pred outn=NNlist; ..............[1]

coordinate xc=X1 yc=X2; ..........................................[2]
predict var=Y NP=5; .................................................[3]
model scale=10 range=1E-7 form=exp; ......................[4]
grid griddata=toscore xc=X1 yc=X2; ...........................[5]

 run;

Statement [1] invokes \verb|KRIGE2D| procedure and ask SAS to operate on data named \emph{train}, and output the prediction to data named \emph{pred} via option OUTEST=. The option OUTN= ask SAS to output information about the selected nearest neighbors of each data points in the data specified in GRID statement. KRIGE2D will directly output the prediction of dependent variable in the data from 'OUTEST= ', but only output the corresponding value of dependent variable in the selected closest data points as well as the coordinates. This is a disadvantage in the sense that when cross validation is conducted, the data from 'OUTN=' does not contain the actual value of dependent variable and this data set has to be merged back to the validation part to obtain actual values so that errors can be computed.

Statement [2] tells \verb|KRIGE2D| that the names of the 2 coordinates to be used, or in the language of data mining, the 2 features to be used. In this example, it is X1 and X2.

Statement [3] tells \texttt{KRIGE2D} to use 5 nearest neighbors to predicts dependent variables's value at current grid position. As to be discussed, the number of nearest neighbors is a tuning parameter in this algorithm, and should be determined using data driven approaches, such as Cross Validation.

Statement [4] tells \texttt{KRIGE2D} what model to be used in kriging using 'FORM=' option and the associated parameters, scale and range. Note that, in order to conduct a kNN regression, it is required to specify a large scale value and a small range value. For example, in example above, scale is set to 10 while range is set to 1E-7. Typically, a scale value larger than 1 and a range value smaller than 1E-7 work well. The type of model is irrelevant. The reason why this works is beyond the scope of this paper and we will not discuss on this issue here.

Statement [5] tells \texttt{KRIGE2D} that the name of the data to be scored is \emph{tosore} in option \texttt{GRIDDATA=} and the corresponding names of features (coordinates) to be used for scoring.

\verb|PROC KRIGE2D| was designed to perform ordinary kriging in two dimensions, therefore it is not fully functional in the standard kNN regression practice. There are two major issues with \verb|PROC KRIGE2D| when it is applied in a kNN regression analysis.

The first problem is that \verb|PROC KRIGE2D| accept two and only two inputs as features because it treat kNN regression as a two dimensional kriging. The second problem is due to the nature of semivariogram model used in kriging. In cases where the data are dense and highly overlapped, \verb|PROC KRIGE2D| will report singularity in kriging system. However, these issues should not bother analysis in practice and there are good reasons which will be detailed below.

Overcome Fixed Number of Features in KRIGE2D:

In a real kNN regression application, the number of features ranges from 1 to hundreds if necessary, but \verb|PROC KRIGE2D| requires exactly two features as inputs. When only 1 feature is used, or when more than 2 features are required, the analysts can take a detour by pre-processing the data up front. We will discuss these two cases separately below.

Only one feature. When only one feature is used, we can set up a dummy input, which is a constant number, in both the training data and scoring data. The dummy variable will cancel out in calculating, thus it has no effect virtually. The code below shows a demonstration:

data trainv/view=trainv;
set train;
dummy=1;
run;
data toscorev/view=toscorev;
set toscore;
dummy=1;
run;
proc krige2d data=trainv outest=pred outn=knn;
coordinate xc=dummy yc=feature1;
predict var=Y NP=5;
model scale=10 range=1E-7 form=exp;
grid gdata=toscorev xc=dummy yc=feature1;
run;

In the example code above, we used SAS data view to avoid re-write the data sets and save computing resources. This technique is explained in \cite{view}.

More than two features. When more than two features are used, there are two approaches to handle this problem. The first approach conduct dimension reduction upfront and \verb|PROC KRIGE2D| take two derived new dimensions as features to perform kNN regression; the second approach follows the ensemble principle \cite{ESL}, and randomly select two features out of the pool of available features. The final kNN regression prediction is an average of predictions from each of the ensemble regressions. The first approach can also use the ensemble principle, too, in practice. There are many dimension reduction techniques, and among all, the most widely used is Singular Value Decomposition (SVD), or its equivalent counterpart Principal Components Analysis (PCA). Details of SVD or PCA is beyond the scope of this paper and interested readers can refer to \cite{PCA}, \cite{Matrix}.

The example code in Appendix 1 demonstrates the first method where we conduct dimension reduction using PCA upfront and use the leading two orthogonal bases for kNN regression. Note that in the example code, we first combine both train and score data sets together because PCA is sample sensitive \cite{PCA}, and using records from both training data and scoring data will stabilize the generated orthogonal bases.

The second method is demonstrated below. It is worth noting that when the number of features is small, say 10, it is a good practice to select all pairwise combination as input features and ensemble all results at the final stage. On the other hand, when the number of such combinations blows as the number of features increases, a random sample of pairwise combination works fine, just like in a RandomForest algorithm. In the code below, we assume 100 features which are stored in a table called FEATURES and the ensemble iterates 50 times, each with a pair of randomly selected features.



%macro knnreg(train_dsn, score_dsn, feature1 , feature2 , loopid);

    proc krige2d data=&train_dsn outest=pred&loopid outn=nn&loopid;

            coord xc=&feature1 yc=&feature2;

            pred var=Y NP=10;

            model scale=2 range=1E-8 form=exp;

            grid gdata=&score_dsn xc=&feature1 yc=&feature2;

    run;

%mend;

%let nloops=50;

proc surveyselect data=features out=feature2 reps=&nloops sampsize=2 method=srs;

run;

data _null_;

        set feature2;

        by replicate;

        retain train_dsn 'train' score_dsn 'toscore';

        retain _id 0; array _f{2} $ _temporary_;

        if first.replicate then do;

           call missing(of _f[*]); 

           _id=1;

        end;

        _f[_id]=varname; _id+1;

       if last.replicate then do;

          call execute('%knnreg(' || compress(train_dsn) || ','

                                                  || compress(score_dsn) || ','

                                                  || compress(_f[1]) || ','

                                                  || compress(_f[2]) || ')' );

       end;

run;

The 'call execute' statement calls the pre-defined macro \%kddreg within a data step. For details about 'call execute' command, consult \cite{CallExecute}.

Singularity in kriging process. When the scoring point is extremely close to at least one of the training record, local ordinary kriging process will be singular and the estimation for the scoring point will be skipped. One solution is to set a large number for scale and small number for range in MODEL statement, for example:

MODEL SCALE=10 RANGE=1e-8 FORM=EXP;

In the extreme cases, while inconvenient, it is necessary to manually process the data containing nearest neighboring points from OUTN= to generate prediction. It turns out to be fairly simple, because ordinary kNN regression only requires average every k observations. This technique is demonstrated in Appendix 2. In such case, it is unnecessary to output predictions from KRIGE2D directly, and the output can be suppressed by specifying OUTEST=_NULL_

Practical issues on implementation. kNN is a proven technique in many classification and regression problems. To fully appreciate its power in real analysis, there are a couple of issues need to be addressed appropriately. These problems touch the pain points of efficiency, accuracy, robustness and more. In the sections below, we will first discuss efficiency around storage and computation issues and see how parallelization can improve the efficiency. Then we will discuss model selection that balances accuracy and robustness. Lastly, we will discuss other miscellaneous issues when applying KNN in practice.

Optimize Storage and Computation. kNN algorithm is expensive in both storage and computation. Therefore, it is important to the factors that contribute to complexity in time and space of this algorithm, thus optimize the implementation in SAS.

Optimize Storage. kNN classification requires a lot of storage because this is a in-memory algorithm, and all the training information to score a new data will be load into the memory to build a kd-Tree when conduct scoring. Therefore, when the size of the training data is large, and when the number of variables used is large, the required memory storage will be non-trivial. According to SAS/STAT Manual, the required storage for a kNN classification is at least c(32v+3l+128)+8v^2+104v+4l ; bytes, where v is the number of variables, c is the number of classes levels of the dependent variable, and l is the length of CLASS variable. Apparently, the storage requirement is linear in the latter two factors which are given by the problem at hand, while quadratic in the number of variables. To optimize storage, it is necessary to reduce the number of variables involved in the analysis.

1. Only choose the necessary variables in kNN analysis. In many cases, 2 or 3 features is enough, actually. On the other hand, like in RandomForest algorithm, multiple runs of kNN algorithm were conducted, each with a small subset of variables, and these classifiers or regression results will be ensembeled to be the final solution.

2. When selected number of variables is still relatively large, one commonly used technique for dimension reduction is Singular Value Decomposition (SVD), or equivalently, the Principal Component Analysis (PCA). Because kNN algorithm is highly adaptive to local information, SVD or PCA won't suffer from common critics that the information in dependent variable is not used when variable reduction is used. This technique is demonstrated in the Digits Recognition example below.

kNN algorithm used in classification and regression suffers a lot from curse of dimensionality, therefore, choosing only necessary and more useful features for the problem in hand is critical for the success of application and implementation.

Optimize Computation. kNN classification is also computationally intensive. Unlike many classification methods, such as logistic regression, discriminant analysis and other parametric regressions, kNN has no condense mathematical form to represent the information in training data, and lack of explicit model training. When scoring is needed, the training data is simply used to populate the sample search space with known results, and scoring is no different from a searching and comparison process, which is very time consuming. In the most naive way of searching, the algorithm compares scoring current data point to all data in the training sample, calculate desired distance metrics and select the top k nearest neighbors for their response information. This is computationally very redundant and intensive and doesn't leverage the distribution information in training sample well. Modern algorithm relies on heuristic searching algorithms and one of the most significant one is kd-Tree \cite{kd-tree}. kd-Tree stands for k-dimensional tree, and it is a space-partitioning data structure to organize data points in k-dimensional space for efficient search. In SAS, PROC DISCRIM ; stores data in kd-Tree structure. The time required to organize the observations into this specific tree structure is proportional to $$vn\ln(n)$$. The time for performing each tree search is proportional to $$\ln(n)$$. Apparently, the way to reduce computational time is reduce the number of data points n. On the other hand, it is advised to keep enough data points in training sample. This is because a k-d trees are not suitable for efficiently finding the nearest neighbor in high dimensional spaces. As a general rule, if the dimensionality is k, the number of points in the data, N, should be N >> 2k. Otherwise, when k-d trees are used with high-dimensional data, most of the points in the tree will be evaluated and the efficiency is no better than exhaustive search, and approximate nearest-neighbor methods are used instead \cite{SAS-STAT}. \subsection{Select the Optimal Values for Tuning Parameters} kNN algorithms have two basic tuning parameters, the number of nearest neighbors k, and the dimension of feature space d. The typical way to find the right values for these two parameters is brutal force computation, where program searches through a grid of different combinations of these two parameters and decide the best set based on cross validation results. The k in kNN refers to the number of closest data points, which is an important tuning parameter in kNN algorithm. There is no set rule for choosing the right k value, but usually when the data is dense, k is set to be small while when data is sparse, a relatively higher k should be used to reduce the variance. In many cases, the dimension is pre-determined. However, kNN is sensitive to curse of dimensionality and more often than not, the features are derived values, such as Principle Component Scores. The analysts have to find which subset of features is enough for the problem at hand. According to Friedman, et al, an objective data driven approach to choose the right value for k is the way to go. Typically, a Cross-Validation (CV) based approach is used, usually 5-fold CV or 10-fold CV. When using an M-fold CV to choose the tuning parameter k and d, it is customized to check the mean value and standard deviation of certain model accuracy measurement. For example, for Classification problem, both error rate and AUC can be used; for regression problems, RMSE is a good candidate measurement. The code in the two examples below illustrate these techniques.

Parallelization. kNN classification and regression can be easily parallelized. The most significant applications are in Cross-validation based tuning parameter evaluation and scoring. In Cross-Validation process, the analyst is able to open M concurrent sessions, each overs mutually exclusive set of tuning parameters. For example, the simplest case is issuing 2 concurrent SAS sessions, one examines k=1, ..., 5 while the other one examines k=6, ... , 11. The actual implementation follows the idea of \cite{MultiThreading}. In scoring, because each data point is assumed to be independent from the others, the scoring program can split the scoring data set into several pieces, depending on the number of available CPUs and the master control module will simultaneously call the scoring module which only operates on each piece at a time.

kNN in Action. In this section, we demonstrate kNN algorithms in real applications. Both examples will go through steps such as data introduction, data analysis and modeling using kNN algorithm, as well as selecting the tuning parameter k. The first example examines handwritten digits recognition, a 10-label classification problem. The second example examines a regression problem where the researcher wanted to estimate the phenolic records.

kNN Classification in Optical Digits. One successful application of kNN classifiers is in digits recognition projects, see \cite{LeCun1998} for examples. In this example, I will demonstrate a complete data analysis process using kNN classifier to identity hand written digits from Optical Digits images. The data is from the pre-processed Optical Digits images available at UC Irving Machine Learning data set repository and can be downloaded at the link here .

The data is extracted using the NIST preprocessing routines \cite{NIST} from normalized bitmaps of handwritten digits from a preprinted form. The original bitmaps are 32X32 but divided into non-overlapping blocks of 4X4 and the the number of on pixels are counted in each block. This generates an input matrix of 8x8 where each element is an integer in the range 0..16. This reduces dimensionality to 64 instead of 1024, and gives invariance to small distortions. In this project, the following steps are performed to fine tune the kNN classifier.

First, the training and testing data in CSV format are imported at \emph{Step 1}.

Then a baseline 1-NN classifier using all 64 features is built for comparison purpose at \emph{Step 2}. This baseline classifier return an error rate of 3.73\%.

In \emph{Step 3}, a fine tuned kNN classifier is built, where the tuning parameter k, number of closest points, and d, number features, are selected based on 10-fold Cross Validation.

Besides, instead of using the original 64 features, a Singular Value Decomposition (SVD) was conducted on the training and testing data sets combined to obtain derived features. This data pre-processing technique brings in several advantages.

First, the derived features are ranked by how much variance in the feature space each explained, so that the selection of features becomes simpler using CV. Second, SVD changes original sparse feature space into a dense one, and in general fewer closest data points are required to obtain optimal performance in this case.

In \emph{Step 4}, the selected kNN classifier, with k=4 closest data points and d=36 derived features, is applied to the testing data set, obtaining an error rate of 1.89\%.

As an alternate to 10-fold CV, in \emph{Step 5}, Leave-One-Out CV was used to fine tune the kNN classifier, resulting a set of parameters of k=6 and d=4.

In \emph{Step 6}, the new classifier was applied to the testing data set, resulting in an error rate of 1.4\%. \end{enumerate} The 10-fold Cross Validation process is pretty computationally intensive even for such a small data set. Scanning across 900 different combinations of k, d pairs took almost 25 minutes on a main stream PC. The pay off worth it, however, as the error rate drops to 1.89\% from 3.73\%, a 49\% improvement. On the other hand, using Leave-One-Out Cross Validation process took 22 minutes on the same PC, and the selected parameter pair is k=6, d=38, with an error rate of 2.06\% from the same scoring data set, which is not as good as the result from 10-fold CV but still way better than baseline result.

The 10-fold CV is preferred to Leave-One-Out cross validation from pure computational perspective, as the computing time for the latter one increases quadratically and will become infeasible when the number of observation is more than 100K. Figure \ref{fig:C1} shows the error rate as a function of k and d, while figure \ref{fig:C2} zooms in for a closer look at area where number of features is between 20 and 40 and closest data points is smaller than 5. We can see that globally, the error rate is a relatively smooth function of number of features over all and the number of neighbors does not influence the error rate as much as the number of features does, while locally there is significant variation as shown in figure \ref{fig:C2}. This implies that one heuristic way to scan for the best possible numbers of features is bi-section searching.

The complete SAS code that conducts brute force searching for tuning parameters can be found in Appendix 1. kNN Regression Example. As a demonstration of kNN regression, we used the \emph{Galaxy} data set from the book \emph{Elements of Statistical Learning 2}, which is available at \cite{galaxy}. The galaxy data was used to demonstrate smoothing techniques in the book and was visualized in figure 6.8, where local linear regression with nearest neighbor window of 15\% of the data in $\textrm{R}^2$ was applied to predict velocity data of galaxies.

In our example, a k-Nearest Neighbor prediction is applied, where two values of k are used: k=5 and k=48. k=48 is equivalent to using 15\% of the data as neighbor while k=5 is arbitrarily chosen to demonstrate the fact that smaller number of neighbors will produce more accurate but higher variance prediction while larger value of k will product less accurate but smaller variance hence smoother prediction.

Combine the two predictions in many cases may product better results among all. In figure \ref{fig:galaxy}, 3-D views of actual data and smoothed surface using different smoothing parameters are shown. The upper left panel shows the raw data. Obviously, the galaxy has a star like spatial shape with a total of about 10 arms extending outside from the center. The upper left panel shows smoothed 3-D surface using 5-nearest neighbors while the lower left panel shows smoothed 3-D surface using 48-nearest neighbors. It is noticeable that using less neighbors produce a rougher surface but adaptive better to local where there are more data, such as along the galaxy arms, but in more empty spaces, the prediction is wilder. On the other hand, using 48 neighbors produces much smoother surface yet the galaxy arms are less visible. The lower right panel shows the average of these two surface maps. This averaged out map keeps the galaxy arms clearly visible yet smooth out empty space pretty well.

The complete SAS code is in Appendix 2.

kNN algorithm is an old but very powerful technique if used appropriately. It has been shown that both kNN Classification and kNN regression can be carried out within SAS/STAT very easily. The functions are basic but more than enough to handle many real world problems. Both the theory and practice are discussed, with example SAS code demonstrated. The discussion here, however, is only the tip of iceberg. There are many advanced topics that can't be fitted into this short article. These topics include but not limited to: customized distance function, tangent invariance, etc. Interested readers are encouraged to read the related books and papers.

Appendix 1.

/*@@@@ --- Step 1 --- @@@@*/
libname digits "C:\_Data\Digits";

data digits.train;
infile "C:\_Data\Digits\optdigits.tra" lrecl=512 dlm=',';
input x1-x64 y;
run;

data digits.test;
infile "C:\_Data\Digits\optdigits.tes" lrecl=512 dlm=',';
input x1-x64 y;
run;

/*@@@@ --- Step 2 --- @@@@*/
proc discrim data=digits.train testdata=digits.test
testout=testout
method=npar k=1 noprint;
class y;
var x1-x64;
run;

proc freq data=testout;
table y*_INTO_;
run;

/*@@@@ --- Step 3 --- @@@@*/
data all;
set digits.train(in=_1) digits.test(in=_2);
if _1 then Group='T';
else Group='S';
run;

proc princomp data=all out=pca noprint noint cov;
var x1-x64;
run;
proc standard data=pca out=pcaout mean=0 std=1;
var Prin1-Prin62;
run;

/*@@@@@ 10-fold CV approach @@@@@@@@@@*/
data pcaout;
set pcaout;
array _xp{10} p1-p10 (10*0.1);
cv=rantbl(2342342, of _xp[*]);
run;
proc freq data=pcaout; table cv; run;

%macro loopy;
ods select none;
options nonotes nosource;
%do k=1 %to 15;
%do n=2 %to 60;
%put k=&k n=&n;
%do cv=1 %to 10;
proc discrim data=pcaout(where=(group='T' & cv^=&cv))
testdata=pcaout(where=(group='T' & cv=&cv))
testout=testout
method=npar k=&k noprint;
class y;
var Prin1 - Prin&n;
run;
proc freq data=testout;
table y*_INTO_ /out=result;
run;
proc sql;
create table _tmp&cv as
select sum(COUNT) as total,
sum((y=_INTO_)*COUNT) as correct,
&k as k, &n as n,
1-calculated correct / calculated total as error format percent7.3
from result
;
quit;
%end;
data _tmp;
set %do cv=1 %to 10; _tmp&cv %end;;
run;
proc means data=_tmp noprint;
var error;
output out=_tmp mean(error)=error std(error)=error_std;
run;
data _tmp;
length k n 8;
set _tmp;
k=&k; n=&n; drop _TYPE_;
run;
%if %eval(&k+&n=3) %then %do;
data iter; set _tmp; run;
%end;
%else %do;
proc append base=iter data=_tmp force; run;
%end;

%end;
%end;
ods select all; options notes source;
%mend;

dm log 'clear';
%let time0=%sysfunc(time());
%loopy;
%let time1=%sysfunc(time());

data _null_;
dt=(&time1-&time0)/60; put dt= 'Min';
run;

data iter;
set iter; format error percent7.4;
run;

proc transpose data=iter out=itert name=n prefix=n;
by k;
var error;
id n;
run;

data itert;
set itert;
array _n{*} n2-n60;
min=min(of _n[*]); max=max(of _n[*]); format min max percent7.4;
run;

/*@@@@ --- Step 4 --- @@@@*/
/*@@@@ use selected parameter to score the testing data @@@@*/
proc discrim data=pcaout(where=(group='T')) testdata=pcaout(where=(group='S'))
testout=testout(drop=x: Prin:)
method=npar k=4 noprint;
class y;
var Prin1-Prin36;
run;
proc freq data=testout;
table y*_INTO_ /out=result;
run;
proc sql;
create table temp as
select sum(COUNT) as total,
sum((y=_INTO_)*COUNT) as correct,
4 as k, 36 as n,
1-calculated correct / calculated total as error format percent7.3
from result
;
quit;

/*@@@@ --- Step 5 --- @@@@*/
/*@@@@@ Use leave-one-out cross-validation approach @@@@@@*/
%macro loop;
%do k=1 %to 15;
%do i=1 %to 62;
proc discrim data=pcaout(where=(group='T'))
method=npar k=&k noprint crossvalidate
outcross=outcross out=out;
class y;
var Prin1-Prin&i;
run;

proc freq data=outcross noprint;
table y*_INTO_ /out=table;
run;

proc sql;
create table _tmp as
select &k as K, &i as I, sum(PERCENT) as correct
from table
where y=_INTO_
;
quit;
%if %eval(&i+&k)=2 %then %do;
data result; set _tmp; run;
%end;
%else %do;
proc append base=result data=_tmp; run;
%end;
%end;
%end;
%mend;

options nonotes nosource;
%let time0=%sysfunc(time());
dm log 'clear';
%loop;
options notes source;
%let time1=%sysfunc(time());
data _null_;
dt=(&time1-&time0)/60;
put dt= 'Min';
run;

proc transpose data=result out=result_t name=i prefix=n;
by k;
var correct;
id i;
run;

data result_t;
set result_t;
array _n{*} n2-n62;
min=min(of _n[*]); max=max(of _n[*]);
loc_min=whichN(min, of _n[*]); loc_max=whichN(max, of _n[*]);
format min max best.;
run;

proc discrim data=pcaout(where=(group='T'))
testdata=pcaout(where=(group='S'))
testout=testout(drop=x: Prin:)
method=npar k=6 noprint;
class y;
var Prin1-Prin38;
run;
proc freq data=testout;
table y*_INTO_ /out=result;
run;
proc sql;
create table temp as
select sum(COUNT) as total,
sum((y=_INTO_)*COUNT) as correct,
6 as k, 38 as n,
1-calculated correct / calculated total as error format percent7.3
from result
;
quit;

Appendix 2.

proc import datafile="c:\galaxy.csv" dbms=csv out=galaxy replace; run;
run;

proc g3d data=galaxy;
scatter east_west*north_south=velocity /rotate=25 size=1 noneedle tilt=45;
format _numeric_ 8.2;
run;
quit;

%macro knnreg(k);
ods select none;
proc krige2D data=galaxy outest=knn_pred&k outn=neighbor;
coordinates xc=east_west yc=north_south;
predict var=velocity NP=&k;
model scale=10 range=1E-8 form=exp;
grid gdata=a(drop=velocity) xc=east_west yc=north_south ;
run;
ods select all;

proc g3d data=knn_pred;
plot gxc*gyc=estimate / grid rotate=25 tilt=45
ctop=cx993366
cbottom=cx9999ff
caxis=black;
label gxc="east_west"
gyc="north_south";
format _numeric_ 8.2;
run;

%mend;
%knnreg(5);
%knnreg(48);

data knn_pred;
merge knn_pred5(rename=(estimate=estimate5))
knn_pred48(keep=Estimate rename=(estimate=estimate48));
estimate=(estimate5+estimate48)*0.5;
run;

proc g3d data=knn_pred;
plot gxc*gyc=estimate / grid rotate=25 tilt=45
ctop=cx993366
cbottom=cx9999ff
caxis=black;
label gxc="east_west"
gyc="north_south";
format _numeric_ 8.2;
run;

Finite Mixture Model for Loss Given Default (LGD)

2012-10-04T08:24:00.000-07:00

Loss Given Default (LGD) is a key business metric of risk in financial service. One unique feature of this metric is overdispersion and the other is multi-mode.

Finite mixture model is an effective way to accommodate both. Multi-mode refers to the case where the distribution has more than one peak. Over-dispersion refers to the case where the dispersion of the data is more than what the assumed model implies. [3] provides a good introduction to finite mixture models.

Consider again the loss given default case discussed in the paper of Matt Flynn [1]. In this section, a complete modeling process is presented using finite mixture model that leverage the computing capability of PROC FMM in SAS v9.3 to model the variance in data that is not captured in a unified single model. Notice that because the dependent variable is bounded between 0 and 1, there are several options for modeling purpose. One can use a beta model, as discussed above; or a logit transformation can be applied first to the dependent variable (lgd) and a normal distribution is usually assumed to model the data. In this exhibition, the second approach is used.

First of all, a complete Exploratory Data Analysis (EDA) is conducted to study the distribution and relationship between and among dependent variable and covariates. Here, EDA consists of two exercises. At step one, the distribution and pairwise relationship of variables are studied. Then, non-parametric models are applied to each covariate to study the possible functional forms between dependent variable and independent variables.

The first figure shows the distribution of logit transformed Loss Given Default data. Apparently, it mimics normal distribution much closer than the raw distribution, but still it skews towards left rather than symmetric. Figure 2 below shows a matrix view of pairwise relationship between pairs of the four variables in the data

Figure 1. Histogram of Logit Transformed Loss Given Default Variable

Figure 2. Matrix Scatter Plot of All Variables.

In addition to a simple descriptive analysis using raw variables, it is often useful to conduct a non-parametric univariate regression so that the possible functional relationship between dependent variable and each covariate can be studied. Figure 3 shows the LOESS regression results of logit transformed LGD against the three covariates, namely, the annual average default rate, leverage coefficient by companies and industrial average default rate. It is noticeable that there is slightly quadratic relationship between logit(LGD) and average default rate by year and by industry. The linear relationship between logit(LGD) and leverage coefficient, however, is robust. On the other hand, we suspect that there is over-dispersion that can’t be captured adequately by these variables.

Figure 3. LOESS regression exhibition LGD against average default rate by year, leverage coefficient by firm and industrial average default rate.

In order to examine the sample and their relationship closer, an OLS model is built with both average default rates by year and by firm are modeled using quadratic relationship and leverage coefficient by firm is modeled as a linear term.

Figure 4. Diagnostics from polynomial OLS

Major model fitting diagnostic results are shown in figure 4. Key observations are:

1. The fit is reasonably well, with residuals approximately normally distributed.

2. There is noticeable over-dispersion when the residuals are plotted against predicted values. The variance is smaller at higher predicted values but larger at lower range.

3. Q-Q plot also indicates existence of over-dispersion.

As an alternative approach to a polynomial OLS, a finite mixture model is proposed. As a demonstration, we don’t go through the complete modeling process, especially the testing process to determine appropriate number of latent components; instead, a 2-component finite mixture model is used while only linear terms of the three covariates will enter the model. The purpose of this design is to show that even with a simpler and hence easier to interpret individual model structure, a finite mixture model is able to come up a better model that captures major stochastic effects in the data.

The following code builds a 2-component mixture model and examines the distribution of residuals.


ods html;
ods graphics on;
proc fmm data=lgddata2 plots(unpack);
     model lgd2 = lev lgd_a i_def /k=2 ;
     output out=fmm_out pred resid(component) resid(overall);
run;quit;
ods graphics off;
proc sgplot data=fmm_out;
     density resid_1 /type=kernel;
     density resid_2 /type=kernel;
     histogram resid ;
     density resid /type=normal;
run;

Figure 5. Distribution of residuals from a 2-component normal mixture model

There are two key features from PROC FMM that worth mention. First, the mixture of distributions can be from different families. For example, a normal distribution can be mixed with a T-distribution. Second, the probability each data point belongs to a latent class can be modeled with covariates, so to enhance the model interpretability.

The following code demonstrates a 2-component mixture model with one component is modeled as normally distributed while the other one is modeled as T-distributed.


/* mixture of heterogeneous distributions */

ods graphics on;

proc fmm data=lgddata2;

model lgd2 = lev lgd_a i_def /dist=normal;
     model      + lev lgd_a i_def /dist=t ;
     output out = fmm_out2 pred resid(component) resid(overall);
run;

ods graphics off;

The idea is simple. Outline a group of individual models of common distribution family in one MODEL statement, and the other groups of individual models of the same distribution family in other MODEL statement, without specifying the dependent variable but instead a ‘+’ symbol to indicate it is one layer add-up to existing individual models. Each group of individual models can be a mixture of multiple components. For example, in above code, if ‘K=2’ is specified in one of the MODEL statement, say in the first MODEL statement, then it tells SAS to model the data as a mixture of 2-component normal distribution and 1-component T distribution. This capability of modeling data with heterogeneous mixture distributions is a powerful tool in predictive modeling and advanced analytics. Of course, in this particular sample data, there is no immediate benefit by using a heterogeneous mixture distribution.

The second feature that gives PROC FMM extreme power in advanced analytics is the capability to model probability model of latent classes using covariates, so that analysts are able to study which factors will contribute to classifying data points into different latent classes. This feature gives analyst the power to interpret latent class with insightful factors and even provide better model fit. What an analyst needs to do is to add the following statement:


PROBMODEL &covars;

The probability model can be one of the following four: LOGITISTIC, PROBIT, LOG-LOG, and COMPLEMENTARY LOG-LOG. Please refer to SAS/STAT manual [2] for details.

Reference:
[1]. Matt Flynn, http://www.casact.org/education/spring/2011/handouts/C11-Flynn.pdf
[2]. SAS Institute, Inc., SAS/STAT v9.3 User's Guide
[3]. Geoffrey McLachlan, David Peel, Finite Mixture Models, John Wiley & Sons, 2000
Complete Code:


data lgddata;
   informat lgd  lev  12.9 lgd_a 6.4  i_def 4.3;
   input lgd lev lgd_a i_def;
   label lgd   = 'Real loss given default'
      lev   = 'Leverage coefficient by firm'
      lgd_a = 'Mean default rate by year'
      i_def = 'Mean default rate by industry';
   cards;
0.747573451 0.413989786 0.6261 1.415
0.99 0.413989786 0.6849 1.415
0.06581075 0.230361142 0.4566 1.183
0.351287992 0.541339309 0.6261 2.353
0.25844921 0.541339309 0.4566 2.353
0.01968009 0.812 0.6715 0.743
0.931035513 0.546732229 0.6715 2.353
0.341254925 0.71 0.6715 1.183
0.35075456 0.855339361 0.6715 2.353
0.045826764 0.313983237 0.6261 0.743
0.025754193 0.190648237 0.4566 0
0.759732568 0.490953756 0.6261 2.353
0.757989999 0.910788759 0.6261 1.415
0.6 0.336071518 0.6261 1.183
0.374480256 0.414862374 0.4566 0.967
0.168726407 0.612063995 0.6261 1.183
0.283909643 0.693928717 0.6261 2.353
0.747018382 0.937072431 0.6261 1.415
0.686300059 0.801162532 0.6715 2.353
0.050051313 0.365725066 0.4566 0
/* More data is not shown here to save space */
;
run;

/* dep var is between (0, 1), using logit transformation */
data lgddata2;
      set lgddata;
   lgd2=log(lgd/(1-lgd));
run;

/* check raw dep var distribution */
proc sgplot data=lgddata2;
      histogram lgd2;
   density   lgd2/type=normal;
run;

/* check relationship between raw dep and covariates */
proc sgscatter data=lgddata2;
      matrix  lgd2  lev lgd_a  i_def;
run;


/* study the relationship */
ods graphics on;
proc loess data=lgddata2;
      model lgd2 = lev;
run;

proc loess data=lgddata2;
     model lgd2 = lgd_a;
run;

proc loess data=lgddata2;
      model lgd2 = i_def;
run;
ods graphics off;

data lgddata2;
      set lgddata2;
   lgd_a2=lgd_a**2;
   i_def2=i_def**2/10;
run;

/* build baseline model */
ods graphics on;
proc reg data=lgddata2;
      model lgd2 = lev lgd_a i_def  lgd_a2   i_def2;
   output  out=reg_out  pred  resid;
run;quit;
ods graphics off;


ods html;
ods graphics on;
proc fmm   data=lgddata2  plots(unpack);
      model lgd2 = lev  lgd_a  i_def /k=2 ;
   *probmodel lev  lgd_a  i_def;
   output  out=fmm_out  pred  resid(component)  resid(overall);
   *bayes;
run;quit;
ods graphics off;

proc sgplot data=fmm_out;
      density   resid_1 /type=kernel;
   density   resid_2 /type=kernel;
   *density   resid_3 /type=kernel;
   histogram   resid ;
   density     resid /type=normal;
run;


/* mixture of heterogeneous distributions */
ods graphics on;
proc fmm  data=lgddata2;
     model   lgd2 = lev  lgd_a  i_def  /dist=normal;
  model       + lev  lgd_a  i_def  /dist=t  ;
  output  out=fmm_out2  pred resid(component)  resid(overall);
  PROBMODEL &covars;
run;
ods graphics off;


proc sgplot data=fmm_out2;
      density   resid_1 /type=kernel;
   density   resid_2 /type=kernel;
   *density   resid_3 /type=kernel;
   histogram   resid ;
   density     resid /type=normal;
run;

SAS functions for computing parameters in Erlang-C model

2012-07-12T12:51:00.001-07:00

Call center management is both Arts and Sciences. While driving moral and setting up strategies is more about Arts, staffing and servicing level configuration based on call load is in the domain of Sciences.

The science part of call center management is based on Queueing Theory, which studies "the Phenomena of standing, waiting and serving" [1], [2].The queueing model abstracts the complex phenomena of call center acticities into an idealized mathematical framework, so that it is possible to study the key characteristics and dynamics of the real world. There are many different models that could be applied to Call Center management, one simple yet popular model is called Erlang-C model. Some other highly related models are Erlang-B and Erlang-A models, but they won't be covered in this post.

In this article, Erlang-C model will be briefly described and the computation of key parameters of this model will be discussed, followed by a complete suite of SAS code pieces. It is interesting that Erlang-C model is rarely touched in the SAS community. Searching through SAS-L archives, there were only two discussions relevant, and they are set apart by 10 years [3], [4] ! The responders includes legendary David Cassell, SAS guru Daniel J. Nordlund, SAS computing expert Rick Wicklin (his SAS blog), and of course me.

Erlang-C model is a commonly used mathematical framework for many Call center staffing and management softwares. Based on Kendall notation, this model can be expressed as M/M/C/\inf, where the first M indicates a Poisson distribution for Arrival process, while the second M indicates an Exponential distribution for servicing time. Both 'M' refers markovian the 'C' indicates number of servers which is a non-negative integer, while the last \inf refers infinite number of places in the system. Erlang-C system has extra assumption that it is First In First Out (FIFO) with infinite calling population. Therefore, Erlang-C model describe a Call Center that serves customers in the order they arrive in, while customers arrives following a Poisson process and the servicing time is Exponentially distributed. When a customer is not immediately served, it is put on hold until next available server (CSR). The customers are assumed to be patient so that they won't drop out but instead wait how-long-it-takes until next CSR. The system can take infinite number of waiting customers.

In Erlang-C model, the traffic to a Call Center is measured in a so-called term of Offered Traffic Load, which is defined to be the ratio between the parameter \lambda of Poisson process and the parameter \mu of Exponential servicing time. Offered Traffic Load is therefore interpreted as the average number of services that are trying to happen simultaneously. It is free of time unit, but does has its own unit, called Erlangs.

At the core of Erlang-C model is the Erlang-C function. The Erlang C function computes the probability that an arriving customer in the Erlang C queueing model will find that all servers are busy. This is the same as the fraction of arriving customers that are delayed (i.e., must wait) before beginning service.

Erlang-C function is expressed as:

       P_c(n, x) = dPoisson(n, x)/(dPoisson(n, x) - (1-x/n)*pPoisson(n, x))

where dPoisson(n, x) is the density function of Poisson distribution with support 0,1, .., n and shape parameter x. pPoisson(n, x) is the corresponding cumulation distribution.

In Erlang-C function, there are three key elements:

    1. Probability of delay, P_c.
    2. # of servers, n.
    3. Offered Traffic Load, x.

Given any two, the last parameter can be calculated or approximated. For example, in the code pieces below:

    1. Function ErlcFractionDelay(Nsrv, traffErlang) calculates Probability of delay given number of servers and offered traffic Load in Erlangs units.

    2. Function ErlcNsrv(traffErlang, Erlc_max) estimates the minimum number servers to reach below the max probability of delay given by Erlc_max for a given offered traffic load of traffErlang.

    3. Function ErlcTrafFromFractionDelay(Nsrv, fractionDelay) Compute the max traffic load that Nserv servers can handle with queuing probability less or equal to fractionDelay.

In addition to the three key numbers in Erlang-C model, there are some other paremeters of particular interests to Call Center management. For example, the probability an arriving call will be put in queue for no more than certain waiting time given the average handling time. Function ErlcGoS(Nsrv, traffErlang, AHT, AWT) will calcualte this quantity.

The results of this set of functions has been calibrated against the data from [5], specifically Table 1. of [5] is reproduced exactly. The DATA STEP at the end of the program below will reproduce the results, and here is the result my own run.

For a complete list of quantities that Erlang-C model will work with and compute for, the webpage of Abstract Micro Systems @ http://abstractmicro.com/erlang/helppages/mod-c.htm provides a very good introduction.

Reference:
[1]. Kleinrock, Leonard (January 1975). Queueing Systems: Volume I – Theory. New York: Wiley Interscience. ISBN 978-0471491101.

[2]. Kleinrock, Leonard (April 1976). Queueing Systems: Volume II – Computer Applications. New York: Wiley Interscience. ISBN 978-0471491118.

[3]. http://www.listserv.uga.edu/cgi-bin/wa?A2=ind0210C&L=sas-l&P=R10763

[4]. http://www.listserv.uga.edu/cgi-bin/wa?A2=ind1206C&L=sas-l&P=R45472

[5]. Tibor Mišuth, Erik Chromý, Ivan Baroňák. Method for Fast Estimation of Contact Centre Parameters Using Erlang-C Model, in Proceedings of 2010 Third International Conference on Communication Theory, Reliability, and Quality of Service, DOI 10.1109/CTRQ.2010.38

All functions below have relatively detailed comments. The functions are provided as is without any implied and explicated warranties, and are for research purpose only.



proc fcmp outlib=sasuser.funcs.Erlc;   

     function ErlcFractionDelay(Nsrv, traffErlang); 
  /*
        Entry        : ErlcFractionDelay
     Purpose      : Compute the Probability of Delay when an arriving call is not
                    answered immediately and put in queue
     Usage        : P=ErlcFractionDelay(Nsrv, traffErlang)
     Inputs       : Nsrv=Number of Server (CSR)
                    traffErlang=Traffic in Erlang
     Handle Error : Return -1 if traffic in Erlang is larger than the number of 
                           servers;
                    Return -2 if either traffic in Erlang or Number of servers is 
                           negative;
                    Return -3 if number of servers is not integer;                          
  */
         if traffErlang>Nsrv then do;
            *put "Erlang amount should be smaller than Number of Servers";
   Ec=-1; return(Ec);
         end; 
   else if (traffErlang<0) or (Nsrv<0) then do;
            *put "Both Erlang amount and Number of Servers should be larger than 0";
   Ec=-2; return(Ec);
   end;
   else if abs(int(Nsrv)-Nsrv)>1E-8 then do;
            *put "Number of Servers should be integer";
   Ec=-3; return(Ec);
   end;  
         else do; 
            _p= PDF('POISSON', Nsrv, traffErlang);
         denom = _p
                  + (1-traffErlang/Nsrv)*CDF('POISSON', Nsrv-1, traffErlang);
            Ec = _p / denom;
            return (Ec);
   end;
     endsub;

  function ErlcGoS(Nsrv, traffErlang, AHT, AWT);
  /*
        Entry        : ErlcGoS
     Purpose      : Compute the Probability an arriving call will be put in
                    queue and wait for no more than AWT when the average
                    handling time is AHT. This term is also called Grade of 
                    Service (GoS). Note that the probability that an
                    arriving call will NOT be put in queue is:
                      P=1-ErlcFractionDelay();
     Usage        : GoS=ErlcGoS(Nsrv, traffErlang, AHT, AWT)
     Inputs       : Nsrv=Number of Server (CSR)
                    traffErlang=Traffic in Erlang
                    AHT=Average Handling Time
                    AWT=Average Waiting Time
     Handle Error : Return -1 if traffic in Erlang is larger than the number of 
                           servers;
                    Return -2 if either traffic in Erlang or Number of servers is 
                           negative;
                    Return -3 if number of servers is not integer;     
                       Return -4 if either of two time parameters is negative; 
  */
         if traffErlang>Nsrv then do;
            *put "Erlang amount should be smaller than Number of Servers";
      GoS=-1; return(GoS);
         end; 
   else if (traffErlang<0) or (Nsrv<0) then do;
            *put "Both Erlang amount and Number of Servers should be larger than 0";
   GoS=-2; return(GoS);
   end;
   else if  abs(int(Nsrv)-Nsrv)>1E-8 then do;
            *put "Number of Servers should be integer";
   GoS=-3; return(GoS);
   end;   
         else  if AHT<0 or AWT<0 then do;
            put "Both Average Holding Time and Average Waiting Time should be >0"; 
   GoS=-4; return(GoS);
   end;
   else do;
      GoS = 1-ErlcFractionDelay(Nsrv, traffErlang)*exp(-1*(Nsrv-traffErlang)*AWT/AHT);
      return(GoS);
    end;
  endsub;

  function ErlcNsrv(traffErlang, Erlc_max);
  /*
        Entry        : ErlcNsrv
     Purpose      : Compute the required number of servers (CSRs) to handle
                    given traffic load with queuing probability no more than
                    Erlc_max
     Usage        : N=ErlcNsrv(traffErlang, Erlc_max)
     Inputs       : traffErlang=Traffic load in Erlang
                    Erlc_max=max Queuing probability
     Handle Error : Return -1 if max queuing probability is not within (0, 1) 
                           range so to violate probability definition;
                    Return -2 if traffic in Erlang is negative;
  */

  /*REF:
       Algorithm is based on that from paper "Method for Fast 
             Estimation of Contact Centre Parameters Using 
             Erlang C Model" by Tibor Mišuth, Erik Chromý 
             and Ivan Baronák from Slovak University of 
             Technology. This paper appeared at "The Third 
             International Conference on Communication 
             Theory, Reliability, and Quality of Service", 
             2010
  */
        if Erlc_max>=1 or Erlc_max<=0 then do;
        *put "Max Erlang-C value should be between 0 and 1";
     N=-1; return(N);
     end;
     else if traffErlang<0 then do;
              *put "Traffic in Erlang should be >0";
     N=-2; return(N);
     end;
     else do;
           _k=floor(traffErlang)+1; s=0;
           do N=1 to _k; 
        s=(1+s)*N/traffErlang;
              end;
              _threshold=(1-Erlc_max)/((1-traffErlang/N)*Erlc_max);
              do while (s<=_threshold);
                 N=N+1;
        s=(1+s)*N/traffErlang;
        _threshold=(1-Erlc_max)/((1-traffErlang/N)*Erlc_max);
     end;
              return(N);
     end;
  endsub;

  function ErlcNsrvFromGoS( traffErlang, AHT, AWT, GoS);
  /*
        Entry        : ErlcNsrvFromGoS
     Purpose      : Compute the required number of servers (CSRs) to handle
                    given traffic load with given service level (GoS):
                           Average Handling Time /
                           the Probability an arriving call is put in queue
                           and wait no more than Average Waiting Time
     Usage        : N=ErlcNsrvFromGoS( traffErlang, AHT, AWT, GoS)
     Inputs       : traffErlang=Traffic load in Erlang
                    AHT=Average Handling Time
                    AWT=Average Waiting Time
                    GoS=probability an arriving Call in put in queue
                        but wait no more than AWT
     Handle Error : Return -1 if either time parameter is negative;
                    Return -2 if queuing probability is not in (0, 1) range;
                    Return -3 if traffic in Erlang is negative
  */
         if AHT<0 or AWT<0 then do;
            *put "Both Average Holding Time and Average Waiting Time should be >0"; 
      N0=-1; return(N0);
   end; 
   else if GoS<0 or GoS>1 then do;
      *put "GoS should be between 0 and 1";
   N0=-2; return(N0);
   end;
   else if traffErlang<0 then do;
            *put "Traffic in Erlang must be positive real number";
   N0=-3; return(N0);
   end;
         else do; 
            N0=floor(traffErlang)+1;
         _GoS=ErlcGoS(N0, traffErlang, AHT, AWT);
         do while (_GoS<=GoS);
            N0=N0+1;
            _GoS=ErlcGoS(N0, traffErlang, AHT, AWT);
         end;
         return(N0);
   end;
  endsub;

  function ErlcNsrvFromWait( traffErlang, AHT, AWT);
  /*
        Entry        : ErlcNsrvFromWait
     Purpose      : Compute the required number of servers (CSRs) to handle
                    given traffic load with given Average Handling Time and
                    no more than given Average Waiting Time
     Usage        : N=ErlcNsrvFromWait( traffErlang, AHT, AWT)
     Inputs       : traffErlang=Traffic load in Erlang
                    AHT=Average Handling Time
                    AWT=Average Waiting Time                    
     Handle Error : Return -1 if either time parameter is negative;                   
                    Return -2 if traffic in Erlang is negative
  */
         if AHT<0 or AWT<0 then do;
            *put "Both Average Holding Time and Average Waiting Time should be >0"; 
   Nsrv=-1; return(Nsrv);
   end;
   else if traffErlang<0 then do;
            *put "Traffic in Erlang must be positive real number";
   N0=-2; return(N0);
   end;
   else do;
           Nsrv=floor(traffErlang)+1;
           c=ErlcWait (Nsrv, traffErlang, AHT);
           do while(c>AWT);
                 Nsrv=Nsrv+1;
           c=ErlcWait (Nsrv, traffErlang, AHT);;
            end;
            return (Nsrv);
   end;
  endsub;

  function ErlcNwaiting(Nsrv, traffErlang);
  /*
        Entry        : ErlcNwaiting
     Purpose      : Compute the average number of calls waiting 
                    in queue, given number of servers and traffic
                    in Erlang
     Usage        : N=ErlcNwaiting(Nsrv, traffErlang)
     Inputs       : Nsrv=Number of Server
                    traffErlang=Traffic load in Erlang       
     Handle Error : Return -1 if traffic in Erlang is larger than number
                              of servers;                   
                    Return -2 if number of servers is negative
                    Return -3 if traffic in Erlang is negative
                    Return -4 if number of servers is not integer
  */
         if traffErlang>Nsrv then do;
            *put "Erlang amount should be smaller than Number of Servers";
   Q=-1; return(Q);
         end; 
   else if (Nsrv<0) then do;
            *put "Both Erlang amount and Number of Servers should be larger than 0";
   Q=-2; return(Q);
   end;
   else if traffErlang<0 then do;
            *put "Traffic in Erlang must be positive real number";
   Q=-3; return(Q);
   end;
   else if  abs(int(Nsrv)-Nsrv)>1E-8 then do;
            *put "Number of Servers should be integer";
   Q=-4; return(Q);
   end;  
         else do; 
            Q=traffErlang/(Nsrv-traffErlang)*ErlcFractionDelay(Nsrv, traffErlang);
      return(Q);
   end;
  endsub;

  function ErlcTrafFromFractionDelay(Nsrv, fractionDelay);
  /*
        Entry        : ErlcTrafFromFractionDelay
     Purpose      : Compute the max traffic load that N servers
                    can handle with queuing probability less or equal
                    to fractionDelay. 
     Usage        : A=ErlcTrafFromFractionDelay(Nsrv, fractionDelay)
     Inputs       : Nsrv=Number of Server
                    fractionDelay=Probability of waiting in queue
                                  when a call arrives                    
     Handle Error : Return -1 if probability of waiting in Queue is not
                              in range of (0, 1) 
                    Return -2 if number of servers is negative                   
                    Return -3 if number of servers is not integer
  */
  /*  NOTE:
         Using Bisection Method to solve for Erlang 
      given N server and Erlang C number. Bisection 
      is simple and robust. Speed is not of primiary
      concern here.
  */
            if fractionDelay>1 or fractionDelay<0  then do;
                *put "Number of servers must be positive integers";
    c=-1; return(c);
    end;
    else if Nsrv<=0 then do;  
                *put "Fraction of Delay must be real number between 0 and 1";
    c=-2; return(c);
    end;
    else if abs(int(Nsrv)-Nsrv)>1E-8 then do;;
       c=-3; return(c);
    end;
    else do;
                min=1E-8; max=Nsrv-1E-8;
             itermax=500; iter=1;
             do while (iter<=itermax);
                   c=(min+max)/2;
             _x=ErlcFractionDelay(Nsrv, c)-fractionDelay;
             if abs(_x)<1E-8 or abs(max-min)/2<1E-8 then do;                    
                iter=itermax+1;
             end;
             else do;
                      iter=iter+1;
       _tmp=ErlcFractionDelay(Nsrv, max)-fractionDelay;
                if sign(_x)=sign(_tmp) then max=c;
                      else min=c;
             end;
            end;      
      return (c);
   end;
  endsub;

  function ErlcWait (Nsrv, traffErlang, AHT);
  /*
        Entry        : ErlcWait
     Purpose      : Compute the Average Speed of Answer (ASA, also called
                    Average Waiting Time, AWT), given the number of servers, 
                    traffic in Erlang and average handling time. 
     Usage        : N=ErlcWait(Nsrv, traffErlang, AHT)
     Inputs       : Nsrv=Number of Server
                    traffErlang=Traffic load in Erlang       
                    AHT=Average Handling Time
     Handle Error : Return -1 if Average Handling Time is negative           
                    Return -2 if number of servers is negative
                    Return -3 if traffic in Erlang is negative
                    Return -4 if number of servers is not integer
                    Return -5 if traffic in Erlang is larger than the
                           number of servers
  */
      if AHT<0 then do;
            *put "Average Holding Time should be >0"; 
      ASA=-1; return(ASA);
   end;
   else if (Nsrv<0) then do;
            *put "Both Erlang amount and Number of Servers should be larger than 0";
   ASA=-2; return(ASA);
   end;
   else if traffErlang<0 then do;
            *put "Traffic in Erlang must be positive real number";
   ASA=-3; return(ASA);
   end;
   else if  abs(int(Nsrv)-Nsrv)>1E-8 then do;
            *put "Number of Servers should be integer";
   ASA=-4; return(ASA);
   end; 
   else if traffErlang>Nsrv then do;
      ASA=-5; return(ASA);
   end;  
   else do;
            ASA=ErlcFractionDelay(Nsrv, traffErlang)*AHT/(Nsrv-traffErlang);
      return(ASA);
   end;
  endsub;

  function ErlcK (Nsrv, traffErlang);
  /*
        Entry        : ErlcK
     Purpose      : Compute the average number of calls in queue, given 
                    number of servers and traffic in Erlang. It is the 
                    sum of average number of calls being served (traffErlang) 
                    and average number of calls waiting in queue.
     Usage        : N=ErlcK(Nsrv, traffErlang)
     Inputs       : Nsrv=Number of Server
                    traffErlang=Traffic load in Erlang       
     Handle Error : Return -1 if traffic in Erlang is larger than number
                              of servers;                   
                    Return -2 if number of servers is negative
                    Return -3 if traffic in Erlang is negative
                    Return -4 if number of servers is not integer
  */
        if traffErlang>Nsrv then do;
            *put "Erlang amount should be smaller than Number of Servers";
   K=-1; return(K);
         end; 
      else if (Nsrv<0) then do;
            *put "Both Erlang amount and Number of Servers should be larger than 0";
   K=-2; return(K);
   end;
   else if traffErlang<0 then do;
            *put "Traffic in Erlang must be positive real number";
   K=-3; return(K);
   end;
   else if  abs(int(Nsrv)-Nsrv)>1E-8 then do;
            *put "Number of Servers should be integer";
   K=-4; return(K);
   end;  
   else do;
            K = traffErlang*(1+ ErlcFractionDelay(Nsrv, traffErlang)/(Nsrv-traffErlang));
      return (K);
   end;
  endsub;

  function ErlcT( Nsrv, traffErlang, AHT);
  /*
        Entry        : ErlcT
     Purpose      : Compute the average time a call will spend in the system,
                    which is the sum of average handling time and average 
                       waiting time (ASA). Three relavent parameters are
                    number of servers, traffic in Erlang and average handling
                    time. 
     Usage        : N=ErlcT(Nsrv, traffErlang, AHT)
     Inputs       : Nsrv=Number of Server
                    traffErlang=Traffic load in Erlang   
                       AHT=Average Handling Time 
     Handle Error : Return -1 if average handling time is negative;                   
                    Return -2 if number of servers is negative
                    Return -3 if traffic in Erlang is negative
                    Return -4 if number of servers is not integer
                    Return -5 if traffic in Erlang is larger than number 
                           of servers
  */
      if AHT<0 then do;
            *put "Average Holding Time should be >0"; 
      T=-1; return(T);
   end;
   else if (Nsrv<0) then do;
            *put "Both Erlang amount and Number of Servers should be larger than 0";
   T=-2; return(T);
   end;
   else if traffErlang<0 then do;
            *put "Traffic in Erlang must be positive real number";
   T=-3; return(T);
   end;
   else if  abs(int(Nsrv)-Nsrv)>1E-8 then do;
            *put "Number of Servers should be integer";
   T=-4; return(T);
   end; 
   else if traffErlang>Nsrv then do;
      T=-5; return(T);
   end;
         else do; 
            T=AHT + AHT*ErlcFractionDelay(Nsrv, traffErlang)/((Nsrv-traffErlang));
      return(T);
   end;
  endsub;

run;
options cmplib=sasuser.funcs;


data _null_;
     lambda=667; AHT=150; AWT=20;
  A=lambda*150/3600;
  
  do j=1 to 10;
     N=floor(A)+j;
     Pc=ErlcFractionDelay(N, A);
  GoS=ErlcGoS(N, A, AHT, AWT); 
        Q=ErlcNwaiting(N, A);
        T=ErlcT(N, A, AHT); 
  K=ErlcK(N, A);
  rho=A/N;  
  put N=  @8  Pc= percent.4  @16 K= bestd8.1    @26  T= bestd8.1 
                @38 Q= bestd8.1    @50 Gos= bestd8.1  @62  rho= percent8.4;
  end;
run;

Stochastic Gradient Decending Logistic Regression in SAS

2012-05-24T05:19:00.001-07:00

Test the Stochastic Gradient Decending Logistic Regression in SAS. The logic and code follows the code piece of Ravi Varadhan, Ph.D from this discussion of R Help. The blog SAS Die Hard also has a post about SGD Logistic Regression in SAS.



filename foo url "http://www.biostat.jhsph.edu/~ririzarr/Teaching/754/lbw.dat" ;

data temp;
   infile foo length=len;
   input low age lwt race smoke ptl ht ui ftv bwt;
   put low age lwt race smoke ptl ht ui ftv bwt;
   if _n_>1;
run;
proc standard data=temp out=temp mean=0  std=1;
      var  age lwt smoke ht ui;
run;

proc contents data=temp  out=vars(keep=varnum  name type) noprint; run;

proc sql noprint;
      select name into :covars separated by " "
   from   vars
   where  lowcase(name) in ("age", "lwt", "smoke", "ht", "ui")
   ;
   select cats("b_", name) into :covars2 separated by " "
   from   vars
   where  lowcase(name) in ("age", "lwt", "smoke", "ht", "ui")
   ; 
     select count(*)+1 into :nparms
     from   vars
  where  lowcase(name) in ("age", "lwt", "smoke", "ht", "ui")
  ;
quit;
%put &covars2;

%lr_sgd(temp, beta, low, &covars, 
        alpha=0.01, decay=0.98, 
        criterion=0.00001, maxiter=1000);


options fullstimer;
proc logistic data=temp  outest=_beta  desc noprint;
      model low = age lwt smoke ht ui;
run;

Execution log shows:
******************************************************************************************
606

607 %lr_sgd(temp, beta, low, &covars, alpha=0.01, decay=0.98, criterion=0.00001, maxiter=1000);

MLOGIC(LR_SGD): Beginning execution.

MLOGIC(LR_SGD): Parameter DSN has value temp

MLOGIC(LR_SGD): Parameter OUTEST has value beta

MLOGIC(LR_SGD): Parameter RESPONSE has value low

MLOGIC(LR_SGD): Parameter COVARS has value age ht lwt smoke ui

MLOGIC(LR_SGD): Parameter ALPHA has value 0.01

MLOGIC(LR_SGD): Parameter DECAY has value 0.98

MLOGIC(LR_SGD): Parameter CRITERION has value 0.00001

MLOGIC(LR_SGD): Parameter MAXITER has value 1000

MPRINT(LR_SGD): options nosource nonotes;

MPRINT(LR_SGD): options nomlogic nomprint;

Iteration 1, time used 0.046, converge criteria is 0.355

Iteration 2, time used 0.032, converge criteria is 0.1886732527

Iteration 3, time used 0.14, converge criteria is 0.110486903

Iteration 4, time used 0.031, converge criteria is 0.0694367942

Iteration 5, time used 0.031, converge criteria is 0.0458381768

Iteration 6, time used 0.032, converge criteria is 0.0313858459

Iteration 7, time used 0.031, converge criteria is 0.0221141453

Iteration 8, time used 0.047, converge criteria is 0.0159488055

Iteration 9, time used 0.031, converge criteria is 0.0125402184

Iteration 10, time used 0.031, converge criteria is 0.0102548103

Iteration 11, time used 0.031, converge criteria is 0.008416673

Iteration 12, time used 0.047, converge criteria is 0.0069325467

Iteration 13, time used 0.032, converge criteria is 0.0057299824

Iteration 14, time used 0.031, converge criteria is 0.0047522323

Iteration 15, time used 0.031, converge criteria is 0.0039546116

Iteration 16, time used 0.031, converge criteria is 0.0033017935

Iteration 17, time used 0.047, converge criteria is 0.0027657536

Iteration 18, time used 0.031, converge criteria is 0.0023241894

Iteration 19, time used 0.032, converge criteria is 0.0019592976

Iteration 20, time used 0.046, converge criteria is 0.0016568231

Iteration 21, time used 0.032, converge criteria is 0.0014053171

Iteration 22, time used 0.031, converge criteria is 0.0011955575

Iteration 23, time used 0.031, converge criteria is 0.0010200936

Iteration 24, time used 0.047, converge criteria is 0.0008728877

Iteration 25, time used 0.031, converge criteria is 0.0007490333

Iteration 26, time used 0.032, converge criteria is 0.0006445313

Iteration 27, time used 0.031, converge criteria is 0.0005561129

Iteration 28, time used 0.047, converge criteria is 0.0004810987

Iteration 29, time used 0.031, converge criteria is 0.0004172864

Iteration 30, time used 0.031, converge criteria is 0.0003628606

Iteration 31, time used 0.031, converge criteria is 0.0003163211

Iteration 32, time used 0.032, converge criteria is 0.0002764247

Iteration 33, time used 0.047, converge criteria is 0.0002421384

Iteration 34, time used 0.031, converge criteria is 0.0002126018

Iteration 35, time used 0.047, converge criteria is 0.0001870963

Iteration 36, time used 0.046, converge criteria is 0.0001650204

Iteration 37, time used 0.032, converge criteria is 0.0001458692

Iteration 38, time used 0.031, converge criteria is 0.000129218

Iteration 39, time used 0.031, converge criteria is 0.0001147086

Iteration 40, time used 0.047, converge criteria is 0.0001020383

Iteration 41, time used 0.031, converge criteria is 0.0000909506

Iteration 42, time used 0.032, converge criteria is 0.0000812278

Iteration 43, time used 0.031, converge criteria is 0.0000726846

Iteration 44, time used 0.047, converge criteria is 0.0000651629

Iteration 45, time used 0.031, converge criteria is 0.0000585277

Iteration 46, time used 0.031, converge criteria is 0.0000526634

Iteration 47, time used 0.031, converge criteria is 0.0000474707

Iteration 48, time used 0.047, converge criteria is 0.0000428642

Iteration 49, time used 0.063, converge criteria is 0.0000387705

Iteration 50, time used 0.031, converge criteria is 0.0000351261

Iteration 51, time used 0.031, converge criteria is 0.000031876

Iteration 52, time used 0.032, converge criteria is 0.0000289726

Iteration 53, time used 0.031, converge criteria is 0.0000263748

Iteration 54, time used 0.031, converge criteria is 0.0000240465

Iteration 55, time used 0.031, converge criteria is 0.0000219565

Iteration 56, time used 0.031, converge criteria is 0.0000200774

Iteration 57, time used 0.047, converge criteria is 0.0000183854

Iteration 58, time used 0.172, converge criteria is 0.0000168595

Iteration 59, time used 0.031, converge criteria is 0.0000154815

Iteration 60, time used 0.032, converge criteria is 0.0000142352

Iteration 61, time used 0.031, converge criteria is 0.0000131064

Iteration 62, time used 0.031, converge criteria is 0.0000120826

Iteration 63, time used 0.031, converge criteria is 0.0000111528

Iteration 64, time used 0.032, converge criteria is 0.0000103072

Iteration 65, time used 0.031, converge criteria is 9.5372542E-6

Total Time is 2.56 sec.

Total Iteration is 65, convergence status is Converged.

At Final Iteration, max difference is 9.5372542E-6

MPRINT(LR_SGD): source;

MLOGIC(LR_SGD): Ending execution.

******************************************************************************************

The macro %LR_SGD.


/*
  SAS macro:
     Logistic Regression using Stochastic Gradient Descent.    
  Name: 
     %ls_sgd();
  Copyright (c) 2009, Liang Xie (Contact me @ xie1978 at gmail dot com)
  
  
  The SAS macro is a demonstration of an implementation of logistic 
  regression modelstrained by Stochastic Gradient Decent (SGD).This 
  program reads a training set specified as &dsn_in, trains a logistic 
  regression model, and outputs the estimated coefficients to &outest. 
  Example usage:

  %let inputdata=train_data;
  %let beta=coefficient;
  %let response=Event;
  %lr_sgd(&inputdata, &beta, &response, &covars, 
          alpha=0.008, decay=0.8, 
          criterion=0.00001, maxiter=1000);


  The following topics are not covered for simplicity:    
      - bias term    
      - regularization    
      - multiclass logistic regression (maximum entropy model)         
      - calibration of learning rate

  Distributed under GNU Affero General Public License version 3. This 
  program is free software: you can redistribute it and/or modify
  it under the terms of the GNU Affero General Public License as
  published by the Free Software Foundation, only version 3 of the
  License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU Affero General Public License for more details. 

*/

%macro logistic(dsn_in, outest, response, alpha=0.0005);
proc score data=&dsn_in  score=&outest  type=parms  out=score(keep=score);
     var intercept &covars;
run;

data _xtemp/view=_xtemp;
     merge &dsn score ;
     _w=&response - 1/(1+exp(-score));
  /*
  array x{*} intercept &covars;
  _w=&response - 1/(1+exp(-score));
  do i=1 to dim(x); x[i]=x[i]*_w; end;    
  */
run;

data _x&outest;
  array x{*} intercept &covars;
  array _x{*} b_intercept &covars2;
  retain b_intercept &covars2;
  retain logneg logpos 0;
  modify _x&outest;
  do i=1 to dim(x); x[i]=_x[i]; end;
  
  do until (eof);
        set _xtemp  end=eof;
     do  i=1 to dim(x);
         _x[i]=_x[i]+&alpha*x[i]*_w;
     end;
  end;
     replace;
run;

%mend;

%macro compare(dsn1, dsn2);
data _null_;
     merge &dsn1  &dsn2;
  array _x1{*} intercept &covars;
  array _x2{*} b_intercept &covars2;
  retain maxdiff 0;
  do i=1 to dim(_x1);      
     maxdiff=max(maxdiff, abs(_x1[i]-_x2[i]));  
  *put _x1[*]=;
  *put _x2[*]=;
  end;
  call symput('maxdiff', maxdiff);
run;
%mend;




%macro lr_sgd(dsn, outest, response, covars, 
              alpha=0.0005, decay=0.9, 
              criterion=0.00001, maxiter=1000);
options nosource nonotes;
options nomlogic nomprint;
%local i t0 t1 dt maxdiff status  stopiter;

%let t00=%sysfunc(datetime());

data &dsn;
     set &dsn;
  intercept=1; _w=1;
run;

data &outest;
     retain _TYPE_ "PARMS"  _NAME_ "SCORE";
     array x{*} intercept &covars;
  do i=1 to dim(x); 
     x[i]=0;
  end;
  drop i;
  output;
run;

data _x&outest;
     retain _TYPE_ "PARMS"  _NAME_ "SCORE";
     array bx{*} b_intercept &covars2;
  array x{*}  intercept &covars;
  set &outest;
  do j=1 to dim(x); bx[j]=x[j]; end;
  keep b_intercept &covars2 _TYPE_  _NAME_;
  drop j;
run;

sasfile _x&outest load;
%let stopiter=&maxiter;
%let status=Not Converged.;
%do i=1 %to &maxiter;
    %let t0=%sysfunc(datetime());

    %logistic(&dsn, &outest, &response, alpha=&alpha);
    %compare(&outest, _x&outest);    
    data &outest;
         retain _TYPE_ "PARMS"  _NAME_ "SCORE";
         array bx{*} b_intercept &covars2;
      array x{*}  intercept &covars;
      set _x&outest;
      do j=1 to dim(x); x[j]=bx[j]; end;
      keep intercept &covars _TYPE_  _NAME_;
      drop j;
    run;
 %let alpha=%sysevalf(&alpha * &decay);
 %let alpha=%sysfunc(max(0.00005, &alpha));
    %let t1=%sysfunc(datetime());
    %let dt=%sysfunc(round(&t1-&t0, 0.001));
    %put Iteration &i, time used &dt, converge criteria is &maxdiff; 
 %if %sysevalf(&maxdiff<&criterion) %then %do;
     %let stopiter=&i;
        %let i=%eval(&maxiter+1);
  %let status=Converged.;
 %end;
%end;
sasfile _x&outest close;
%let t11=%sysfunc(datetime());
%let dt=%sysfunc(round(&t11-&t00, 0.01));
%put Total Time is &dt sec.;
%put Total Iteration is &stopiter, convergence status is &status;
%put At Final Iteration, max difference is &maxdiff;
options mlogic mprint notes source;
%mend;