Introduction to Neural Networks

Supervised learning: neural networks, statistics, & machine learning

Initialize:

Off[General::spell1];

SetOptions[ListPlot,ImageSize->Small];
SetOptions[DensityPlot,ImageSize->Small];
SetOptions[ContourPlot,ImageSize->Small];

Overview

Previously...

Introduced the notions of loss and risk for inference.

Rationale for Bayesian from the point of view of supervised learning that minimizes empirical risk.

We did this by starting from a goal to learn a decision rule (or estimator) to minimize average loss over a set of training pairs (the empirical risk), and then showed that this is equivalent (in the limit) to making a Bayes optimal decision. This provides a tie between empirical methods such as support vector machines discussed in this lecture, and Bayesian decision theory.

Today

Set neural network supervised learning  in the context of various statistical/machine learning methods.

Rationale: The relationship of brain to behavior is complicated. Good to understand bottom-up, from neurons to behavior. But also good to understand top-down, from behavior to quantitative models with as few free parameters as possible. But with a view to plausible neural network functioning. 

In other words, always a good idea to try the simplest model first. 

Preview of “big data” and scientific python

Nearest neighbor demonstration: “random Xor”

Just because a problem isn’t linearly separable, doesn’t mean the algorithm has to be complicated. The nearest neighbor classifier is one of the simplest algorithms for supervised learning. If one has a prior rationale for a representation, and measure of distance (or similarity) that reflects how similar patterns cluster based on category, the nearest neighbor algorithm is a good thing to try.

k nearest neighbors or k-NN for short. Has no problem with Xor.

Define 𝒟 to be a mixture distribution of bivariate normals

𝒟=MixtureDistribution[{1,1,1,1},{MultinormalDistribution[{0,0.},{{1.,0},{0,1.}}],MultinormalDistribution[{0,10.},{{1.,0},{0,1.}}],MultinormalDistribution[{10,0.},{{1.,0},{0,1.}}],MultinormalDistribution[{10,10.},{{1.,0},{0,1.}}]}];

Here’s a sample:

RandomVariate[𝒟]

We can generate n = 100 samples:

n=100;

somedata=Table[RandomVariate[𝒟],{n}];

And then artificially assign labels, either 0 or 1 to each of the 100 points.

labels=Which[EuclideanDistance[#,{0,0}]<5,0, EuclideanDistance[#,{0,10}]<5,1,EuclideanDistance[#,{10,0}]<5,1,EuclideanDistance[#,{10,10}]<5,0]&  /@ somedata;

ListPlot[somedata]

We can color the plots with the labels, red for 1’s and blue for 0’s.

colorf=Blend[{{0,Blue},{1,Red}},#]&
pl=Graphics[MapThread[{colorf[#1],Point[#2]}&,{labels,somedata}],Axes->True,AspectRatio->1]

Dimensions[somedata]

Estimate nearest neighbor boundaries

OK, so we have our artificial data and a pretty plot. Let’s classify.

Run through lots of points inside {{0,10},{0,10}}. For each point find the k nearest neighbors, and count how many have labels of 1 (i.e. how many dots are “red”). Plot up the count proportions--i.e. plot  the proportion of 0s or 1s falling within a circle containing k=5 neighbors, for each (x,y) point. (“circle” is implicit in the default distance function used by Nearest[], which has EuclideanDistance as the default distance measure).

trainingdata= Thread[somedata->labels];

Clear[nf2];
nf2[x_,k_]:=Total[Nearest[trainingdata,x,k]]/k;

k=5;
Show[{DensityPlot[nf2[{x,y},k],{x,-2,12},{y,-2,12},PlotPoints->n,ColorFunction->"LightTerrain",Exclusions->None],pl}]

Make the decision: if the proportion is more than 1/2 decide 1, otherwise 0:

k=5;
ContourPlot[If[nf2[{x,y},k]>1/2,1,0],{x,0,10},{y,0,10},PlotPoints->60,ColorFunction->"LightTerrain",Exclusions->None]

k=1 and Voronoi diagrams

How to pick K? Note that small values of k produce high variance predictors--variations in future samples produce big differences in the decision boundary. If k includes all the training data, the prediction is the most frequently occurring label--i.e. behaves like it is basing its decision on the prior. A solution is to use cross-validation, where one divides the training pairs into a training set and a test or validation set.
And then compare the results of training with various values of k on the training set with respect to the ability to classify members of the test or validation set.

Nearest neighbor regression

“If something is similar to something else in one respect, it is likely to be similar in another respect.” Can be OK for smooth interpolation, but can have problem for extrapolation.

Return to linear separability

Discriminant functions

Before looking at traditional (Fisher) and more modern discriminants (linear and non-linear SVM), let’s build our geometric intuitions of what a simple threshold logic or perceptron unit does by viewing it from a more formal point of view. Perceptron learning is an example of nonparametric statistical learning, because it doesn't require knowledge of the underlying probability distributions generating the data (such distributions are characterized by a relatively small number of "parameters", such as the mean and variance of a Gaussian distribution). Of course, how well it does will depend on the generative structure of the data. Much of the material below is covered in Duda and Hart (1978).

Linear discriminant functions: Two category case

Let x be a feature vector. A discriminant function, g(x) divides input space into two category regions depending on whether g(x)>0 or g(x)<0. (We've switched notation, x=f). This is the hard threshold case, i.e. a step function.

This linear case then corresponds to the simple perceptron unit:

g(x)=w.x+(Subscript[w, 0])

where w is the weight vector and (Subscript[w, 0]) is the (scalar) threshold (which we’ve called bias, although this "bias" has nothing to do with statistical "bias"), with the decision: g(x)>0 or g(x)<0, corresponding to a step-function non-linearity.

Discriminant functions can be generalized, for example to quadratic decision surfaces:

g(x)=(Subscript[w, 0])+(∑_(i=1))(Subscript[w, i])(Subscript[x, i])+(∑_(i=1))(∑_(j=1))(Subscript[w, ij])(Subscript[x, i])(Subscript[x, j])

where x = {(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])...}. g(x)=0 defines a decision surface. For now we restrict ourselves to the linear case where the decision surface is  a hyperplane. 

Suppose (Subscript[x , 1])and (Subscript[x, 2]) are vectors, with endpoints sitting on the hyperplane (or in the simple case of the red line below), then their difference is a vector lying in the hyperplane

g((Subscript[x, 1])) = g((Subscript[x, 2])) =w.(Subscript[x, 1])+(Subscript[w, 0])=w.(Subscript[x, 2])+(Subscript[w, 0])
			w.((Subscript[x, 1])-(Subscript[x, 2]))=0

so because the dot product is zero, the weight vector w is normal to any vector lying in the hyperplane. Thus w determines how the plane is oriented. The normal vector w points into the region for which g(x)>0, and -w points into the region for which g(x)<0.

w.x+(Subscript[w, 0])=0, so w.x =-(Subscript[w, 0]).

If we project x onto the normalized weight vector x.w/|w|, we have the normal distance of the hyperplane from the origin equal to:

w.x/|w|=-(Subscript[w, 0])/|w|

Thus, the threshold determines the position of the hyperplane.

One can also show that the normal distance of an arbitrary vector x to the hyperplane is given by:

g(x)/|w|

Let xp be the vector to the point of projection on the discriminant. Then (x-xp) is the vector from x to xp.  If we project (x-xp) on to w/|w| (which is the unit vector perpendicular to the decision plane), then (x-xp).w/|w| is the distance we want. 
		 x.w - xp.w = x.w - g(xp) + (Subscript[w, 0]) = x.w + (Subscript[w, 0]) = g(x) 
So the distance is g(x)/|w|.

To summarize: 
	1) discriminant function divides the input space by a hyperplane decision surface; 
	2) The orientation of the surface is determined by the weight vector w; 
	3) the location is determined by the threshold (Subscript[w, 0]); 
	4) the discriminant function gives a measure of how far an input vector is from the hyperplane.
	
The figure summarizes the basic properties of the linear discriminant.

Clear[w1,w2,w,w0,x2,x1];
Manipulate[
(* x0={2,1};*)
w={w1,w2};wn=w/Norm[w];
g[x_]:={w1,w2}.x+w0;
gg=Plot[Tooltip[x2 /. Solve[{w1,w2}.{x1,x2}+w0==0,x2],"discriminant"],{x1,-1,3}];
ggg=Graphics[g[Dynamic[MousePosition["Graphics"]]]];
Show[{gg,Graphics[Inset["g[x]=",{1.6,2}]],Graphics[Inset[ToString[g[x0]],{2,2}]],Graphics[{Tooltip[Arrow[{{0,0},w}],"w"],Tooltip[Arrow[{{0,0},(-w0/Norm[w])*wn}],"-(Subscript[w, 0])/|w|"],Tooltip[{Arrow[{{0,0},x0}]},"x"],Tooltip[{Arrow[{x0,x0-wn*g[x0]/Norm[w]}]},"g(x)/|w|"]}]},PlotRange->{{-1,3},{-1,3}},AxesOrigin->{0,0},Axes->True,AspectRatio->1,ImageSize->Medium],{{w0,-2.5},-6,3},{{w1,1},0,3},{{w2,2},0,3},{{x0,{2,1}},Locator}]

Multiple classes

Suppose there are c classes. There are a number of ways to define multiple class discriminant rules.  One way that avoids undefined regions is:

(Subscript[g, i])(x)=(Subscript[w, i]).x+(Subscript[w, i0]), i=1,...,c

Assign x to the ith class if: (Subscript[g, i])(x)>(Subscript[g, j])(x) for all j≠i.

(Adapted from Duda & Hart, 1973)

It can be shown that this classifier partitions the input space into simply connected convex regions. This means that if you connect any two feature vectors belonging to the same class by a line, all points on the line are in the same class. Thus this linear classifier won't be able to handle problems for which there are disconnected clusters of features that all belong to the same class. Also, from a probabilistic perspective, if the underlying generative probability model for a given class has multiple peaks, this linear classifier won't do a good job either.

Fisher linear discriminant

Basic idea: find a projection that minimizes the within-class variance, while maximizing the between-class variance.

Task-dependent Dimensionality Reduction

Motivation

Later, when we consider the problem of dimensionality reduction, we will look at unsupervised methods, in particular principal components, PCA. But here the idea will be to find hyperplanes onto which we can project our input data, and from there divide up the hyperplane into decision regions. The idea is that the original input space may be impractically huge, but if we can find a subspace (hyperplane) that preserves the distinctions between categories as well as possible, we can make our decisions in smaller space. We will derive the Fisher linear "discriminant".

This is closely related to the psychology idea of finding "distinctive" features. E.g. consider bird identification. If I want to discriminate cardinals from other birds in my backyard, I can make use of the fact that (males) cardinals may be the only birds that are red. So even though the image of a bird can have lots of dimensions, if I project the image on to the "red" axis, I can do fairly well with just one number. How about male vs. female human faces?

A conceptually similar idea was used to learn diagnostic features in complex images (cf. Ullman et al., Hegde et al.).

We’ll look at this more later when we learn about how to discover feature hierarchies.

Fisher's linear discriminant: building intuitions

Generative model: two nearby gaussian classes

Define two bivariate base distributions

(ar = {{1, 0.99}, {0.99, 1}};
ndista = MultinormalDistribution[{0, -1}, ar];)
(br = {{1, .9}, {.9, 2}};
ndistb = MultinormalDistribution[{0, 1}, br];)

Find the expression for the probability distribution function of ndista

pdf = PDF[ndista, {x1, x2}]

Use Mean[ ] and CovarianceMatrix[ndista] to verify the population mean and the covariance matrix of ndistb

Mean[ndistb]

Covariance[ndista]

pdfa=PDF[ndista,{x1,x2}];
(*ContourPlot[pdfa,{x1,-3,3},{x2,-3,3},PlotPoints->64,PlotRange->All]*)

nsamples=500;
a=Table[RandomVariate[ndista],{nsamples}];
ga=ListPlot[a,PlotRange->{{-4,4},{-4,4}},AspectRatio->1,PlotStyle->Hue[0.2`]];
b=Table[RandomVariate[ndistb],{nsamples}];
gb=ListPlot[b,PlotRange->{{-4,6},{-4,4}},AspectRatio->1,PlotStyle->Hue[0.6`]];
Show[ga,gb]

Use Mean[ ] to find the sample mean of  b. Whats is the sample covariance of b?

Mean[b]

Covariance[b]

Try out different projections of the data by varying the slope (m) of the discriminant line

The basic idea is fairly straight forward, and if you are familiar with signal detection theory and the definition of d’, it should look familiar. We want to find a projection for which the squared difference between the means of the projections relative to the variance

(((((Subscript[μ, b])-(Subscript[μ, 1])))^2)/((Subscript[σ, a]^2)+(Subscript[σ, b]^2)))

is biggest.

Manipulate[
wnvec={1,m}/Sqrt[1+m^2];
aproj=(#1 wnvec&)/@(a.wnvec);
gaproj=ListPlot[aproj,AspectRatio->1,PlotStyle->Hue[0.3]];
bproj=(#1 wnvec&)/@(b.wnvec);
gbproj=ListPlot[bproj,AspectRatio->1,PlotStyle->Hue[0.7]];
Show[ga,gb,gbproj,gaproj,AspectRatio->Automatic],
{{m,-.66},-2,2}]

Calculate the "signal-to-noise" ratio along the projection line.  Do this by taking the difference between the means divided by the square root of the product of the standard deviations along the line.

Theory and program demo for simple 2-class case

Here’s the derivation of the Fisher discriminant for the 2-class case (see Duda and Hart for general case).

A measure of the separation between the projections is the difference between the means:

|w.((Subscript[m, a])-(Subscript[m, b]))|
		and 
		(Subscript[m, a])=(Underoverscript[(1/N)∑, (i=1), N])x, summed over the N  x's from class a
		(Subscript[m, b])=(1/M)∑x, summed over the M  x's from class b

where w (wnvec) is the unknown unit vector along the discriminant line .

In our case above, the vector difference between the means is:

Mean[a]-Mean[b]

and the difference between the means projected onto a discriminant line is:

wnvec={1,2/3}/Sqrt[1+2/3^2];

wnvec.(Mean[a]-Mean[b])

To improve separation, we can't just scale w, because the noise scales too.

We'd like the difference between the means to be large relative to the variation for each class. We can define a measure of the scatter for the projected samples in say class a (a==1), by:

(∑_(y∈class a))(((y-(Subscript[(OverTilde[m]), a])))^2)

where (Subscript[(OverTilde[m]), a])is the sample mean of the points from class a projected onto discriminant line and y=w.x 
Or in terms of the Mathematica example:

Apply[Plus,aproj-wnvec.Mean[a]];

The total scatter S is defined by the sum of the scatters for both classes (a and b).

S=(∑_(y∈class a))(((y-(Subscript[(OverTilde[m]), a])))^2)+(∑_(y∈class b))(((y-(Subscript[(OverTilde[m]), b])))^2)

If we divide the above number by the total number of points, we have an estimate of the variance of the combined data along the projected axis.

We now have the basic ingredients behind the intuition for the Fisher linear discriminant. We'd like to find that w for which J:

J(w)=((|(Subscript[(OverTilde[m]), a])-(Subscript[(OverTilde[m]), b])((|)^2))/S)=((|w.((Subscript[m, a])-(Subscript[m, b]))((|)^2))/S)

is biggest. We want to maximize the difference between the projected class means, while minimizing the dispersion of the data on the projected line.

One can show that S = (w^T).(Subscript[S, W]).w, where

(Subscript[S, W]) is measure of within-class variation called the within-class scatter matrix:

(Subscript[S, W])=(∑_(i=1)^2)(∑_(x∈Class i))(x-(Subscript[m, i]))(((x-(Subscript[m, i])))^T)

For the numerator, a measure of between class variation is the between-class scatter matrix:

(Subscript[S, B])=((Subscript[m, 1])-(Subscript[m, 2])).((((Subscript[m, 1])-(Subscript[m, 2])))^T)

and the difference between the projected means can be show to be:

|(Subscript[(OverTilde[m]), a])-(Subscript[(OverTilde[m]), b])((|)^2)=(w^T).(Subscript[S, B]).w

Find w (corresponding to slope) to maximize the criterion function

J(w)=(((w^T).(Subscript[S, B]).w)/((w^T).(Subscript[S, W]).w))

The answer is given by:

w=((Subscript[S, W])^(-1)).((Subscript[m, a])-(Subscript[m, b]))

Demo: Program to finding Fisher's linear discriminant

normalize[x_]:=x/Sqrt[x.x];

ma=Mean[a];
mb=Mean[b];

Sa=Sum[Outer[Times,a[[i]]-ma,a[[i]]-ma],{i,1,nsamples}];
Sb=Sum[Outer[Times,b[[i]]-mb,b[[i]]-mb],{i,1,nsamples}];
Sw=Sa+Sb;
wldf = normalize[Inverse[Sw].(ma-mb)]

aproj=(#1 wldf&)/@(a.wldf);
gaproj=ListPlot[aproj,AspectRatio->1,PlotStyle->Hue[0.3`]];
bproj=(#1 wldf&)/@(b.wldf);
gbproj=ListPlot[bproj,AspectRatio->1,PlotStyle->Hue[0.7`]];
Show[ga,gb,gbproj,gaproj]

We started off with a 2-dimensional input problem and turned it into a 1-D problem. For the n-dimensional case, see Duda and Hart.

Compare with the principal component axes

Later we will look at principal components analysis as an unsupervised method to reduce dimensionality. But it is useful at this point to get an intuition to compare the two different ways to reduce dimensionality. The first component of PCA finds the dimension (axis) for which if you project the data on to it, the projections have the biggest variance.

c=Join[a,b];
ListPlot[c,PlotRange->{{-4,4},{-4,4}}]

g1=ListPlot[c,PlotRange->{{-4,4},{-4,4}},AspectRatio->1];

auto=Covariance[c];
eigvalues=Eigenvalues[auto];
eigauto=Eigenvectors[auto];

gPCA = Plot[{eigauto[[1,2]]/eigauto[[1,1]] x,
	eigauto[[2,2]]/eigauto[[2,1]] x},
		{x,-4,4}, AspectRatio->1,
		PlotStyle->{RGBColor[.2,0,1]}];

Show[g1,gPCA]

If the main information for discrimination is in the variance, rather than the means, PCA can be better.

ndistc = MultinormalDistribution[{0, 0},{{.05,0},{0,.05}}];
c=Table[RandomVariate[ndistc],{nsamples/2}];
gc=ListPlot[c,PlotRange->{{-3,3},{-3,3}},AspectRatio->1,PlotStyle->Hue[0.9]];
Show[ga,gc]

Support vector machines (SVMs) and kernel methods

Initialize

Read in SVM package:

Make sure the SVM package is downloaded in the default directory

SetDirectory[NotebookDirectory[]]

<<MathSVMv7`

Sidenote: Some other useful directory functions:

FileNames[];
Directory[];

We assume a simple perceptron TLU to classify vector data x into one of two classes depending on the sign of g(x):

decision(x)=sign(w.x+b).

Given g(x)=w.x+b, recall that g(x)/||w|| is the distance of a data point x from a plane defined by g(x) = 0. In support vector machines, the goal is to find the separating plane (i.e. find w and b) that is as far as possible from any of the data points. The intuition is that this will minimize the probability of making wrong classifications when given new data at some future point.

Let (Subscript[Π, w,b]) represent the points on a potential decision line, and d((Subscript[Π, w,b]),(Subscript[x, i]))is the distance between data point  (Subscript[x, i]) and the line. Formally, we want to solve:

(max_((w,b)))((min_i)d((Subscript[Π, w,b]),(Subscript[x, i]))),

where d((Subscript[Π, w,b]),(Subscript[x, i]))= g((Subscript[x, i]))/||w||, i.e.

d((Subscript[Π, w,b]),(Subscript[x, i]))= g((Subscript[x, i]))/||w||=|w.(Subscript[x, i])+b|/||w||

Two-class data (black and grey dots), their optimal separating hyperplane (continuous line), and support vectors (circled in blue). This is an example output of the SVMPlot function in MathSVM (Nilsson et al., 2006). The width of the “corridor” defined by the two dotted lines connecting the support vectors is the called the margin of the optimal separating hyperplane. If the data are not linearly separable, the margin can be “softened” to allow some points to fall inside the corridor. The  soft-margin classifier uses a regularizing term (similar to a prior on the parameters) to allow for a trade-off between classification error and tolerance.

The Primal Problem

It can be shown that the optimal separating hyperplane solving (1) can be found as the solution to the equivalent optimization problem

(min_(w,b)) (1/2)||w((||)^2)
			subject to (Subscript[y, i])((w^T)(Subscript[x, i])+b)≥1,

Where (Subscript[y, i]) is the class label, either 1 or -1, for the ith point. Typically, equality will hold for a relatively small number of the data vectors. These data are termed support vectors. The solution (w,b) depends only on these specific points, and in effect contain all the information for the decision rule. The "dual problem".

A simple linear SVM Example

Here's a demo of a simple SVM problem. It uses the add on package MathSVMv7 written by Nilsson et al. (2006)

We use the simple SVM formulation  provided in MathSVM by the SeparableSVM function.

τ=0.01;
α=SeparableSVM[X,y,τ]

In the output figure below, the solid line marks the optimal hyperplane, and dotted lines mark the width of the corridor that joins support vectors (highlighted in blue).

SVMPlot[α,X,y]

A Nonlinear Example: Using Kernels

What if the data are not linearly separable? The essential idea is to map the data (through some non-linear mapping, e.g. polynomial) to a higher-dimensional "feature" space to find the optimal hyperplane separating the data. The dot product gets replaced by a non-linear kernel function. For example, the polynomial kernel is given by:

If d = 1, we have the standard dot product, but for d = 2, 3, etc.. we have polynomial functions of the elements of the vectors x. See Nilsson et al (2009), and paper by Jäkel (2009) for more information on kernels.

A demo for an application for nonlinear classification.

We’ll use second-degree kernel, e.g.:

PolynomialKernel[{xx1,xx2,xx3},{1,-1,-1},2]

Some synthetic data which is not linearly separable.

len=50;
X=Join[
RandomReal[NormalDistribution[0,0.03],{len/2,2}],
Table[
{RandomReal[NormalDistribution[i/len-1/4,0.01]],Random[NormalDistribution[(((2i/len-1/2))^2)-1/6,0.01]]},
{i,len/2}]];
y=Join[Table[1,{len/2}],Table[-1,{len/2}]];
SVMDataPlot[X,y,PlotRange->All]

We use the KernelFunction to specify the kernel type and run SeparableSVM[].

τ=0.01;
pk=PolynomialKernel[#1,#2,2]&;
α=SeparableSVM[X,y,τ,KernelFunction->pk]

When visualizing the results, SVMPlot can use the kernel functions to draw any nonlinear decision curves.

SVMPlot[α,X,y,KernelFunction->pk]

More information

The wiki entry for SVMs has a good introduction and overview.

To go to the source see Vapnik (1995)

And for applications of kernel methods more generally, to cognitive and neuroscience see reviews by Jäkel et al. (2006; 2009).

Logistic regression, the generalized linear model, and the generic neuron model

This section is not about learning per se, but about how to motivate the class of functions one might like to fit. 

The single (weight) layer perceptron is a special case of logistic regression, which in turn is a special case of a generalized linear model. Logistic regression has been studied by statisticians since the 1950s. And the generalized linear model is a basic tool in theoretical modeling of spiking neurons, in particular the linear-Nonlinear-Poisson model (LNP) is a special case, the first part of which should look familiar. In the figure below, the set of weights specifying the receptive field of a neuron are given by k:

(Figure due to Jonathan Pillow, lecture in Methods in Computational Neuroscience, 2011.)

Here’s a general mathematical motivation for modeling the output as a probability of binary events. 

The starting point is that we want to model a binary response Y (think action potential on or off) by the conditional probability:

	P(Y=1 | X=x ), 
	
where X are the input features. Let's write

	P(Y=1 | X=x ) = p(x;w),
	
where p(x;w) is a likelihood function of x for some parameters w (“k” in the above figure). How to model this function p? 

As we've done before, let's assume the simplest: try a linear function. But over what function of p? 

p(x) = w.x + b has problems because the left side is bounded between 0 and 1, and the right side is unbounded. 

Log p(x) = w.x + b may seem nice because adding input features multiplies probabilities, but is bounded on only one side. A solution is to model log odds (of a spike vs. no spike) as a linear function:

	Log ((p(x))/(1-p(x)))= w.x + b

and after solving for p(x),

	p(x) =(( 1)/(1+(e^(-(w.x + b)))))

Looking familiar? This says that the probability of Y=1 is determined by taking a linear weighted sum of the inputs plus b (think bias) and deciding “yes” if this is bigger than zero. This is a linear classifier. This is also the probabilistic update rule we used for the Boltzmann machine. 

To make things look even more familiar let’s calculate the expected value of  Y=1 (of a neuron firing). First note that we can be a bit clever and write the probability of Y= (Subscript[y, i]), whose values are 1 or 0 in a concise summary expression as:

		p((Subscript[y, i]) | x) =((  (e^((Subscript[y, i]) (w.x + b))))/(1+(e^((w.x + b)))))
		
Then the probability of firing, i.e. (Subscript[y, i]) = 1 is:

		p((Subscript[y, i]) = 1| x) =  (( 1)/(1+(e^(-(w.x + b))))) = σ(w.x + b), where as we’ve seen before, σ is sigmoidal.

By the definition of expectation, the average rate is:

		(∑_((Subscript[y, i])=0)^1)(Subscript[y, i]) p((Subscript[y, i]) | x) = 0 × p((Subscript[y, i]) =0| x) + 1×p((Subscript[y, i]) = 1| x) = σ(w.x + b)

This is stage 1 and 2 of our original generic neuron model.

The logistic regression model can be extended to multiple class decisions.

And as introduced above, these models in turn can be viewed as special cases of generalized linear models. 

Note: The general linear model is different! It really is linear, and can be viewed as a special case. The generalized linear model is in general not. For applications of the generalized linear model to computational neuroscience,  see Pillow (2007).

Naive Bayes and other methods

Naive Bayes. Assume we have a class variable C that takes on discrete values corresponding to labels. And assume a set of features {(Subscript[F, i])}. Naive Bayes assumes that the probabilities of features are all conditionally independent of the class they came from:

(You can view this as an extension to the cue integration problem studied earlier.) To do classification, we need  a decision rule, such as MAP, which chooses the C that maximizes the posterior.

(Underscript[argmax, C]) p(C | (Subscript[F, 1]),...,(Subscript[F, n]))

http://en.wikipedia.org/wiki/Naive_Bayes_classifier

And others, AdaBoost, http://en.wikipedia.org/wiki/AdaBoost,  and Random Forest, This Random Forest link also has a concise description of k-NN. http://en.wikipedia.org/wiki/Random_forest

Mathematica’s Classify[ ] function demo

Mathematica tries to make it easy for non-experts to use its machine learning functions. This is convenient for learning concepts, and for many problems. The downside is that it can be difficult to find out exactly what the routines are doing.

Here we adapt code from one of Mathematica’s Classify[] examples and apply it to the noisy Xor problem. Define  clusters sampled from normal distributions:

Set up training pairs

std = 0.05;
sampledata[center_]:=RandomVariate[MultinormalDistribution[center,std*IdentityMatrix[2]],30];

Generate clusters, assign labels to each and plot;

clustercenters = {{-1,-1},{-1,1},{1,-1},{1,1}}
clusters=sampledata/@clustercenters;
colors={Blue, Yellow,Yellow, Blue};
pl2=ListPlot[clusters, PlotStyle->Darker@colors]

Dimensions[clusters]

The finagling in the next line shows how to assign all of the 2D points to their corresponding colors

trainingdata0 = Flatten[Thread[Transpose[clusters][[#]]->colors]& /@Range[Dimensions[clusters][[2]]],1];

Short[trainingdata0]

Set up the classifier, and classify

Now run the training pair associations through the Classifier, to automatically find the method, and do classifications on test data.

c30=Classify[trainingdata0]

What method did it pick? With the default, setting Mathematica decides.

How to calculate probabilities and make decisions:

c3[{-2,3},"Probabilities"]

ClassifierInformation[c3, "MethodDescription"];
ClassifierInformation[c3, #]& /@ {"NeighborsNumber", "DistanceFunction"};
ClassifierInformation[c3];
ClassifierInformation[c3,"Properties"];

ctemp=Classify[trainingdata0, Method->{"NearestNeighbors","NeighborsNumber"->4}]

Visualize decision probabilities

plotprob[method_]:=Module[{},
c3=Classify[Flatten[Thread[Transpose[clusters][[#]]->colors]& /@Range[30],1],Method->method];
DensityPlot[Normal@c3[{x,y},"Probabilities"][[1]],{x,-5,5},{y,-5,5},Exclusions->None,ColorFunction->"LightTerrain"]
];

plotprob/@{"LogisticRegression","NearestNeighbors","RandomForest","NeuralNetwork","SupportVectorMachine","NaiveBayes"}

plotprobabilities[method_]:=Module[{c, proba},
c3=Classify[Flatten[Thread[Transpose[clusters][[#]]->colors]& /@Range[30],1],Method->method];proba=Table[Append[{x,y},#]&/@Lookup[c3[{x,y},"Probabilities"],colors],{x,-5,5,.5},{y,-5,5,.5}];proba=Flatten[#,1]&/@Transpose[proba,{3,2,1,4}];
Show[
ListPlot3D[proba,PlotRange->{0,1},PlotLabel->method],
ListPointPlot3D[Map[Append[#,1]&,clusters,{2}],
PlotStyle->({Opacity[.9],#}&/@colors)], ImageSize-> 150
]
];

Visualize decision probabilities using 3D plots

plotprobabilities/@{"LogisticRegression","NearestNeighbors","RandomForest","NeuralNetwork","SupportVectorMachine","NaiveBayes"}

Next time

Doing inference by sampling

Transitioning to Python

Why?

Pros:
	Free

	Rapid growth in packages for scientific computation, including spiking neural networks (http://briansimulator.org), machine learning, image processing, computer vision, and psychophysics (PsychoPy).
	
	Lots of packages, examples, tutorials. More specialized packages than Mathematica. For example, deep convolutional networks: http://deeplearning.net/tutorial/lenet.html.

Cons: 
	Lots and lots of packages, examples, tutorials. Packages are a set of moving targets that are may live or fade away and die. Not as many specialized packages as matlab. Installation can be a pain.

Scientific python: numpy, scipy, ...

Raw python is too basic. We need add on modules. So in order of importance, we will always start by loading numpy for creating and manipulating numerical arrays, scipy for scientific functions (e.g. optimization, ..), and matplotlib for graphics:

We’ll also use IPython which brings the notebook functionality of Mathematica to python.

In addition the modules, pandas,  scikit-learn,  and scikit-image, provide a wide range of machine learning and an image processing functions. Lots of others, such as neurolab.

In this course, the main purpose for python for is for the lecture material on MCMC sampling where we’ll be using pymc. 

In order to get started, we’ll first go to a terminal window and see if python is there. Execute a demo program. Then
we’ll open a new IPython notebook and bring the python program into the notebook.

References

Duda, R. O., Hart, P. E., & Stork, D. G. (2001). Pattern classification (2nd ed.). New York: Wiley.
(Amazon.com)
Hegdé, J., Bart, E., & Kersten, D. (2008). Fragment-based learning of visual object categories. Current Biology, 18(8), 597–601. doi:10.1016/j.cub.2008.03.058
Jakel, F., Schölkopf, B., & Wichmann, F. A. (2013). Generalization and similarity in exemplar models of categorization: Insights from machine learning. Psychonomic Bulletin & Review, 15(2), 256–271. doi:10.3758/PBR.15.2.256
Vapnik, V. N. (1995). The nature of statistical learning. New York: Springer-Verlag.
http://neuron.eng.wayne.edu/software.html
Ullman, S., Vidal-Naquet, M., & Sali, E. (2002). Visual features of intermediate complexity and their use in classification. Nature Neuroscience, 5(7), 682–687. doi:10.1038/nn870
Yuille, A., Coughlan J., Kersten D. (1998) (pdf)