Skip to main content

probability or statistics - Multinomial logistic regression


Has anyone done multinomial logistic regression in Mathematica?


The binomial case is essentially done on the LogitModelFit documentation page and works fine.


I am following this formulation to obtain the outcome probabilities for k classes with k-1 calls to LogitModelFit.


I am getting normalized probabilities out, but they depend on how the classes are encoded for each independent regression. I can post what I have on request but it's a lot of code for something that isn't working.


Either multiple-binary or true multinomial would be fine for my purposes, but bonus for both implementations.



Answer



Multiple-Binary Approach


Lets bring out the famous iris data set to set this up and split it into 80% training and 20% testing data to be used in fitting the model.


iris = ExampleData[{"Statistics", "FisherIris"}];

n = Length[iris]

rs = RandomSample[Range[n]];
cut = Ceiling[.8 n];
train = iris[[rs[[1 ;; cut]]]];
test = iris[[rs[[cut + 1 ;;]]]];

speciesList = Union[iris[[All,-1]]];

First we need to run a separate regression for each species of iris. I'm creating an encoder for the response for this purpose which gives 1 if the response matches a particular species and 0 otherwise.



encodeResponse[data_, species_] := Block[{resp},
  resp = data[[All, -1]];
  Join[data[[All, 1 ;; -2]]\[Transpose], {Boole[# === species] & /@ 
      resp}]\[Transpose]
  ]

Now we fit a model for each of the three species using this encoded training data.


logmods = 
  Table[LogitModelFit[encodeResponse[train, i], Array[x, 4], 
    Array[x, 4]], {i, speciesList}];


Each model can be used to estimate a probability that a particular set of features yields its class. For classification, we simply let them all "vote" and pick the category with highest probability.


mlogitClassify[mods_, speciesList_][x_] :=
 Block[{p},
  p = #[Sequence @@ x] & /@ mods;
  speciesList[[Ordering[p, -1][[1]]]]
  ]

So how well did this perform? In this case we get 14 out of 15 correct on the testing data (pretty good!).


class = mlogitClassify[logmods, speciesList][#] & /@ test[[All, 1 ;; -2]];


Mean[Boole[Equal @@@ Thread[{class, test[[All, -1]]}]]]

(* 14/15 *)

It is interesting to see what this classifier actually does visually. To do so, I'll reduce the number of variables to 2.


logmods2 = 
  Table[LogitModelFit[encodeResponse[train[[All, {3, 4, 5}]], i], 
    Array[x, 2], Array[x, 2]], {i, speciesList}];


Show[Table[
  RegionPlot[
   mlogitClassify[logmods2, speciesList][{x1, x2}] === i, {x1, .5, 
    7}, {x2, 0, 2.5}, 
   PlotStyle -> 
    Directive[Opacity[.25], 
     Switch[i, "setosa", Red, "virginica", Green, "versicolor", 
      Blue]]], {i, {"setosa", "virginica", "versicolor"}}], 
 ListPlot[Table[
   Tooltip[Pick[iris[[All, {3, 4}]], iris[[All, -1]], i], 

    i], {i, {"setosa", "virginica", "versicolor"}}], 
  PlotStyle -> {Red, Darker[Green], Blue}]]

enter image description here


The decision boundaries are more interesting if we allow more flexibility in the basis functions. Here I allow all of the possible quadratic terms.


logmods2 = 
  Table[LogitModelFit[
    encodeResponse[train[[All, {3, 4, 5}]], i], {1, x[1], x[2], 
     x[1]^2, x[2]^2, x[1] x[2]}, Array[x, 2]], {i, speciesList}];


Show[Table[
  RegionPlot[
   mlogitClassify[logmods2, speciesList][{x1, x2}] === i, {x1, .5, 
    7}, {x2, 0, 2.5}, 
   PlotStyle -> 
    Directive[Opacity[.25], 
     Switch[i, "setosa", Red, "virginica", Green, "versicolor", 
      Blue]], PlotPoints -> 100], {i, {"setosa", "virginica", 
    "versicolor"}}], 
 ListPlot[Table[

   Tooltip[Pick[iris[[All, {3, 4}]], iris[[All, -1]], i], 
    i], {i, {"setosa", "virginica", "versicolor"}}], 
  PlotStyle -> {Red, Darker[Green], Blue}]]

enter image description here


Comments

Post a Comment

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more