# Calculating the probability of hitting, Fisher linear discriminant

I have 4 data matrix, which are described two parameters. (M1, M2, M3, M4). So I made transition from base parameters to principal components and displayed them in the space of new signs. I need help in calculating the probability of hitting a value in an area, example M1 (red triangles on picture). And how to calculate the straight line separating the plane? I want to separate the M1 from the other matrices, so only two classes can be considered, example M1 and M2. (M1 - values with no defects, others - with defects). Example, I want say, that if y1 and y2 will 0 and -0.5, then the probability that the object has no defect is 90%.

M1={{0.07, 0.09}, {0.09, 0.04}, {0.06, 0.5}, {0.05, 0.1}, {0.08,
0.07}, {0.1, 0.05}, {0.12, 0.14}, {0.09, 0.12}, {0.1, 0.07}, {0.04,
0.1}, {0.06, 0.1}, {0.09, 0.3}}
M2={{0.330347, 0.66736}, {0.283692, 0.59829}, {0.230995,
0.549643}, {0.353132, 0.711974}, {0.19679, 0.360885}, {0.268241,
0.505476}, {0.625252, 0.87299}, {0.732484, 0.939535}, {0.499508,
0.728917}, {0.646746, 0.870734}, {0.684303, 0.818975}, {0.448433,
0.644385}}
M3={{0.569952, 0.838573}, {0.451936, 0.801501}, {0.228608,
0.547101}, {0.161797, 0.345203}, {0.116077, 0.235871}, {0.124553,
0.249579}, {0.115536, 0.233847}, {0.114552, 0.23632}, {0.126844,
0.269214}, {0.158276, 0.340853}, {0.204667, 0.474522}, {0.263137,
0.576838}}
M4={{0.136822, 0.317039}, {0.198275, 0.482787}, {0.301204,
0.606273}, {0.424842, 0.540176}, {0.222032, 0.355069}, {0.12944,
0.2534}, {0.192282, 0.487613}, {0.230306, 0.517964}, {0.155553,
0.30778}, {0.205009, 0.476323}, {0.211225, 0.496847}, {0.146949,
0.255088}}
data={{0.330347, 0.66736}, {0.283692, 0.59829}, {0.230995,
0.549643}, {0.353132, 0.711974}, {0.19679, 0.360885}, {0.268241,
0.505476}, {0.625252, 0.87299}, {0.732484, 0.939535}, {0.499508,
0.728917}, {0.646746, 0.870734}, {0.684303, 0.818975}, {0.448433,
0.644385}, {0.569952, 0.838573}, {0.451936, 0.801501}, {0.228608,
0.547101}, {0.161797, 0.345203}, {0.116077, 0.235871}, {0.124553,
0.249579}, {0.115536, 0.233847}, {0.114552, 0.23632}, {0.126844,
0.269214}, {0.158276, 0.340853}, {0.204667, 0.474522}, {0.263137,
0.576838}, {0.136822, 0.317039}, {0.198275, 0.482787}, {0.301204,
0.606273}, {0.424842, 0.540176}, {0.222032, 0.355069}, {0.12944,
0.2534}, {0.192282, 0.487613}, {0.230306, 0.517964}, {0.155553,
0.30778}, {0.205009, 0.476323}, {0.211225, 0.496847}, {0.146949,
0.255088}, {0.07, 0.09}, {0.09, 0.04}, {0.06, 0.5}, {0.05,
0.1}, {0.08, 0.07}, {0.1, 0.05}, {0.12, 0.14}, {0.09, 0.12}, {0.1,
0.07}, {0.04, 0.1}, {0.06, 0.1}, {0.09, 0.3}}
builtinscores = PrincipalComponents[data, Method -> "Correlation"];
manscores =
Standardize[data].Transpose[Eigenvectors[Correlation[data]]];
ListPlot[Partition[builtinscores[[All, 1 ;; 2]], 12], Frame -> True,
PlotMarkers -> {Automatic, 20}, LabelStyle -> Directive[Black]]


• Sorry, my bad. I added the numerical data. – Aleksei Zhidkov Feb 10 '19 at 16:30
• @HenrikSchumacher I agree, sorry, my bad. I edited text. – Aleksei Zhidkov Feb 10 '19 at 16:41
• It seems you might want quadratic discriminant analysis as the covariance matrices of each of the 4 groups are different. So given that there are only two predictor variables the only thing that discriminant analysis does for you is to reduce the white space in a plot of the data. Also, what are the prior probabilities for each of the 4 groups? And even then the calculated probabilities are constructed on the assumption that each group has a bivariate normal distribution and that assumption does not look likely given the plot of the data. – JimB Feb 10 '19 at 17:52
• @JimB In fact, I want to separate the M1 from the other matrices, so only two classes can be considered, example M1 and M2. (M1 - values with no defects, others - with defects). But I do not know where to start. For example, if you set the probability of hitting the value range M1 or calculate it. Example, I want say, that if y1 and y2 will 0 and -0.5, then the probability that the object has no defect is 90%. – Aleksei Zhidkov Feb 10 '19 at 18:04
• While Mathematica can certainly do all the necessary calculations for linear and quadratic discriminant analysis (and nonparametric discriminant analysis), I am not aware of a specific function dedicated to doing so (such as NonlinearModelFit is for nonlinear regression). Given that so many other packages (such as R) have such built-in functionality, I'd use one of those packages rather than Mathematica especially if you are a relatively new user of Mathematica. (You should also add into your question - rather than in a comment - what you say about M1 vs the other groups.) – JimB Feb 10 '19 at 18:13

I still think that using other software that explicitly handles discriminant analysis (both linear and quadratic) is best as there are usually options for unequal priors, measures of precision, bias-removal, etc. However, the basic commands in Mathematica can do much of the work.

M1 = {{0.07, 0.09}, {0.09, 0.04}, {0.06, 0.5}, {0.05, 0.1}, {0.08,
0.07}, {0.1, 0.05}, {0.12, 0.14}, {0.09, 0.12}, {0.1,
0.07}, {0.04, 0.1}, {0.06, 0.1}, {0.09, 0.3}};
M2 = {{0.330347, 0.66736}, {0.283692, 0.59829}, {0.230995,
0.549643}, {0.353132, 0.711974}, {0.19679, 0.360885}, {0.268241,
0.505476}, {0.625252, 0.87299}, {0.732484, 0.939535}, {0.499508,
0.728917}, {0.646746, 0.870734}, {0.684303, 0.818975}, {0.448433,
0.644385}};
M3 = {{0.569952, 0.838573}, {0.451936, 0.801501}, {0.228608,
0.547101}, {0.161797, 0.345203}, {0.116077, 0.235871}, {0.124553,
0.249579}, {0.115536, 0.233847}, {0.114552, 0.23632}, {0.126844,
0.269214}, {0.158276, 0.340853}, {0.204667, 0.474522}, {0.263137,
0.576838}};
M4 = {{0.136822, 0.317039}, {0.198275, 0.482787}, {0.301204,
0.606273}, {0.424842, 0.540176}, {0.222032, 0.355069}, {0.12944,
0.2534}, {0.192282, 0.487613}, {0.230306, 0.517964}, {0.155553,
0.30778}, {0.205009, 0.476323}, {0.211225, 0.496847}, {0.146949,
0.255088}};

data = Join[M1, M2, M3, M4];
groups = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};


Find the log of the likelihood for each group:

logL = LogLikelihood[MultinormalDistribution[Mean[#], Covariance[#]], {{y1, y2}}] & /@ {M1, M2, M3, M4};


If the priors are equal, then the table of predicted posterior probabilities for the sample points is given by

logTotal = Log[Exp[logL[[1]]] + Exp[logL[[2]]] + Exp[logL[[3]]] + Exp[logL[[4]]]];
p = Exp[# - logTotal] & /@ logL;
Style[TableForm[p /. {y1 -> #[[1]], y2 -> #[[2]]} & /@ data,
TableHeadings -> {groups, {1, 2, 3, 4}}]]


A plot of the regions:

RegionPlot[{p[[1]] > Max[p[[{2, 3, 4}]]], p[[2]] > Max[p[[{1, 3, 4}]]],
p[[3]] > Max[p[[{1, 2, 4}]]], p[[4]] > Max[p[[{1, 2, 3}]]]},
{y1, 0, 1}, {y2, 0, 1}, PlotPoints -> 100, PlotRangePadding -> None,
PlotStyle -> ColorData[97, "ColorList"][[{1, 2, 3, 4}]],
PlotLegends -> {"1", "2", "3", "4"}, AspectRatio -> 1]


• Thank you! When I have time, I will definitely look other programs for quadratic discriminant analysis. If I find other solutions, i will show them. – Aleksei Zhidkov Feb 12 '19 at 7:22

1. You are responsible for knowing and vetting the necessary assumptions.
2. One can only expect so much from just 12 data points per group (irrespective as to how expensive it was just to get 12 observations).
3. Some of the 4 groups have large deviations from the assumption of bivariate normal distributions. (You can potentially get around this assumption by performing nonparametric discriminant analysis but you're still limited by just having 12 observations per group.)
4. The spread and orientation of the data points differ by group so even if the bivariate normal distribution held, you'd need to perform quadratic discriminant analysis as opposed to linear discriminant analysis.
5. Predictions should only be made in the range of the observed data (unless maybe you have some associated measure of precision that shows less precision as you move farther away form the observed data).

OK. Enough sermonizing.

Mathematica does not have a specialized function to perform linear or quadratic discriminant analysis (as there is for nonlinear regression with NonlinearModelFit). So I would suggest using R directly or using Mathematica's RLink functionality.

Here's how to get the posterior probabilities (assuming that bivariate normal distributions and equal prior probabilities are appropriate assumptions):

First give the sample data and groupings:

M1 = {{0.07, 0.09}, {0.09, 0.04}, {0.06, 0.5}, {0.05, 0.1}, {0.08,
0.07}, {0.1, 0.05}, {0.12, 0.14}, {0.09, 0.12}, {0.1,
0.07}, {0.04, 0.1}, {0.06, 0.1}, {0.09, 0.3}};
M2 = {{0.330347, 0.66736}, {0.283692, 0.59829}, {0.230995,
0.549643}, {0.353132, 0.711974}, {0.19679, 0.360885}, {0.268241,
0.505476}, {0.625252, 0.87299}, {0.732484, 0.939535}, {0.499508,
0.728917}, {0.646746, 0.870734}, {0.684303, 0.818975}, {0.448433,
0.644385}};
M3 = {{0.569952, 0.838573}, {0.451936, 0.801501}, {0.228608,
0.547101}, {0.161797, 0.345203}, {0.116077, 0.235871}, {0.124553,
0.249579}, {0.115536, 0.233847}, {0.114552, 0.23632}, {0.126844,
0.269214}, {0.158276, 0.340853}, {0.204667, 0.474522}, {0.263137,
0.576838}};
M4 = {{0.136822, 0.317039}, {0.198275, 0.482787}, {0.301204,
0.606273}, {0.424842, 0.540176}, {0.222032, 0.355069}, {0.12944,
0.2534}, {0.192282, 0.487613}, {0.230306, 0.517964}, {0.155553,
0.30778}, {0.205009, 0.476323}, {0.211225, 0.496847}, {0.146949,
0.255088}};
data = Join[M1, M2, M3, M4];

groups = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};


Here's a plot of the data:

ListPlot[{M1, M2, M3, M4}, PlotStyle -> PointSize[0.015], PlotLegends -> {"1", "2", "3", "4"},
AxesLabel -> {"\!$$\*SubscriptBox[\(y$$, $$1$$]\)", "\!$$\*SubscriptBox[\(y$$, $$2$$]\)"},
PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> 1]


Needs["RLink"]
InstallR["RHomeLocation" -> "c:\\Program Files\\R\\R-3.2.2"]


The version that R you can use depends on the version of Mathematica you have. (And I think setting the particular version of R only works on Windows - at least if you are using Mathematica 10.4.)

RSet["groups", groups];
RSet["d", data];


Run the quadratic discriminant analysis function from the R package MASS:

REvaluate["library(MASS)"];
REvaluate["x = matrix(d, ncol=2)"];
REvaluate["results = qda(x, groups)"];
samplePredictions =
Style[TableForm[REvaluate["predict(results, x)"][[1, 2, 1]],
TableHeadings -> {groups, {1, 2, 3, 4}}], Tiny]


If you want a specific prediction say at y1 = 0.25 and y2 = 0.5:

samplePredictions = REvaluate["predict(results, matrix(c(0.25,0.5), nrow=1, ncol=2))"][[1, 2, 1]]
(* {0.16935928707934308,0.371194179362576,0.45944653355807963,1.2147055565174368*^-15}} *)


Making predictions at a fine grid of values can get you a map as to which group has the maximum posterior probability:

grid = Flatten[Table[{y1, y2}, {y1, 0, 1, 0.005}, {y2, 0, 1, 0.005}], 1];
RSet["grid", grid];
gridPredictions = Transpose[{grid, REvaluate["predict(results, grid)"][[1, 2, 1]]}];
r1 = Select[gridPredictions, #[[2, 1]] > Max[#[[2, {2, 3, 4}]]] &][[All, 1]];
r2 = Select[gridPredictions, #[[2, 2]] > Max[#[[2, {1, 3, 4}]]] &][[All, 1]];
r3 = Select[gridPredictions, #[[2, 3]] > Max[#[[2, {1, 2, 4}]]] &][[All, 1]];
r4 = Select[gridPredictions, #[[2, 4]] > Max[#[[2, {1, 2, 3}]]] &][[All, 1]];
ListPlot[{r1, r2, r3, r4}, AspectRatio -> 1, PlotStyle -> PointSize[0.018],
PlotLegends -> {1, 2, 3, 4}, PlotLabel -> "Group with largest posterior probability",
AxesLabel -> {"\!$$\*SubscriptBox[\(y$$, $$1$$]\)", "\!$$\*SubscriptBox[\(y$$, $$2$$]\)"}]
`

• Wow, thank you very much. Can i apply this algorithm for 3 parameters (M={x1,x2,x3})? – Aleksei Zhidkov Feb 11 '19 at 19:58