<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/classification-from-scratch-linear-discrimination-8-8/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://freakonometrics.hypotheses.org/53021"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Eighth post of our series on <a rel="nofollow" href="https://freakonometrics.hypotheses.org/52731" target="_blank">classification from scratch</a>. The latest one was on the <a rel="nofollow" href="https://freakonometrics.hypotheses.org/527744" target="_blank">SVM</a>, and today, I want to get back on very old stuff, with here also a linear separation of the space, using Fisher’s linear discriminent analysis.</p>
<h2>Bayes (naive) classifier</h2>
<p>Consider the follwing naive classification rule<img alt="m^\star(\mathbf{x})=\text{argmin}_y\{\mathbb{P}[Y=y\vert\mathbf{X}=\mathbf{x}]\}" title="m^\star(\mathbf{x})=\text{argmin}_y\{\mathbb{P}[Y=y\vert\mathbf{X}=\mathbf{x}]\}" src="https://s0.wp.com/latex.php?latex=m%5E%5Cstar%28%5Cmathbf%7Bx%7D%29%3D%5Ctext%7Bargmin%7D_y%5C%7B%5Cmathbb%7BP%7D%5BY%3Dy%5Cvert%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%5D%5C%7D&bg=ffffff&fg=000&s=0" class="latex">or<img alt="m^\star(\mathbf{x})=\text{argmin}_y\left\{\frac{\mathbb{P}[\mathbf{X}=\mathbf{x}\vert Y=y]}{\mathbb{P}[\mathbf{X}=\mathbf{x}]}\right\}" title="m^\star(\mathbf{x})=\text{argmin}_y\left\{\frac{\mathbb{P}[\mathbf{X}=\mathbf{x}\vert Y=y]}{\mathbb{P}[\mathbf{X}=\mathbf{x}]}\right\}" src="https://s0.wp.com/latex.php?latex=m%5E%5Cstar%28%5Cmathbf%7Bx%7D%29%3D%5Ctext%7Bargmin%7D_y%5Cleft%5C%7B%5Cfrac%7B%5Cmathbb%7BP%7D%5B%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%5Cvert+Y%3Dy%5D%7D%7B%5Cmathbb%7BP%7D%5B%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%5D%7D%5Cright%5C%7D&bg=ffffff&fg=000&s=0" class="latex">(where <img alt="\mathbb{P}[\mathbf{X}=\mathbf{x}]" title="\mathbb{P}[\mathbf{X}=\mathbf{x}]" src="https://s0.wp.com/latex.php?latex=%5Cmathbb%7BP%7D%5B%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%5D&bg=ffffff&fg=000&s=0" class="latex"> is the density in the continuous case).</p>
<p>In the case where <img alt="y" title="y" src="https://s0.wp.com/latex.php?latex=y&bg=ffffff&fg=000&s=0" class="latex"> takes two values, that will be standard <img alt="\{0,1\}" title="\{0,1\}" src="https://s0.wp.com/latex.php?latex=%5C%7B0%2C1%5C%7D&bg=ffffff&fg=000&s=0" class="latex"> here, one can rewrite the later as<img alt="m^\star(\mathbf{x})=\begin{cases}1\text{ if }\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})&gt;\displaystyle{\frac{1}{2}}\\0\text{ otherwise}\end{cases}" title="m^\star(\mathbf{x})=\begin{cases}1\text{ if }\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})&gt;\displaystyle{\frac{1}{2}}\\0\text{ otherwise}\end{cases}" src="https://s0.wp.com/latex.php?latex=m%5E%5Cstar%28%5Cmathbf%7Bx%7D%29%3D%5Cbegin%7Bcases%7D1%5Ctext%7B+if+%7D%5Cmathbb%7BE%7D%28Y%5Cvert+%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%29%3E%5Cdisplaystyle%7B%5Cfrac%7B1%7D%7B2%7D%7D%5C%5C0%5Ctext%7B+otherwise%7D%5Cend%7Bcases%7D&bg=ffffff&fg=000&s=0" class="latex">and the set<img alt="\mathcal{D}_S =\left\{\mathbf{x},\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})=\frac{1}{2}\right\}" title="\mathcal{D}_S =\left\{\mathbf{x},\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})=\frac{1}{2}\right\}" src="https://s0.wp.com/latex.php?latex=%5Cmathcal%7BD%7D_S+%3D%5Cleft%5C%7B%5Cmathbf%7Bx%7D%2C%5Cmathbb%7BE%7D%28Y%5Cvert+%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%29%3D%5Cfrac%7B1%7D%7B2%7D%5Cright%5C%7D&bg=ffffff&fg=000&s=0" class="latex">is called the decision boundary.</p>
<p>Assume that<img alt="\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})" title="\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7BX%7D%5Cvert+Y%3D0%5Csim%5Cmathcal%7BN%7D%28%5Cmathbf%7B%5Cmu%7D_0%2C%5Cmathbf%7B%5CSigma%7D%29&bg=ffffff&fg=000&s=0" class="latex">and<img alt="\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})" title="\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7BX%7D%5Cvert+Y%3D1%5Csim%5Cmathcal%7BN%7D%28%5Cmathbf%7B%5Cmu%7D_1%2C%5Cmathbf%7B%5CSigma%7D%29&bg=ffffff&fg=000&s=0" class="latex">then explicit expressions can be derived.<img alt="m^\star(\mathbf{x})=\begin{cases}1\text{ if }r_1^2&lt; r_0^2+2\displaystyle{\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}+\log\frac{\vert\mathbf{\Sigma}_0\vert}{\vert\mathbf{\Sigma}_1\vert}}\\0\text{ otherwise}\end{cases}[/latex]where [latex]r_y^2[/latex] is the Manalahobis distance, [latex display=&quot;true&quot;]r_y^2 = [\mathbf{X}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[\mathbf{X}-\mathbf{\mu}_y][/latex]&lt;/p&gt; &lt;p&gt;Let [latex]\delta_y" title="m^\star(\mathbf{x})=\begin{cases}1\text{ if }r_1^2&lt; r_0^2+2\displaystyle{\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}+\log\frac{\vert\mathbf{\Sigma}_0\vert}{\vert\mathbf{\Sigma}_1\vert}}\\0\text{ otherwise}\end{cases}[/latex]where [latex]r_y^2[/latex] is the Manalahobis distance, [latex display=&quot;true&quot;]r_y^2 = [\mathbf{X}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[\mathbf{X}-\mathbf{\mu}_y][/latex]&lt;/p&gt; &lt;p&gt;Let [latex]\delta_y" src="https://s0.wp.com/latex.php?latex=m%5E%5Cstar%28%5Cmathbf%7Bx%7D%29%3D%5Cbegin%7Bcases%7D1%5Ctext%7B+if+%7Dr_1%5E2%3C+r_0%5E2%2B2%5Cdisplaystyle%7B%5Clog%5Cfrac%7B%5Cmathbb%7BP%7D%28Y%3D1%29%7D%7B%5Cmathbb%7BP%7D%28Y%3D0%29%7D%2B%5Clog%5Cfrac%7B%5Cvert%5Cmathbf%7B%5CSigma%7D_0%5Cvert%7D%7B%5Cvert%5Cmathbf%7B%5CSigma%7D_1%5Cvert%7D%7D%5C%5C0%5Ctext%7B+otherwise%7D%5Cend%7Bcases%7D%26%2391%3B%2Flatex%26%2393%3Bwhere+%26%2391%3Blatex%26%2393%3Br_y%5E2%26%2391%3B%2Flatex%26%2393%3B+is+the+Manalahobis+distance%2C+%26%2391%3Blatex+display%3D%22true%22%26%2393%3Br_y%5E2+%3D+%26%2391%3B%5Cmathbf%7BX%7D-%5Cmathbf%7B%5Cmu%7D_y%26%2393%3B%5E%7B%5Ctext%7B%7BT%7D%7D%7D%5Cmathbf%7B%5CSigma%7D_y%5E%7B-1%7D%26%2391%3B%5Cmathbf%7BX%7D-%5Cmathbf%7B%5Cmu%7D_y%26%2393%3B%26%2391%3B%2Flatex%26%2393%3B%3C%2Fp%3E+%3Cp%3ELet+%5Blatex%5D%5Cdelta_y&bg=ffffff&fg=000&s=0" class="latex">be defined as<img alt="\delta_y(\mathbf{x})=-\frac{1}{2}\log\vert\mathbf{\Sigma}_y\vert-\frac{1}{2}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]+\log\mathbb{P}(Y=y)" title="\delta_y(\mathbf{x})=-\frac{1}{2}\log\vert\mathbf{\Sigma}_y\vert-\frac{1}{2}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]+\log\mathbb{P}(Y=y)" src="https://s0.wp.com/latex.php?latex=%5Cdelta_y%28%5Cmathbf%7Bx%7D%29%3D-%5Cfrac%7B1%7D%7B2%7D%5Clog%5Cvert%5Cmathbf%7B%5CSigma%7D_y%5Cvert-%5Cfrac%7B1%7D%7B2%7D%5B%7B%5Ccolor%7Bblue%7D%7B%5Cmathbf%7Bx%7D%7D%7D-%5Cmathbf%7B%5Cmu%7D_y%5D%5E%7B%5Ctext%7B%7BT%7D%7D%7D%5Cmathbf%7B%5CSigma%7D_y%5E%7B-1%7D%5B%7B%5Ccolor%7Bblue%7D%7B%5Cmathbf%7Bx%7D%7D%7D-%5Cmathbf%7B%5Cmu%7D_y%5D%2B%5Clog%5Cmathbb%7BP%7D%28Y%3Dy%29&bg=ffffff&fg=000&s=0" class="latex">the decision boundary of this classifier is <img alt="\{\mathbf{x}\text{ such that }\delta_0(\mathbf{x})=\delta_1(\mathbf{x})\}" title="\{\mathbf{x}\text{ such that }\delta_0(\mathbf{x})=\delta_1(\mathbf{x})\}" src="https://s0.wp.com/latex.php?latex=%5C%7B%5Cmathbf%7Bx%7D%5Ctext%7B+such+that+%7D%5Cdelta_0%28%5Cmathbf%7Bx%7D%29%3D%5Cdelta_1%28%5Cmathbf%7Bx%7D%29%5C%7D&bg=ffffff&fg=000&s=0" class="latex">which is quadratic in <img alt="{\color{blue}{\mathbf{x}}}" title="{\color{blue}{\mathbf{x}}}" src="https://s0.wp.com/latex.php?latex=%7B%5Ccolor%7Bblue%7D%7B%5Cmathbf%7Bx%7D%7D%7D&bg=ffffff&fg=000&s=0" class="latex">. This is the quadratic discriminant analysis. This can be visualized bellow.<br>
<img src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-05-a%CC%80-21.59.58.png?resize=300%2C266&ssl=1" alt="" width="300" data-recalc-dims="1" height="266"><br>
The decision boundary is here<br>
<img src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-05-a%CC%80-22.00.10.png?resize=300%2C266&ssl=1" alt="" width="300" data-recalc-dims="1" height="266"><br>
But that can’t be the linear discriminant analysis, right? I mean, the frontier is not linear… Actually, in Fisher’s <a rel="nofollow" href="https://digital.library.adelaide.edu.au/dspace/bitstream/2440/15227/1/138.pdf" target="_blank">seminal paper</a>, it was assumed that <img alt="\mathbf{\Sigma}_0=\mathbf{\Sigma}_1" title="\mathbf{\Sigma}_0=\mathbf{\Sigma}_1" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5CSigma%7D_0%3D%5Cmathbf%7B%5CSigma%7D_1&bg=ffffff&fg=000&s=0" class="latex">.</p>
<p>In that case, actually, <img alt="\delta_y(\mathbf{x})={\color{blue}{\mathbf{x}}}^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y-\frac{1}{2}\mathbf{\mu}_y^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y+\log\mathbb{P}(Y=y)" title="\delta_y(\mathbf{x})={\color{blue}{\mathbf{x}}}^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y-\frac{1}{2}\mathbf{\mu}_y^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y+\log\mathbb{P}(Y=y)" src="https://s0.wp.com/latex.php?latex=%5Cdelta_y%28%5Cmathbf%7Bx%7D%29%3D%7B%5Ccolor%7Bblue%7D%7B%5Cmathbf%7Bx%7D%7D%7D%5E%7B%5Ctext%7BT%7D%7D%5Cmathbf%7B%5CSigma%7D%5E%7B-1%7D%5Cmathbf%7B%5Cmu%7D_y-%5Cfrac%7B1%7D%7B2%7D%5Cmathbf%7B%5Cmu%7D_y%5E%7B%5Ctext%7BT%7D%7D%5Cmathbf%7B%5CSigma%7D%5E%7B-1%7D%5Cmathbf%7B%5Cmu%7D_y%2B%5Clog%5Cmathbb%7BP%7D%28Y%3Dy%29&bg=ffffff&fg=000&s=0" class="latex"> and the decision frontier is now linear in <img alt="{\color{blue}{\mathbf{x}}}" title="{\color{blue}{\mathbf{x}}}" src="https://s0.wp.com/latex.php?latex=%7B%5Ccolor%7Bblue%7D%7B%5Cmathbf%7Bx%7D%7D%7D&bg=ffffff&fg=000&s=0" class="latex">. This is the linear discriminant analysis. This can be visualized bellow<br>
<img src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-05-a%CC%80-21.59.28.png?resize=300%2C266&ssl=1" alt="" width="300" data-recalc-dims="1" height="266"><br>
Here the two samples have the same variance matrix and the frontier is<br>
<img src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-05-a%CC%80-21.59.45.png?resize=300%2C266&ssl=1" alt="" width="300" data-recalc-dims="1" height="266"></p>
<h2>Link with the logistic regression</h2>
<p>Assume as previously that<img alt="\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})" title="\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7BX%7D%5Cvert+Y%3D0%5Csim%5Cmathcal%7BN%7D%28%5Cmathbf%7B%5Cmu%7D_0%2C%5Cmathbf%7B%5CSigma%7D%29&bg=ffffff&fg=000&s=0" class="latex">and<img alt="\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})" title="\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7BX%7D%5Cvert+Y%3D1%5Csim%5Cmathcal%7BN%7D%28%5Cmathbf%7B%5Cmu%7D_1%2C%5Cmathbf%7B%5CSigma%7D%29&bg=ffffff&fg=000&s=0" class="latex">then<img alt="\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}" title="\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}" src="https://s0.wp.com/latex.php?latex=%5Clog%5Cfrac%7B%5Cmathbb%7BP%7D%28Y%3D1%5Cvert+%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%29%7D%7B%5Cmathbb%7BP%7D%28Y%3D0%5Cvert+%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%29%7D&bg=ffffff&fg=000&s=0" class="latex">is equal to <img alt="\mathbf{x}^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_y]-\frac{1}{2}[\mathbf{\mu}_1-\mathbf{\mu}_0]^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]+\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}" title="\mathbf{x}^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_y]-\frac{1}{2}[\mathbf{\mu}_1-\mathbf{\mu}_0]^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]+\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7Bx%7D%5E%7B%5Ctext%7B%7BT%7D%7D%7D%5Cmathbf%7B%5CSigma%7D%5E%7B-1%7D%5B%5Cmathbf%7B%5Cmu%7D_y%5D-%5Cfrac%7B1%7D%7B2%7D%5B%5Cmathbf%7B%5Cmu%7D_1-%5Cmathbf%7B%5Cmu%7D_0%5D%5E%7B%5Ctext%7B%7BT%7D%7D%7D%5Cmathbf%7B%5CSigma%7D%5E%7B-1%7D%5B%5Cmathbf%7B%5Cmu%7D_1-%5Cmathbf%7B%5Cmu%7D_0%5D%2B%5Clog%5Cfrac%7B%5Cmathbb%7BP%7D%28Y%3D1%29%7D%7B%5Cmathbb%7BP%7D%28Y%3D0%29%7D&bg=ffffff&fg=000&s=0" class="latex">which is linear in <img alt="\mathbf{x}" title="\mathbf{x}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7Bx%7D&bg=ffffff&fg=000&s=0" class="latex"><img alt="\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}=\mathbf{x}^{\text{{T}}}\mathbf{\beta}" title="\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}=\mathbf{x}^{\text{{T}}}\mathbf{\beta}" src="https://s0.wp.com/latex.php?latex=%5Clog%5Cfrac%7B%5Cmathbb%7BP%7D%28Y%3D1%5Cvert+%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%29%7D%7B%5Cmathbb%7BP%7D%28Y%3D0%5Cvert+%5Cmathbf%7BX%7D%3D%5Cmathbf%7Bx%7D%29%7D%3D%5Cmathbf%7Bx%7D%5E%7B%5Ctext%7B%7BT%7D%7D%7D%5Cmathbf%7B%5Cbeta%7D&bg=ffffff&fg=000&s=0" class="latex">Hence, when each groups have Gaussian distributions with identical variance matrix, then LDA and the logistic regression lead to the same classification rule.</p>
<p>Observe furthermore that the slope is proportional to <img alt="\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]" title="\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5CSigma%7D%5E%7B-1%7D%5B%5Cmathbf%7B%5Cmu%7D_1-%5Cmathbf%7B%5Cmu%7D_0%5D&bg=ffffff&fg=000&s=0" class="latex">, as stated in Fisher’s article. But to obtain such a relationship, he observe that the ratio of between and within variances (in the two groups) was<img alt="\frac{\text{variance between}}{\text{variance within}}=\frac{[\mathbf{\omega}\mathbf{\mu}_1-\mathbf{\omega}\mathbf{\mu}_0]^2}{\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_1\mathbf{\omega}+\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_0\mathbf{\omega}}" title="\frac{\text{variance between}}{\text{variance within}}=\frac{[\mathbf{\omega}\mathbf{\mu}_1-\mathbf{\omega}\mathbf{\mu}_0]^2}{\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_1\mathbf{\omega}+\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_0\mathbf{\omega}}" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B%5Ctext%7Bvariance+between%7D%7D%7B%5Ctext%7Bvariance+within%7D%7D%3D%5Cfrac%7B%5B%5Cmathbf%7B%5Comega%7D%5Cmathbf%7B%5Cmu%7D_1-%5Cmathbf%7B%5Comega%7D%5Cmathbf%7B%5Cmu%7D_0%5D%5E2%7D%7B%5Cmathbf%7B%5Comega%7D%5E%7B%5Ctext%7BT%7D%7D%5Cmathbf%7B%5CSigma%7D_1%5Cmathbf%7B%5Comega%7D%2B%5Cmathbf%7B%5Comega%7D%5E%7B%5Ctext%7BT%7D%7D%5Cmathbf%7B%5CSigma%7D_0%5Cmathbf%7B%5Comega%7D%7D&bg=ffffff&fg=000&s=0" class="latex">which is maximal when <img alt="\mathbf{\omega}" title="\mathbf{\omega}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Comega%7D&bg=ffffff&fg=000&s=0" class="latex"> is proportional to <img alt="\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]" title="\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5CSigma%7D%5E%7B-1%7D%5B%5Cmathbf%7B%5Cmu%7D_1-%5Cmathbf%7B%5Cmu%7D_0%5D&bg=ffffff&fg=000&s=0" class="latex">, when <img alt="\mathbf{\Sigma}_0=\mathbf{\Sigma}_1" title="\mathbf{\Sigma}_0=\mathbf{\Sigma}_1" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5CSigma%7D_0%3D%5Cmathbf%7B%5CSigma%7D_1&bg=ffffff&fg=000&s=0" class="latex">.</p>
<h2>Homebrew linear discriminant analysis</h2>
<p>To compute vector <img alt="\mathbf{\omega}" title="\mathbf{\omega}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Comega%7D&bg=ffffff&fg=000&s=0" class="latex"></p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
</pre></td><td><pre>m0 <span>=</span> <span>apply</span><span>(</span>myocarde<span>[</span>myocarde$PRONO<span>==</span><span>&quot;0&quot;</span>,<span>1</span><span>:</span><span>7</span><span>]</span>,<span>2</span>,<span>mean</span><span>)</span>
m1 <span>=</span> <span>apply</span><span>(</span>myocarde<span>[</span>myocarde$PRONO<span>==</span><span>&quot;1&quot;</span>,<span>1</span><span>:</span><span>7</span><span>]</span>,<span>2</span>,<span>mean</span><span>)</span>
Sigma <span>=</span> <span>var</span><span>(</span>myocarde<span>[</span>,<span>1</span><span>:</span><span>7</span><span>]</span><span>)</span>
omega <span>=</span> <span>solve</span><span>(</span>Sigma<span>)</span><span>%*%</span><span>(</span>m1<span>-</span>m0<span>)</span>
omega
                 <span>[</span>,<span>1</span><span>]</span>
FRCAR <span>-</span><span>0.012909708542</span>
INCAR  <span>1.088582058796</span>
INSYS <span>-</span><span>0.019390084344</span>
PRDIA <span>-</span><span>0.025817110020</span>
PAPUL  <span>0.020441287970</span>
PVENT <span>-</span><span>0.038298291091</span>
REPUL <span>-</span><span>0.001371677757</span></pre></td></tr></table></div>
<p>For the constant – in the equation <img alt="\omega^T\mathbf{x}+b=0" title="\omega^T\mathbf{x}+b=0" src="https://s0.wp.com/latex.php?latex=%5Comega%5ET%5Cmathbf%7Bx%7D%2Bb%3D0&bg=ffffff&fg=000&s=0" class="latex"> – if we have equiprobable probabilities, use</p>
<div><table><tr><td><pre>1
</pre></td><td><pre>b <span>=</span> <span>(</span><span>t</span><span>(</span>m1<span>)</span><span>%*%</span><span>solve</span><span>(</span>Sigma<span>)</span><span>%*%</span>m1<span>-</span><span>t</span><span>(</span>m0<span>)</span><span>%*%</span><span>solve</span><span>(</span>Sigma<span>)</span><span>%*%</span>m0<span>)</span><span>/</span><span>2</span></pre></td></tr></table></div>
<h2>Application (on the small dataset)</h2>
<p>In order to visualize what’s going on, consider the small dataset, with only two covariates,</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
</pre></td><td><pre>x <span>=</span> <span>c</span><span>(</span>.4,.55,.65,.9,.1,.35,.5,.15,.2,.85<span>)</span>
y <span>=</span> <span>c</span><span>(</span>.85,.95,.8,.87,.5,.55,.5,.2,.1,.3<span>)</span>
z <span>=</span> <span>c</span><span>(</span><span>1</span>,<span>1</span>,<span>1</span>,<span>1</span>,<span>1</span>,<span>0</span>,<span>0</span>,<span>1</span>,<span>0</span>,<span>0</span><span>)</span>
<span>df</span> <span>=</span> <span>data.<span>frame</span></span><span>(</span>x1<span>=</span>x,x2<span>=</span>y,y<span>=</span><span>as.<span>factor</span></span><span>(</span>z<span>)</span><span>)</span>
m0 <span>=</span> <span>apply</span><span>(</span><span>df</span><span>[</span><span>df</span>$y<span>==</span><span>&quot;0&quot;</span>,<span>1</span><span>:</span><span>2</span><span>]</span>,<span>2</span>,<span>mean</span><span>)</span>
m1 <span>=</span> <span>apply</span><span>(</span><span>df</span><span>[</span><span>df</span>$y<span>==</span><span>&quot;1&quot;</span>,<span>1</span><span>:</span><span>2</span><span>]</span>,<span>2</span>,<span>mean</span><span>)</span>
Sigma <span>=</span> <span>var</span><span>(</span><span>df</span><span>[</span>,<span>1</span><span>:</span><span>2</span><span>]</span><span>)</span>
omega <span>=</span> <span>solve</span><span>(</span>Sigma<span>)</span><span>%*%</span><span>(</span>m1<span>-</span>m0<span>)</span>
omega
         <span>[</span>,<span>1</span><span>]</span>
x1 <span>-</span><span>2.640613174</span>
x2  <span>4.858705676</span></pre></td></tr></table></div>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-06-a%CC%80-05.26.42.png?w=456&ssl=1" data-recalc-dims="1"><br>
Using R regular function, we get</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
</pre></td><td><pre><span>library</span><span>(</span>MASS<span>)</span>
fit_lda <span>=</span> lda<span>(</span>y ~x1<span>+</span>x2 , <span>data</span><span>=</span><span>df</span><span>)</span>
fit_lda
 
Coefficients of linear discriminants<span>:</span>
            LD1
x1 <span>-</span><span>2.588389554</span>
x2  <span>4.762614663</span></pre></td></tr></table></div>
<p>which is the same coefficient as the one we got with our own code. For the constant, use</p>
<div><table><tr><td><pre>1
</pre></td><td><pre>b <span>=</span> <span>(</span><span>t</span><span>(</span>m1<span>)</span><span>%*%</span><span>solve</span><span>(</span>Sigma<span>)</span><span>%*%</span>m1<span>-</span><span>t</span><span>(</span>m0<span>)</span><span>%*%</span><span>solve</span><span>(</span>Sigma<span>)</span><span>%*%</span>m0<span>)</span><span>/</span><span>2</span></pre></td></tr></table></div>
<p>If we plot it, we get the red straight line</p>
<div><table><tr><td><pre>1
2
</pre></td><td><pre><span>plot</span><span>(</span><span>df</span>$x1,<span>df</span>$x2,pch<span>=</span><span>c</span><span>(</span><span>1</span>,<span>19</span><span>)</span><span>[</span><span>1</span><span>+</span><span>(</span><span>df</span>$y<span>==</span><span>&quot;1&quot;</span><span>)</span><span>]</span><span>)</span>
<span>abline</span><span>(</span>a<span>=</span>b<span>/</span>omega<span>[</span><span>2</span><span>]</span>,b<span>=-</span>omega<span>[</span><span>1</span><span>]</span><span>/</span>omega<span>[</span><span>2</span><span>]</span>,<span>col</span><span>=</span><span>&quot;red&quot;</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-06-a%CC%80-16.11.32.png?w=456&ssl=1" data-recalc-dims="1"><br>
As we can see (with the blue points), our red line intersects the middle of the segment of the two barycenters</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre><span>points</span><span>(</span>m0<span>[</span><span>&quot;x1&quot;</span><span>]</span>,m0<span>[</span><span>&quot;x2&quot;</span><span>]</span>,pch<span>=</span><span>4</span><span>)</span>
<span>points</span><span>(</span>m1<span>[</span><span>&quot;x1&quot;</span><span>]</span>,m1<span>[</span><span>&quot;x2&quot;</span><span>]</span>,pch<span>=</span><span>4</span><span>)</span>
<span>segments</span><span>(</span>m0<span>[</span><span>&quot;x1&quot;</span><span>]</span>,m0<span>[</span><span>&quot;x2&quot;</span><span>]</span>,m1<span>[</span><span>&quot;x1&quot;</span><span>]</span>,m1<span>[</span><span>&quot;x2&quot;</span><span>]</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span><span>)</span>
<span>points</span><span>(</span>.5<span>*</span>m0<span>[</span><span>&quot;x1&quot;</span><span>]</span><span>+</span>.5<span>*</span>m1<span>[</span><span>&quot;x1&quot;</span><span>]</span>,.5<span>*</span>m0<span>[</span><span>&quot;x2&quot;</span><span>]</span><span>+</span>.5<span>*</span>m1<span>[</span><span>&quot;x2&quot;</span><span>]</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span>,pch<span>=</span><span>19</span><span>)</span></pre></td></tr></table></div>
<p>Of course, we can also use R function</p>
<div><table><tr><td><pre>1
2
3
</pre></td><td><pre>predlda <span>=</span> <span>function</span><span>(</span>x,y<span>)</span> <span>predict</span><span>(</span>fit_lda, <span>data.<span>frame</span></span><span>(</span>x1<span>=</span>x,x2<span>=</span>y<span>)</span><span>)</span>$class<span>==</span><span>1</span>
vv<span>=</span><span>outer</span><span>(</span>vu,vu,predlda<span>)</span>
<span>contour</span><span>(</span>vu,vu,vv,add<span>=</span>TRUE,lwd<span>=</span><span>2</span>,<span>levels</span> <span>=</span> .5<span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-06-a%CC%80-05.36.38.png?w=456&ssl=1" data-recalc-dims="1"><br>
One can also consider the quadratic discriminent analysis since it might be difficult to argue that <img alt="\mathbf{\Sigma}_0=\mathbf{\Sigma}_1" title="\mathbf{\Sigma}_0=\mathbf{\Sigma}_1" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5CSigma%7D_0%3D%5Cmathbf%7B%5CSigma%7D_1&bg=ffffff&fg=000&s=0" class="latex"></p>
<div><table><tr><td><pre>1
</pre></td><td><pre>fit_qda <span>=</span> qda<span>(</span>y ~x1<span>+</span>x2 , <span>data</span><span>=</span><span>df</span><span>)</span></pre></td></tr></table></div>
<p>The separation curve is here</p>
<div><table><tr><td><pre>1
2
3
4
5
</pre></td><td><pre><span>plot</span><span>(</span><span>df</span>$x1,<span>df</span>$x2,pch<span>=</span><span>19</span>,
<span>col</span><span>=</span><span>c</span><span>(</span><span>&quot;blue&quot;</span>,<span>&quot;red&quot;</span><span>)</span><span>[</span><span>1</span><span>+</span><span>(</span><span>df</span>$y<span>==</span><span>&quot;1&quot;</span><span>)</span><span>]</span><span>)</span>
predqda<span>=</span><span>function</span><span>(</span>x,y<span>)</span> <span>predict</span><span>(</span>fit_qda, <span>data.<span>frame</span></span><span>(</span>x1<span>=</span>x,x2<span>=</span>y<span>)</span><span>)</span>$class<span>==</span><span>1</span>
vv<span>=</span><span>outer</span><span>(</span>vu,vu,predlda<span>)</span>
<span>contour</span><span>(</span>vu,vu,vv,add<span>=</span>TRUE,lwd<span>=</span><span>2</span>,<span>levels</span> <span>=</span> .5<span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-06-a%CC%80-05.31.30.png?w=456&ssl=1" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://freakonometrics.hypotheses.org/53021"> R-english – Freakonometrics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/jVMyh7UXvDU" height="1"></a></div>