How to plot logistic regression decision boundary? I am running logistic regression on a small dataset which looks like this:



After implementing gradient descent and the cost function, I am getting a 100% accuracy in the prediction stage, However I want to be sure that everything is in order so I am trying to plot the decision boundary line which separates the two datasets.

Below I present plots showing the cost function and theta parameters. As can be seen, currently I am printing the decision boundary line incorrectly.



Extracting data

clear all; close all; clc;

alpha = 0.01;
num_iters = 1000;

%% Plotting data
x1 = linspace(0,3,50);
mqtrue = 5;
cqtrue = 30;
dat1 = mqtrue*x1+5*randn(1,50);

x2 = linspace(7,10,50);
dat2 = mqtrue*x2 + (cqtrue + 5*randn(1,50));

x = [x1 x2]'; % X

subplot(2,2,1);
dat = [dat1 dat2]'; % Y

scatter(x1, dat1); hold on;
scatter(x2, dat2, '*'); hold on;
classdata = (dat&gt;40);


Computing Cost, Gradient and plotting

%  Setup the data matrix appropriately, and add ones for the intercept term
[m, n] = size(x);

% Add intercept term to x and X_test
x = [ones(m, 1) x];

% Initialize fitting parameters
theta = zeros(n + 1, 1);
%initial_theta = [0.2; 0.2];

J_history = zeros(num_iters, 1); 

plot_x = [min(x(:,2))-2,  max(x(:,2))+2]

for iter = 1:num_iters 
% Compute and display initial cost and gradient
    [cost, grad] = logistic_costFunction(theta, x, classdata);
    theta = theta - alpha * grad;
    J_history(iter) = cost;

    fprintf('Iteration #%d - Cost = %d... \r\n',iter, cost);


    subplot(2,2,2);
    hold on; grid on;
    plot(iter, J_history(iter), '.r');  title(sprintf('Plot of cost against number of iterations. Cost is %g',J_history(iter)));
    xlabel('Iterations')
    ylabel('MSE')
    drawnow

    subplot(2,2,3);
    grid on;
    plot3(theta(1), theta(2), J_history(iter),'o')
    title(sprintf('Tita0 = %g, Tita1=%g', theta(1), theta(2)))
    xlabel('Tita0')
    ylabel('Tita1')
    zlabel('Cost')
    hold on;
    drawnow

    subplot(2,2,1);
    grid on;    
    % Calculate the decision boundary line
    plot_y = theta(2).*plot_x + theta(1);  % &lt;--- Boundary line 
    % Plot, and adjust axes for better viewing
    plot(plot_x, plot_y)
    hold on;
    drawnow

end

fprintf('Cost at initial theta (zeros): %f\n', cost);
fprintf('Gradient at initial theta (zeros): \n');
fprintf(' %f \n', grad);


The above code is implementing gradient descent correctly (I think) but I am still unable to show the boundary line plot. Any suggestions would be appreciated.



logistic_costFunction.m

function [J, grad] = logistic_costFunction(theta, X, y)

    % Initialize some useful values
    m = length(y); % number of training examples

    grad = zeros(size(theta));

    h = sigmoid(X * theta);
    J = -(1 / m) * sum( (y .* log(h)) + ((1 - y) .* log(1 - h)) );

    for i = 1 : size(theta, 1)
        grad(i) = (1 / m) * sum( (h - y) .* X(:, i) );
    end

end


EDIT:

As per the below answer by @Esmailian, now I have something like this:

[m, n] = size(x);

x1_class = [ones(m, 1) x1' dat1'];
x2_class = [ones(m, 1) x2' dat2'];

x = [x1_class ; x2_class]
Your decision boundary is a surface in 3D as your points are in 2D.

With Wolfram Language

Create the data sets.

mqtrue = 5;
cqtrue = 30;
With[{x = Subdivide[0, 3, 50]},
  dat1 = Transpose@{x, mqtrue x + 5 RandomReal[1, Length@x]};
  ];
With[{x = Subdivide[7, 10, 50]},
  dat2 = Transpose@{x, mqtrue x + cqtrue + 5 RandomReal[1, Length@x]};
  ];


View in 2D (ListPlot) and the 3D (ListPointPlot3D) regression space.

ListPlot[{dat1, dat2}, PlotMarkers -&gt; "OpenMarkers", PlotTheme -&gt; "Detailed"]



  


I Append the response variable to the data. 

datPlot =
 ListPointPlot3D[
  MapThread[Append, {#, Boole@Thread[#[[All, 2]] &gt; 40]}] &amp; /@ {dat1, dat2}
  ]



  


Perform a Logistic regression (LogitModelFit). You could use GeneralizedLinearModelFit with ExponentialFamily set to "Binomial" as well.

With[{dat = Join[dat1, dat2]},
 model =
  LogitModelFit[
   MapThread[Append, {dat, Boole@Thread[dat[[All, 2]] &gt; 40]}],
   {x, y}, {x, y}]
 ]



  


From the FittedModel "Properties" we need "Function".

model["Properties"]



  {AdjustedLikelihoodRatioIndex,    DevianceTableDeviances, ParameterConfidenceIntervalTableEntries,
  AIC,    DevianceTableEntries,   ParameterConfidenceRegion,
  AnscombeResiduals,  DevianceTableResidualDegreesOfFreedom,  ParameterErrors,
  BasisFunctions, DevianceTableResidualDeviances, ParameterPValues,
  BestFit,    EfronPseudoRSquared,    ParameterTable,
  BestFitParameters,  EstimatedDispersion,    ParameterTableEntries,
  BIC,    FitResiduals,   ParameterZStatistics,
  CookDistances,  Function,   PearsonChiSquare,
  CorrelationMatrix,  HatDiagonal,    PearsonResiduals,
  CovarianceMatrix,   LikelihoodRatioIndex,   PredictedResponse,
  CoxSnellPseudoRSquared, LikelihoodRatioStatistic,   Properties,
  CraggUhlerPseudoRSquared,   LikelihoodResiduals,    ResidualDeviance,
  Data,   LinearPredictor,    ResidualDegreesOfFreedom,
  DesignMatrix,   LogLikelihood,  Response,
  DevianceResiduals,  NullDeviance,   StandardizedDevianceResiduals,
  Deviances,  NullDegreesOfFreedom,   StandardizedPearsonResiduals,
  DevianceTable,  ParameterConfidenceIntervals,   WorkingResiduals,
  DevianceTableDegreesOfFreedom,  ParameterConfidenceIntervalTable}   


model["Function"]



  


Use this for prediction 

model["Function"][8, 54]



0.0196842



and plot the decision boundary surface in 3D along with the data (datPlot) using Show and Plot3D

modelPlot =
 Show[
  datPlot,
  Plot3D[
   model["Function"][x, y],
   Evaluate[
    Sequence @@ 
     MapThread[Prepend, {MinMax /@ Transpose@Join[dat1, dat2], {x, y}}]],
   Mesh -&gt; None,
   PlotStyle -&gt; Opacity[.25, Green],
   PlotPoints -&gt; 30
   ]
  ]




With ParametricPlot3D and Manipulate you can examine decision boundary curves for values of the variables. For example, keeping x fixed and letting y vary or vice versa.

Manipulate[
 Show[
  modelPlot,
  ParametricPlot3D[
   {x, u, model["Function"][x, u]}, {u, 0, 80}, PlotStyle -&gt; Orange],
  ParametricPlot3D[
   {u, y, model["Function"][u, y]}, {u, 0, 10}, PlotStyle -&gt; Purple],
  PlotLabel -&gt; 
   StringTemplate["model[`1`, `2`] = `3`"] @@ {x, y, model["Function"][x, y]}
  ],
 {{x, 6, Style["x", Orange, Bold]}, 0, 10, Appearance -&gt; "Labeled"},
 {{y, 40, Style["y", Purple, Bold]}, 0, 80, Appearance -&gt; "Labeled"}
 ]




Update

You can also plot contours of the probability in 2D.

plot = ListPlot[{dat1, dat2}, PlotMarkers -&gt; "OpenMarkers", PlotTheme -&gt; "Detailed"];

Manipulate[
 db = y /. First@Quiet@Solve[model["Function"][x, y] == p, y];
 Show[
  plot,
  Plot[db, {x, 0, 10}, PlotStyle -&gt; Red]
  ],
 {{p, .5}, 0, 1, Appearance -&gt; "Labeled"}
 ]




Hope this helps.