% This script will illustrate how Logistic regression works.  Logistic regression is a method for classification: given some data points, we would like to separate them into various distinct categories.  Logistic regression is most often used when the categories are binary: yes/no, dog/cat, etc.  To use logistic regression, we must possess a dataset of some already classified points in order to train our regression (estimate parameters).  We can then use this trained classification method to try to classify some test data where the true category assignments are unknown (or treated as such, for verification purposes).
% First, let's just make a simple dataset to use to train our method.  Here, our two categories will just be called negative (n) and positive (p).  The data is constructed such that each object category has a set of points associated with it -- \{x_n\}, \{x_p\} -- and for this example these sets are generated randomly from two-dimensional multivartiate normal distributions, each with the same (diagonal) covariance matrix, but with different means for the negative and positive groups.

nn=1000; % number of negative and positive points
np=1000;
sd=1; % a value for the standard deviation of the two normals
xn=mvnrnd([-1,0],sd*eye(2),nn); % generate data, using two different means to make the groups somewhat distinct
xp=mvnrnd([1,1],sd*eye(2),np);
scatter(xn(:,1),xn(:,2)) % visualize the data
hold on;
scatter(xp(:,1),xp(:,2))
hold off;

% Looking at the plot to the right, the negative points are blue and the positive points are red.  As we can see, there is a relatively distinct separation of the two groups in terms of their respective datasets, but there are regions where these two groups overlap a fair bit.  These overlapping regions are what lead to potential challenges in classifying.  After working through the script with the default value for "sd", you might go back and make it larger (more overlap) or smaller (less overlap) to see how that affects performance of the logistic regression below.
%
%
% Now we will train our logistic classifier.  The basic assumption underlying logistic regression is that each datapoint i in our full dataset (which includes the positive and negative data) has a probability P(x_i) of belonging to the positive group (the choice of positive group is arbitrary here).  Moreover, the specific functional form for that probability is
% P(x_i|\theta)=\frac{e^{x_i\theta}}{1+e^{x_i\theta}}
% A few words about this formula.  First, x_i here is taken to be a row vector with m+1 entries; the first m entries are just whatever the data provides (so, m=2 for our simulated data here), while the last entry is just an appended "1" to allow for a constant offset for all points (analogous to the y-intercept in linear regression).  Accordingly, \theta is an m+1 dimensional column vector of coefficients, and our task will be to estimate good values for these coefficients given our data (again, this is analogous to linear regression).  The specific functional form of this probability is called the logistic (or logit) model, and has the nice property of being continuous, monotonically increasing, and bounded between 0 and 1 as its argument x_i\theta ranges from negative infinity to infinity.
% It is of use to note that when x_i\theta=0, the point has equal odds of being positive or negative.  Consequently, if x_i\theta>0 the point is more likely to be positive than negative, and if x_i\theta<0 the point is more likely to be negative than positive.  So, if we were to use these probabilities to make educated guesses about which points are positive versus negative, there is a clear discriminator: x_i\theta=0.  This linear structure is referred to as the decision boundary, as one can see that all points on a specific side of this structure (the side where x_i\theta>0) should be classified as positive, while those on the other side should be classified as negative.  One could actually be a bit more general here, and choose to classify a point as positive only if P(x_i|\theta)>\tau, where \tau is some chosen threshold between 0 and 1.  In this case there is still a clear linear structure that divides the data, but it is now x_i\theta=\textrm{logit}(\tau).
% To proceed with logistic regression, assume that the above logistic model is correct for our data.  Then the categories for our data can be interpreted as the output of a series of independent random events, each with probability P(x_i) to yield the value "positive" (and the complementary probability to yield "negative").  Suppose we use y_i to represent the outcome of this random process for datapoint i, with y_i=1 meaning positive and y_i=0 meaning negative.  Then the probabily distribution for y_i is bernoulli:
% p(y_i|x_i,\theta)=P(x_i|\theta)^{y_i}~\left[1-P(x_i|\theta)\right]^{1-y_i}
% Finally, if we use \hat{y}_i to denote the observed category for datapoint i , we can write down the likelihood of our full set of (independent) observations as
% L=\prod_i p(\hat{y}_i|x_i,\theta)
% Then the goal of logistic regression is to find the vector \theta that maximizes this likelihood.  This is generally done by minimizing the negative of the log-likelihood, \mathcal{L}.  By substituting the various defintions above and performing some algebra, one can find that this amounts to minimizing
% \mathcal{L}=\sum_i \ln\left[1+e^{x_i\theta}\right]-y_i~x_i\theta

data=[zeros(nn,4);ones(np,4)]; % represent each datpoint by [x 1 yhat]
data(:,1:2)=[xn;xp];
data(:,3)=ones(nn+np,1);
loglike=@(c) sum(log(1+exp(data(:,1:3)*c)))-data(:,4)'*(data(:,1:3)*c); % define the negative log-likelihood
copt=fmincon(loglike,[0;0;0],[],[]); % have matlab minimize it to find optimal parameters copt
x1max=max([xn(:,1);xp(:,1)]); % for visualization, then show the decision boundary
x1min=min([xn(:,1);xp(:,1)]);
tau=.5;
q=log(tau/(1-tau));
hold on;
fplot(@(x) (-(-q+copt(3)+copt(1)*x))/copt(2),[x1min,x1max])
hold off;

% We see in our image to the right that the decision boundary divides the positives and negatives as best one can given the region of overlap in this data.  If you go back and alter the value of "sd", you might find that some "strange" things happen if the data is 100% non-overlapping.  This is because in that case, there are a lot of different lines that would each peform perfectly on the data, so the minimization might have a hard time settling on a specific one.
%
% The plot above gives a nice way to visualize the performance of logistic regression on this training data.  Another way is to just look through each datapoint and determine how the trained regression would classify that point, then compare to its true classification; this is called a confusion chart.  Note that this is not a completely fair assessment of the method, since we are now using the same data to train it and evaluate it, but we'll do something a little more appropriate later.
% Another way to assess how logistic regression performs on this data is to construct an ROC curve.  The ROC curve is a parametric plot of the true positive rate (what fraction of positive points are classified as positive) vs the false positive rate (what fraction of negative points are classified as positive) while varying the prediction threshold \tau mentioned above.  All ROC curves must certainly contain two points: (0,0), which would occur if no point is ever classified as positive and happens if \tau=1; and (1,1), which would occur if all points are classified as positive and happens if \tau=0.  Every time \tau decreases, by definition there cannot be any fewer points classified as positive, so the false positive and true positive rates can only increase or stay the same.  Hence, the ROC curve non-decreasing between the points (0,0) and (1,1).  Further, if the method is able to perfectly classify the data for some \tau value, then the ROC curve will contain the point (0,1): no false positives and all positives classified correctly (so no false negatives either) at whatever \tau value this occurs at.  Given the non-decreasing nature of the curve, if the data is perfectly classifiable by the method, then the curve is trivial: it is a vertical line connecting (0,0) and (0,1), then a horizontal line connecting (0,1) to (1,1).  For all of these reasons, the ROC curve is often summarized by simply computing the area underneath it, called the AUC (area under curve), which will be 1 for a dataset that is perfectly classifiable by this method, and should be around 0.5 if the data were classified by a totally random classifier (that just assigned each point positive or negative randomly based on the overall fracion of each occurring in the dataset).  For a given dataset, the better the classifier, the closer the AUC will get to 1 and away from 0.5.

q=data(:,1:3)*copt; % find P for each point, then classify off of it
prob=exp(q)./(1+exp(q));
pred=round(prob);
confusionchart(confusionmat(data(:,4),pred,'Order',[0 1]),char('n','p')) % plot of predictions vs truth
acc_train=sum(pred==data(:,4))/(nn+np) % compute overall accuracy
taus=linspace(0,1,1000); % compute ROC curve
rocmat=zeros(1000,2);
for i=1:1000
    pred=(prob>=taus(i));
    tpr=sum((pred==1)&(data(:,4)==1))/np;
    fpr=sum((pred==1)&(data(:,4)==0))/nn;
    rocmat(i,:)=[fpr tpr];
end
scatter(rocmat(:,1),rocmat(:,2))

% Try increasing or decreasing the value of sd at the start of the script to make the data less or more perfectly classifiable by this method to see how this affects the ROC curve.
%
% Now use a test set generated in the same way for a more fair assessment.

xnt=mvnrnd([-1,0],sd*eye(2),nn); % generate test data from the same two multivariate normal distributions
xpt=mvnrnd([1,1],sd*eye(2),np);
data=[zeros(nn,4);ones(np,4)];
data(:,1:2)=[xnt;xpt];
data(:,3)=ones(nn+np,1);
q=data(:,1:3)*copt;
prob=exp(q)./(1+exp(q));
pred=round(prob);
confusionchart(confusionmat(data(:,4),pred,'Order',[0 1]),char('n','p'))
acc_test=sum(pred==data(:,4))/(nn+np)

%
% The below sections perform the same example, but with the two normals having different variances., then illustrates a nonlinear logistic regression that includes quadratic terms

nn=1000;
np=1000;
sd=1;
xn=mvnrnd([-1,0],sd/2*eye(2),nn); % example using data from two multivariate normal distributions with different covariances matrices
xp=mvnrnd([1,1],sd*eye(2),np);
scatter(xn(:,1),xn(:,2)) % visualize the data
hold on;
scatter(xp(:,1),xp(:,2))

data=[zeros(nn,4);ones(np,4)];
data(:,1:2)=[xn;xp];
data(:,3)=ones(nn+np,1);
loglike=@(c) sum(log(1+exp(data(:,1:3)*c)))-data(:,4)'*(data(:,1:3)*c);
copt=fmincon(loglike,[0;0;0],[],[]);
x1max=max([xn(:,1);xp(:,1)]);
x1min=min([xn(:,1);xp(:,1)]);
fplot(@(x) (-(copt(3)+copt(1)*x))/copt(2),[x1min,x1max])
hold off;

q=data(:,1:3)*copt;
prob=exp(q)./(1+exp(q));
pred=round(prob);
confusionchart(confusionmat(data(:,4),pred,'Order',[0 1]),char('n','p'))
acc_train=sum(pred==data(:,4))/(nn+np)
taus=linspace(0,1,1000);
rocmat=zeros(1000,2);
for i=1:1000
    pred=(prob>=taus(i));
    tpr=sum((pred==1)&(data(:,4)==1))/np;
    fpr=sum((pred==1)&(data(:,4)==0))/nn;
    rocmat(i,:)=[fpr tpr];
end
scatter(rocmat(:,1),rocmat(:,2))

xnt=mvnrnd([-1,0],sd/2*eye(2),nn); % example using data from three different multivariate normal distributions
xpt=mvnrnd([1,1],sd*eye(2),np);
data=[zeros(nn,4);ones(np,4)];
data(:,1:2)=[xnt;xpt];
data(:,3)=ones(nn+np,1);
q=data(:,1:3)*copt;
prob=exp(q)./(1+exp(q));
pred=round(prob);
confusionchart(confusionmat(data(:,4),pred,'Order',[0 1]),char('n','p'))
acc_test=sum(pred==data(:,4))/(nn+np)



%
% The below section perform the same example as the one directly above, but illustrates a nonlinear logistic regression that includes quadratic terms.

nn=1000;
np=1000;
sd=1;
xn=mvnrnd([-1,0],sd/2*eye(2),nn); % example using data from two multivariate normal distributions with different covariances matrices
xp=mvnrnd([1,1],sd*eye(2),np);
scatter(xn(:,1),xn(:,2)) % visualize the data
hold on;
scatter(xp(:,1),xp(:,2))

data=[zeros(nn,6);ones(np,6)];
data(:,1:2)=[xn;xp];
data(:,3:4)=[xn.^2;xp.^2];
data(:,5)=ones(nn+np,1);
loglike=@(c) sum(log(1+exp(data(:,1:5)*c)))-data(:,6)'*(data(:,1:5)*c);
copt=fmincon(loglike,[0;0;0;0;0],[],[]);
x1max=max([xn(:,1);xp(:,1)]);
x1min=min([xn(:,1);xp(:,1)]);
x2max=max([xn(:,2);xp(:,2)]);
x2min=min([xn(:,2);xp(:,2)]);
xplot=linspace(x1min,x1max);
yplot=linspace(x2min,x2max);
[X,Y]=meshgrid(xplot,yplot);
Z=copt(1)*X+copt(2)*Y+copt(3)*X.^2+copt(4)*Y.^2+copt(5)*ones(100);
contour(X,Y,Z,[0 0])
hold off
q=data(:,1:5)*copt;
prob=exp(q)./(1+exp(q));
pred=round(prob);
confusionchart(confusionmat(data(:,6),pred,'Order',[0 1]),char('n','p'))
acc_train=sum(pred==data(:,6))/(nn+np)