[Mathematica] Байесовская классификация

[Trius]

Изучаются объекты с несколькими характеристиками. Задача состоит в том, что нужно спроектировать вектор характеристик на какое-то направление единичной длины $u$, построить одномерную ядерную оценку и по ней уже строить баесовский классификатор. Направление $u$ выбирается таким образом, чтобы минимизировать вероятность ошибки: $L=1-\int \max_{m=1,..,M}p_mf_m(x)dx$. Пишу такое:

jdata = List[data1, data2, data3];
tdata = {};
cdata = {};

p[j_] := If[j == 1, p1, 1]*If[j == 2, p2, 1]*If[j == 3, p3, 1];

f[x_, j_, u_] := ( cdata = Switch[j, 1, data1, 2, data2, 3, data3]
     For[i = 1, i <= Length[cdata], i++, 
      temp = u*cdata[[i]][[1]] + Sqrt[1 - u^2]*cdata[[i]][[2]]
         AppendTo[tdata, temp]]; 
   PDF[SmoothKernelDistribution[tdata], x]);

L[u_] := 1 - 
  Integrate[MaxValue[f[x, j, u]*p[j]], {x, -Infinity,    Infinity}]

uopt = NArgMin[{L[u], Abs[u] <= 1}, u]

data1,data2,data3 - части обучающей выборки, отвечающие соответствующей популяции. В общем, оно видимо зацикливается и выдает кучу ошибок -

(Оффтоп)

Part::partd: Part specification Null[[1]] is longer than depth of object. >>

Part::partd: Part specification Null[[2]] is longer than depth of object. >>

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

General::stop: Further output of $RecursionLimit::reclim will be suppressed during this calculation. >>

SmoothKernelDistribution::rctn: A rectangular array of real numbers is expected at position 1 in SmoothKernelDistribution[{{{{{{{{<<1>>}}}}}}},<<6>>,u Null[[1]]+{{{{{{<<1>>}}}}},{Sqrt[1+Times[<<2>>]] Null[[2]] (u Part[<<2>>]+Power[<<2>>] {<<1>>} Part[<<2>>])},Sqrt[1-Power[<<2>>]] (j u+Sqrt[Plus[<<2>>]] {{<<1>>},<<1>>}) Null[[2]],{<<1>>},<<1>>,Sqrt[1-Power[<<2>>]] Null[[2]] <<1>>,Sqrt[1-Power[<<2>>]] Null[[2]] (j u+{{{<<1>>}},{Times[<<2>>]},Power[<<2>>] Plus[<<2>>],{{<<1>>},Times[<<2>>],Times[<<2>>]},Power[<<2>>] Plus[<<2>>],Power[<<2>>] Plus[<<2>>]})}}]. >>

Подскажите пожалуйста, что не так

[Leierkastenmann]

Trius, как минимум надо бы функцию f правильно записать. Если функция представляет собой набор действий, то ее надо оформить при помощи инструкции Block или Module. В функции L(u) где задается j? Это что касается частностей.
А в целом, боюсь, что такая конструкция может не заработать, если математика не сможет взять интеграл. Я бы заменил это на численное интегрирование, хотя не факт что это поможет.

[Trius]

Исправил немного

f[x1_, j1_, u1_] := 
  Module[ {x = x1, j = j1, u = u1}, 
   cdata := Switch[j, 1, data1, 2, data2, 3, data3]
     For[i = 1, i <= Length[cdata], i++, 
      temp := u*cdata[[i]][[1]] + Sqrt[1 - u^2]*cdata[[i]][[2]]
         AppendTo[tdata, temp]]; 
   PDF[SmoothKernelDistribution[tdata], x]];

L[u_] := 1 - 
  NIntegrate[
   MaxValue[f[x, j, u]*p[j], j \[Element] {1, 2, 3}], {x, -Infinity, 
    Infinity}]

uopt = NArgMin[{L[u], Abs[u] <= 1}, u]

Сразу же выдает

SmoothKernelDistribution::rctn: A rectangular array of real numbers is expected at position 1 in SmoothKernelDistribution[{{u Null[[1]]+{{{{{<<1>>}}}}}},<<11>>}]. >>

и дальше что-то якобы считает. Ошибка видимо еще в заполнении tdata
Можно конечно сделать как одногруппник, вручную построить ядерную оценку, но мне не особо это нравиться.
Кстати, можно ли как-то вытащить из SmoothKernelDistribution полученный параметр сглаживания?

[Leierkastenmann]

В L(u) по-прежнему есть неточности - не указана переменная, по отношению к которой ищется максимум. И приведите еще значения для data, чтобы можно было полностью воспроизвести ваш код.

-- Вт май 14, 2013 21:23:44 --

Кстати, и почему вы используете MaxValue, которая ищет максимум функции, когда у вас в формуле вероятности ошибки просто максимум для набора функций? Это же разные вещи, мне кажется вместо MaxValue надо использовать Max. Например вот так:

L[u_] := 1 - 
  NIntegrate[
   Max[f[x, #, u]*p[#] & /@ {1, 2, 3}], {x, -Infinity, Infinity}]

[Trius]

Указано,

j \[Element] {1, 2, 3}

Данные

data1 = {};
data2 = {};
data3 = {};

In[49]:= GoodData

Out[49]= {{1, 14.23, 2.8}, {1, 13.2, 2.65}, {1, 13.16, 2.8}, {1, 
  14.37, 3.85}, {1, 13.24, 2.8}, {1, 14.2, 3.27}, {1, 14.39, 2.5}, {1,
   14.06, 2.6}, {1, 14.83, 2.8}, {1, 13.86, 2.98}, {1, 14.1, 
  2.95}, {1, 14.12, 2.2}, {1, 13.75, 2.6}, {1, 14.75, 3.1}, {1, 14.38,
   3.3}, {1, 13.63, 2.85}, {1, 14.3, 2.8}, {1, 13.83, 2.95}, {1, 
  14.19, 3.3}, {1, 13.64, 2.7}, {1, 14.06, 3.}, {1, 12.93, 2.41}, {1, 
  13.71, 2.61}, {1, 12.85, 2.48}, {1, 13.5, 2.53}, {1, 13.05, 
  2.63}, {1, 13.39, 2.85}, {1, 13.3, 2.4}, {1, 13.87, 2.95}, {1, 
  14.02, 2.65}, {1, 13.73, 3.}, {1, 13.58, 2.86}, {1, 13.68, 
  2.42}, {1, 13.76, 2.95}, {1, 13.51, 2.35}, {1, 13.48, 2.7}, {1, 
  13.28, 2.6}, {1, 13.05, 2.45}, {1, 13.07, 2.4}, {1, 14.22, 3.}, {1, 
  13.56, 3.15}, {1, 13.41, 2.45}, {1, 13.88, 3.25}, {1, 13.24, 
  2.64}, {1, 13.05, 3.}, {1, 14.21, 2.85}, {1, 14.38, 3.25}, {1, 13.9,
   3.1}, {1, 14.1, 2.75}, {1, 13.94, 2.88}, {1, 13.05, 2.72}, {1, 
  13.83, 2.45}, {1, 13.82, 3.88}, {1, 13.77, 3.}, {1, 13.74, 2.6}, {1,
   13.56, 2.96}, {1, 14.22, 3.2}, {1, 13.29, 3.}, {1, 13.72, 3.4}, {2,
   12.37, 1.98}, {2, 12.33, 2.05}, {2, 12.64, 2.02}, {2, 13.67, 
  2.1}, {2, 12.37, 3.5}, {2, 12.17, 1.89}, {2, 12.37, 2.42}, {2, 
  13.11, 2.98}, {2, 12.37, 2.11}, {2, 13.34, 2.53}, {2, 12.21, 
  1.85}, {2, 12.29, 1.1}, {2, 13.86, 2.95}, {2, 13.49, 1.88}, {2, 
  12.99, 3.3}, {2, 11.96, 3.38}, {2, 11.66, 1.61}, {2, 13.03, 
  1.95}, {2, 11.84, 1.72}, {2, 12.33, 1.9}, {2, 12.7, 2.83}, {2, 12., 
  2.42}, {2, 12.72, 2.2}, {2, 12.08, 2.}, {2, 13.05, 1.65}, {2, 11.84,
   2.2}, {2, 12.67, 2.2}, {2, 12.16, 1.78}, {2, 11.65, 1.92}, {2, 
  11.64, 1.95}, {2, 12.08, 2.2}, {2, 12.08, 1.6}, {2, 12., 1.45}, {2, 
  12.69, 1.38}, {2, 12.29, 2.45}, {2, 11.62, 3.02}, {2, 12.47, 
  2.5}, {2, 11.81, 1.6}, {2, 12.29, 2.55}, {2, 12.37, 3.52}, {2, 
  12.29, 2.85}, {2, 12.08, 2.23}, {2, 12.6, 1.45}, {2, 12.34, 
  2.56}, {2, 11.82, 2.5}, {2, 12.51, 2.2}, {2, 12.42, 1.68}, {2, 
  12.25, 1.65}, {2, 12.72, 1.38}, {2, 12.22, 2.36}, {2, 11.61, 
  2.74}, {2, 11.46, 3.18}, {2, 12.52, 2.55}, {2, 11.76, 1.75}, {2, 
  11.41, 2.48}, {2, 12.08, 2.56}, {2, 11.03, 2.46}, {2, 11.82, 
  1.98}, {2, 12.42, 2.}, {2, 12.77, 1.63}, {2, 12., 2.}, {2, 11.45, 
  2.9}, {2, 11.56, 3.18}, {2, 12.42, 2.2}, {2, 13.05, 2.62}, {2, 
  11.87, 2.86}, {2, 12.07, 2.6}, {2, 12.43, 2.74}, {2, 11.79, 
  2.13}, {2, 12.37, 2.22}, {2, 12.04, 2.1}, {3, 12.86, 1.51}, {3, 
  12.88, 1.3}, {3, 12.81, 1.15}, {3, 12.7, 1.7}, {3, 12.51, 2.}, {3, 
  12.6, 1.62}, {3, 12.25, 1.38}, {3, 12.53, 1.79}, {3, 13.49, 
  1.62}, {3, 12.84, 2.32}, {3, 12.93, 1.54}, {3, 13.36, 1.4}, {3, 
  13.52, 1.55}, {3, 13.62, 2.}, {3, 12.25, 1.38}, {3, 13.16, 1.5}, {3,
   13.88, 0.98}, {3, 12.87, 1.7}, {3, 13.32, 1.93}, {3, 13.08, 
  1.41}, {3, 13.5, 1.4}, {3, 12.79, 1.48}, {3, 13.11, 2.2}, {3, 13.23,
   1.8}, {3, 12.58, 1.48}, {3, 13.17, 1.74}, {3, 13.84, 1.8}, {3, 
  12.45, 1.9}, {3, 14.34, 2.8}, {3, 13.48, 2.6}, {3, 12.36, 2.3}, {3, 
  13.69, 1.83}, {3, 12.85, 1.65}, {3, 12.96, 1.39}, {3, 13.78, 
  1.35}, {3, 13.73, 1.28}, {3, 13.45, 1.7}, {3, 12.82, 1.48}, {3, 
  13.58, 1.55}, {3, 13.4, 1.98}, {3, 12.2, 1.25}, {3, 12.77, 
  1.39}, {3, 14.16, 1.68}, {3, 13.71, 1.68}, {3, 13.4, 1.8}, {3, 
  13.27, 1.59}, {3, 13.17, 1.65}}

In[14]:= For[i = 1, i <= Length[GoodData], i++,
     obj = Join[{GoodData[[i]][[2]]}, {GoodData[[i]][[3]]}];
    If[GoodData[[i]][[1]] == 1, AppendTo[data1, obj]];
    If[GoodData[[i]][[1]] == 2, AppendTo[data2, obj]];           
  If[GoodData[[i]][[1]] == 3, AppendTo[data3, obj]]];

[Leierkastenmann]

Саму функцию L вычислить можно так:

f[x_, j_, u_, data_] := Module[
   {
    cdata,
    tdata
    },
   cdata = Rest[#] & /@ Select[data, First[#] == j &];
   tdata = (u*First[#] + Sqrt[1 - u^2]*Last[#]) & /@ cdata;
   Return[PDF[SmoothKernelDistribution[tdata], x]]
   ];
p = {p1, p2, p3};
L[u_] := 1 - 
   NIntegrate[
    Max[f[x, #, u, GoodData]*p[[#]] & /@ {1, 2, 3}], {x, -Infinity, 
     Infinity}];

Не забудьте только, что все p должны быть заданы численно. Правда эта функция очень долго считает и найти ее минимум стандартными способами не удастся. Процедуру для нахождения минимума придется написать самому, но работать она будет очень и очень долго.

[Trius]

С данными, которые в SmoothKernel передаются, не то что-то. Или просто нельзя так делать.

In[39]:= 
uopt = NArgMin[{L[u_], Abs[u] <= 1}, u]

During evaluation of In[39]:= SmoothKernelDistribution::rctn: A rectangular array of real numbers is expected at position 1 in SmoothKernelDistribution[{14.23 u_+2.8 Sqrt[1-Pattern[<<2>>]^2],13.2 u_+2.65 Sqrt[1-Pattern[<<2>>]^2],13.16 u_+2.8 Sqrt[1-Pattern[<<2>>]^2],14.37 u_+3.85 Sqrt[1-Pattern[<<2>>]^2],13.24 u_+2.8 Sqrt[1-Pattern[<<2>>]^2],<<41>>,14.38 u_+3.25 Sqrt[1-Pattern[<<2>>]^2],13.9 u_+3.1 Sqrt[1-Pattern[<<2>>]^2],14.1 u_+2.75 Sqrt[1-Pattern[<<2>>]^2],13.94 u_+2.88 Sqrt[1-Pattern[<<2>>]^2],<<9>>}]. >>

During evaluation of In[39]:= SmoothKernelDistribution::rctn: A rectangular array of real numbers is expected at position 1 in SmoothKernelDistribution[{12.37 u_+1.98 Sqrt[1-Pattern[<<2>>]^2],12.33 u_+2.05 Sqrt[1-Pattern[<<2>>]^2],12.64 u_+2.02 Sqrt[1-Pattern[<<2>>]^2],13.67 u_+2.1 Sqrt[1-Pattern[<<2>>]^2],12.37 u_+3.5 Sqrt[1-Pattern[<<2>>]^2],<<41>>,12.42 u_+1.68 Sqrt[1-Pattern[<<2>>]^2],12.25 u_+1.65 Sqrt[1-Pattern[<<2>>]^2],12.72 u_+1.38 Sqrt[1-Pattern[<<2>>]^2],12.22 u_+2.36 Sqrt[1-Pattern[<<2>>]^2],<<21>>}]. >>

During evaluation of In[39]:= SmoothKernelDistribution::rctn: A rectangular array of real numbers is expected at position 1 in SmoothKernelDistribution[{12.86 u_+1.51 Sqrt[1-Pattern[<<2>>]^2],12.88 u_+1.3 Sqrt[1-Pattern[<<2>>]^2],12.81 u_+1.15 Sqrt[1-Pattern[<<2>>]^2],12.7 u_+1.7 Sqrt[1-Pattern[<<2>>]^2],12.51 u_+2. Sqrt[1-Pattern[<<2>>]^2],<<37>>,14.16 u_+1.68 Sqrt[1-Pattern[<<2>>]^2],13.71 u_+1.68 Sqrt[1-Pattern[<<2>>]^2],13.4 u_+1.8 Sqrt[1-Pattern[<<2>>]^2],13.27 u_+1.59 Sqrt[1-Pattern[<<2>>]^2],13.17 u_+1.65 Sqrt[1-Pattern[<<2>>]^2]}]. >>

During evaluation of In[39]:= General::stop: Further output of SmoothKernelDistribution::rctn will be suppressed during this calculation. >>

During evaluation of In[39]:= NIntegrate::inumr: The integrand Max[47/177 PDF[SmoothKernelDistribution[{12.86 Pattern[<<2>>]+1.51 Power[<<2>>],12.88 Pattern[<<2>>]+1.3 Power[<<2>>],12.81 Pattern[<<2>>]+1.15 Power[<<2>>],12.7 Pattern[<<2>>]+1.7 Power[<<2>>],<<39>>,13.71 Pattern[<<2>>]+1.68 Power[<<2>>],13.4 Pattern[<<2>>]+1.8 Power[<<2>>],13.27 Pattern[<<2>>]+1.59 Power[<<2>>],13.17 Pattern[<<2>>]+1.65 Power[<<2>>]}],x],1/3 PDF[<<1>>,x],71/177 <<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-\[Infinity],0.}}. >>

During evaluation of In[39]:= NIntegrate::inumr: The integrand Max[47/177 PDF[SmoothKernelDistribution[{12.86 Pattern[<<2>>]+1.51 Power[<<2>>],12.88 Pattern[<<2>>]+1.3 Power[<<2>>],12.81 Pattern[<<2>>]+1.15 Power[<<2>>],12.7 Pattern[<<2>>]+1.7 Power[<<2>>],<<39>>,13.71 Pattern[<<2>>]+1.68 Power[<<2>>],13.4 Pattern[<<2>>]+1.8 Power[<<2>>],13.27 Pattern[<<2>>]+1.59 Power[<<2>>],13.17 Pattern[<<2>>]+1.65 Power[<<2>>]}],x],1/3 PDF[<<1>>,x],71/177 <<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-\[Infinity],0.}}. >>

During evaluation of In[39]:= NIntegrate::inumr: The integrand Max[47/177 PDF[SmoothKernelDistribution[{12.86 Pattern[<<2>>]+1.51 Power[<<2>>],12.88 Pattern[<<2>>]+1.3 Power[<<2>>],12.81 Pattern[<<2>>]+1.15 Power[<<2>>],12.7 Pattern[<<2>>]+1.7 Power[<<2>>],<<39>>,13.71 Pattern[<<2>>]+1.68 Power[<<2>>],13.4 Pattern[<<2>>]+1.8 Power[<<2>>],13.27 Pattern[<<2>>]+1.59 Power[<<2>>],13.17 Pattern[<<2>>]+1.65 Power[<<2>>]}],x],1/3 PDF[<<1>>,x],71/177 <<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-\[Infinity],0.}}. >>

During evaluation of In[39]:= General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation. >>

During evaluation of In[39]:= NIntegrate::nlim: x = -1. Compile`$144 is not a valid limit of integration. >>

During evaluation of In[39]:= NArgMin::nnum: The function value 1-NIntegrate[Max[(f[x,Slot[<<1>>],Pattern[<<2>>],{<<177>>}] {<<3>>}[[Slot[<<1>>]]]&)/@{1,2,3}],{x,-\[Infinity],\[Infinity]}] is not a number at {u} = {-0.829053}. >>

During evaluation of In[39]:= NIntegrate::nlim: x = -1. Compile`$144 is not a valid limit of integration. >>

During evaluation of In[39]:= NArgMin::nnum: The function value 1-NIntegrate[Max[(f[x,Slot[<<1>>],Pattern[<<2>>],{<<177>>}] {<<3>>}[[Slot[<<1>>]]]&)/@{1,2,3}],{x,-\[Infinity],\[Infinity]}] is not a number at {u} = {-0.829053}. >>

Out[39]= NArgMin[{1 - 
   NIntegrate[
    Max[(f[x, #1, u_, GoodData] p[[#1]] &) /@ {1, 2, 
       3}], {x, -\[Infinity], \[Infinity]}], Abs[u] <= 1}, u]

Печалька
Повторю вопрос: можно ли как-то вытащить из SmoothKernelDistribution полученный параметр сглаживания?

[Leierkastenmann]

Trius, я же и написал, что стандартные методы минимизации для этой функции работать не будут, поэтому и пытаться даже не стоило искать минимум этим способом. Надо писать свою минимизирующую процедуру.
Насчет SmoothKernelDistribution, судя по тому, что я прочитал о нем в хэлпе, ничего из него вытащить нельзя.

[Trius]

Ладно, я сделал топором Посчитал $L(u)$ для $-1\le u<1$ с шагом 0.01. А тут уже минимум найти нечего делать. Спасибо за помощь