Results for paper about cubic LR

python/logistic/cost_checks.nb

@@ -0,0 +1,588 @@
+(* Content-type: application/vnd.wolfram.mathematica *)
+
+(*** Wolfram Notebook File ***)
+(* http://www.wolfram.com/nb *)
+
+(* CreatedBy='Wolfram 14.3' *)
+
+(*CacheID: 234*)
+(* Internal cache information:
+NotebookFileLineBreakTest
+NotebookFileLineBreakTest
+NotebookDataPosition[       154,          7]
+NotebookDataLength[     24335,        580]
+NotebookOptionsPosition[     22033,        534]
+NotebookOutlinePosition[     22429,        550]
+CellTagsIndexPosition[     22386,        547]
+WindowFrame->Normal*)
+
+(* Beginning of Notebook Content *)
+Notebook[{
+Cell[BoxData[{
+ RowBox[{
+  RowBox[{"Clear", "[", 
+   RowBox[{"GetBeta", ",", " ", "Cost"}], "]"}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{
+    RowBox[{"GetBeta", "[", 
+     RowBox[{"{", 
+      RowBox[{"cc_", ",", "epsilon_", ",", "k1_", ",", "k2_"}], "}"}], "]"}], 
+    ":=", 
+    RowBox[{"Module", "[", 
+     RowBox[{
+      RowBox[{"{", 
+       RowBox[{"beta", ",", "a", ",", "b", ",", "c", ",", "d"}], "}"}], ",", 
+      "\[IndentingNewLine]", 
+      RowBox[{
+       RowBox[{"d", "=", "cc"}], ";", "\[IndentingNewLine]", 
+       RowBox[{"c", "=", 
+        RowBox[{
+         RowBox[{"k2", "^", "2"}], "+", 
+         RowBox[{"epsilon", "^", "2"}]}]}], ";", "\[IndentingNewLine]", 
+       RowBox[{"b", "=", 
+        RowBox[{"k1", "*", "k2"}]}], ";", "\[IndentingNewLine]", 
+       RowBox[{"a", "=", 
+        RowBox[{
+         RowBox[{"k1", "^", "2"}], "/", "3"}]}], ";", "\[IndentingNewLine]", 
+       "\[IndentingNewLine]", 
+       RowBox[{"Return", "[", 
+        RowBox[{"{", 
+         RowBox[{"d", ",", "c", ",", "b", ",", "a"}], "}"}], "]"}]}]}], 
+     "\[IndentingNewLine]", "]"}]}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"Cost", "[", 
+   RowBox[{"pars_", ",", "lam2_", ",", "x_", ",", "y_"}], "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"beta", ",", "X", ",", "nllf", ",", "s", ",", "i"}], "}"}], ",", 
+    "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"beta", " ", "=", " ", 
+      RowBox[{"GetBeta", "[", "pars", "]"}]}], ";", "\[IndentingNewLine]", 
+     "\[IndentingNewLine]", 
+     RowBox[{"s", "=", 
+      RowBox[{
+       RowBox[{"2", "*", "y"}], "-", "1"}]}], ";", "\[IndentingNewLine]", 
+     RowBox[{"X", "=", 
+      RowBox[{"Transpose", "[", 
+       RowBox[{"Table", "[", 
+        RowBox[{
+         RowBox[{"x", "^", "i"}], ",", 
+         RowBox[{"{", 
+          RowBox[{"i", ",", "0", ",", "3"}], "}"}]}], "]"}], "]"}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"nllf", " ", "=", " ", 
+      RowBox[{"Total", "[", 
+       RowBox[{"Log", "[", 
+        RowBox[{"1", "+", 
+         RowBox[{"Exp", "[", 
+          RowBox[{
+           RowBox[{"-", "s"}], "*", 
+           RowBox[{"(", 
+            RowBox[{"X", ".", "beta"}], ")"}]}], "]"}]}], "]"}], "]"}]}], ";",
+      "\[IndentingNewLine]", "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", 
+      RowBox[{"nllf", " ", "+", " ", 
+       RowBox[{"lam2", " ", 
+        RowBox[{"Total", "[", 
+         RowBox[{"pars", "^", "2"}], "]"}]}]}], "]"}]}]}], 
+   "\[IndentingNewLine]", "]"}]}]}], "Input",
+ CellChangeTimes->{{3.980215196159321*^9, 3.98021520112125*^9}, {
+  3.980215465357512*^9, 3.9802155187192497`*^9}, {3.980215561493145*^9, 
+  3.98021556877605*^9}, {3.980215603281454*^9, 3.980215766491699*^9}, {
+  3.980215809381616*^9, 3.980215824734212*^9}, {3.9802159883378687`*^9, 
+  3.980216092393737*^9}, {3.9802179711796503`*^9, 3.980217994933426*^9}, {
+  3.980218037526173*^9, 3.9802180508589983`*^9}, {3.980218114391039*^9, 
+  3.9802182123715487`*^9}, {3.980218668097128*^9, 3.9802186800955877`*^9}, {
+  3.980219900715466*^9, 3.980219900917904*^9}, {3.9802199440804777`*^9, 
+  3.980219950481944*^9}, {3.980220901657138*^9, 
+  3.9802209090949917`*^9}},ExpressionUUID->"db312a24-878f-40c8-962f-\
+6aedbdf27fe4"],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"(*", 
+   RowBox[{"TRIM", ",", " ", "log"}], "*)"}], "\[IndentingNewLine]", 
+  RowBox[{
+   RowBox[{"x", "=", 
+    RowBox[{"{", 
+     RowBox[{
+     "0.7444588568509963", ",", "0.4294306635860587", ",", 
+      "0.2424791244329582", ",", "0.1722166015253219", ",", 
+      "0.5155957608380367", ",", 
+      RowBox[{"-", "0.0044000248522457"}], ",", "0.2619898405306257", ",", 
+      "0.3707014579741952", ",", "0.3458779082346233", ",", 
+      RowBox[{"-", "0.0230997734892581"}], ",", "1.4235506263993014", ",", 
+      "0.1974151389300383", ",", "0.3479428293749495", ",", 
+      "0.038497594143739", ",", 
+      RowBox[{"-", "0.0494408800705622"}], ",", "0.1562561841152657", ",", 
+      "0.2440810426279506", ",", "0.738605475021583", ",", 
+      "0.3388466870335315", ",", "0.2730058878681247", ",", 
+      "0.2349466740981531", ",", "0.0152345391191637", ",", 
+      RowBox[{"-", "0.1306983645886834"}], ",", "0.3863126943346241", ",", 
+      "0.32696785342031", ",", "0.3960496407131636", ",", "0.116543550839155",
+       ",", "0.3745678322136136", ",", "0.3902449371047537", ",", 
+      RowBox[{"-", "0.0224728544974076"}], ",", "0.3120049574170508", ",", 
+      "0.140909818171296", ",", "0.4065884615650837", ",", 
+      "0.1347990203891564", ",", "0.1747544412028135", ",", 
+      "0.8025435019583259", ",", "0.3424702625117166", ",", 
+      RowBox[{"-", "0.2576749555237702"}], ",", "0.2098741187455145", ",", 
+      "0.1827383448113199", ",", "0.4344451063845957", ",", 
+      "0.3628263770555691", ",", "0.3418422883292066", ",", 
+      "0.2690071065545839", ",", "0.3785734614926337", ",", 
+      "0.2619248908950089", ",", "0.232085349459071", ",", 
+      "0.1599028123673964", ",", "0.5800908847567773", ",", 
+      "0.208084169807583", ",", 
+      RowBox[{"-", "0.0239525163747711"}], ",", "0.2010508033108026", ",", 
+      "0.4280507666730073", ",", "0.4184257480672208", ",", 
+      "0.444582693795233", ",", "0.2107665896841349", ",", 
+      "0.3642384221371574"}], "}"}]}], "\[IndentingNewLine]", 
+   RowBox[{"y", "=", 
+    RowBox[{"{", 
+     RowBox[{
+     "1", ",", "0", ",", "0", ",", "0", ",", "1", ",", "0", ",", "0", ",", 
+      "0", ",", "0", ",", "0", ",", "1", ",", "0", ",", "0", ",", "0", ",", 
+      "0", ",", "0", ",", "0", ",", "1", ",", "0", ",", "0", ",", "0", ",", 
+      "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", 
+      "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", 
+      "1", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", 
+      "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", 
+      "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", 
+      "0"}], "}"}]}], "\[IndentingNewLine]", 
+   RowBox[{"p", "=", 
+    RowBox[{"{", 
+     RowBox[{
+      RowBox[{
+       RowBox[{"-", "1.1319747348643707"}], " ", 
+       RowBox[{"10", "^", "1"}]}], ",", 
+      RowBox[{
+       RowBox[{"-", "6.5316161792581431"}], "*", " ", 
+       RowBox[{"10", "^", 
+        RowBox[{"-", "5"}]}]}], ",", "1.6944026199144295", ",", 
+      "3.9357124963564112"}], "}"}]}], "\[IndentingNewLine]", 
+   RowBox[{"lam2", "=", 
+    RowBox[{"10", "^", 
+     RowBox[{"-", "3"}]}]}]}]}]], "Input",
+ CellChangeTimes->{{3.980219550729746*^9, 3.980219552713852*^9}, {
+  3.980219691095812*^9, 3.9802197156948*^9}, {3.98022297307031*^9, 
+  3.980222979562756*^9}},ExpressionUUID->"f22ff9ef-193c-434b-91da-\
+e9e8b3175259"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+  "0.7444588568509963`", ",", "0.4294306635860587`", ",", 
+   "0.2424791244329582`", ",", "0.1722166015253219`", ",", 
+   "0.5155957608380367`", ",", 
+   RowBox[{"-", "0.0044000248522457`"}], ",", "0.2619898405306257`", ",", 
+   "0.3707014579741952`", ",", "0.3458779082346233`", ",", 
+   RowBox[{"-", "0.0230997734892581`"}], ",", "1.4235506263993014`", ",", 
+   "0.1974151389300383`", ",", "0.3479428293749495`", ",", 
+   "0.038497594143739`", ",", 
+   RowBox[{"-", "0.0494408800705622`"}], ",", "0.1562561841152657`", ",", 
+   "0.2440810426279506`", ",", "0.738605475021583`", ",", 
+   "0.3388466870335315`", ",", "0.2730058878681247`", ",", 
+   "0.2349466740981531`", ",", "0.0152345391191637`", ",", 
+   RowBox[{"-", "0.1306983645886834`"}], ",", "0.3863126943346241`", ",", 
+   "0.32696785342031`", ",", "0.3960496407131636`", ",", "0.116543550839155`",
+    ",", "0.3745678322136136`", ",", "0.3902449371047537`", ",", 
+   RowBox[{"-", "0.0224728544974076`"}], ",", "0.3120049574170508`", ",", 
+   "0.140909818171296`", ",", "0.4065884615650837`", ",", 
+   "0.1347990203891564`", ",", "0.1747544412028135`", ",", 
+   "0.8025435019583259`", ",", "0.3424702625117166`", ",", 
+   RowBox[{"-", "0.2576749555237702`"}], ",", "0.2098741187455145`", ",", 
+   "0.1827383448113199`", ",", "0.4344451063845957`", ",", 
+   "0.3628263770555691`", ",", "0.3418422883292066`", ",", 
+   "0.2690071065545839`", ",", "0.3785734614926337`", ",", 
+   "0.2619248908950089`", ",", "0.232085349459071`", ",", 
+   "0.1599028123673964`", ",", "0.5800908847567773`", ",", 
+   "0.208084169807583`", ",", 
+   RowBox[{"-", "0.0239525163747711`"}], ",", "0.2010508033108026`", ",", 
+   "0.4280507666730073`", ",", "0.4184257480672208`", ",", 
+   "0.444582693795233`", ",", "0.2107665896841349`", ",", 
+   "0.3642384221371574`"}], "}"}]], "Output",
+ CellChangeTimes->{3.9802197180368643`*^9, 3.9802229692825193`*^9},
+ CellLabel->"Out[59]=",ExpressionUUID->"a7f74a5e-642c-4b21-9dc6-dc4532fcd2dd"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+  "1", ",", "0", ",", "0", ",", "0", ",", "1", ",", "0", ",", "0", ",", "0", ",",
+    "0", ",", "0", ",", "1", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0",
+    ",", "0", ",", "1", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",",
+    "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0",
+    ",", "0", ",", "0", ",", "0", ",", "0", ",", "1", ",", "0", ",", "0", ",",
+    "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0",
+    ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",", "0", ",",
+    "0", ",", "0", ",", "0", ",", "0"}], "}"}]], "Output",
+ CellChangeTimes->{3.9802197180368643`*^9, 3.980222969285266*^9},
+ CellLabel->"Out[60]=",ExpressionUUID->"6fe63fdd-9529-4ff3-95a7-2b315610cc0a"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+   RowBox[{"-", "11.319747348643707`"}], ",", 
+   RowBox[{"-", "0.00006531616179258143`"}], ",", "1.6944026199144295`", ",", 
+   "3.935712496356411`"}], "}"}]], "Output",
+ CellChangeTimes->{3.9802197180368643`*^9, 3.98022296928763*^9},
+ CellLabel->"Out[61]=",ExpressionUUID->"9174505c-bb9d-4cfc-a204-c6c464fab952"],
+
+Cell[BoxData[
+ FractionBox["1", "1000"]], "Output",
+ CellChangeTimes->{3.9802197180368643`*^9, 3.980222969288666*^9},
+ CellLabel->"Out[62]=",ExpressionUUID->"dfb9c680-45f0-495e-bc68-7c102a014c4a"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"Cost", "[", 
+  RowBox[{"p", ",", "lam2", ",", "x", ",", "y"}], "]"}]], "Input",
+ CellChangeTimes->{{3.9802197403954973`*^9, 3.980219744640079*^9}},
+ CellLabel->"In[17]:=",ExpressionUUID->"73502887-6057-4d13-881a-8932ae266804"],
+
+Cell[BoxData["2.9061537536536557`"], "Output",
+ CellChangeTimes->{{3.980219744932208*^9, 3.9802197515000153`*^9}},
+ CellLabel->"Out[17]=",ExpressionUUID->"603f6210-ea83-425d-9d61-d60dffe5c701"]
+}, Open  ]],
+
+Cell[BoxData[
+ RowBox[{"Needs", "[", "\"\<NumericalCalculus`\>\"", "]"}]], "Input",
+ CellChangeTimes->{{3.9802198372526817`*^9, 3.980219843824862*^9}},
+ CellLabel->"In[18]:=",ExpressionUUID->"417b4445-9108-40f3-a80b-fda472fd9d54"],
+
+Cell[CellGroupData[{
+
+Cell[BoxData["p"], "Input",
+ CellChangeTimes->{3.9802203983793707`*^9},
+ CellLabel->"In[38]:=",ExpressionUUID->"3ca07280-6eb4-41d4-8ab2-dfe9bbf7b8be"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+   RowBox[{"-", "11.319747348643707`"}], ",", 
+   RowBox[{"-", "0.00006531616179258143`"}], ",", "1.6944026199144295`", ",", 
+   "3.935712496356411`"}], "}"}]], "Output",
+ CellChangeTimes->{3.980220398702986*^9},
+ CellLabel->"Out[38]=",ExpressionUUID->"029d324b-6436-44a3-8803-14c8381b5896"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"(*", "gradient", "*)"}], "\[IndentingNewLine]", 
+  RowBox[{
+   RowBox[{
+    RowBox[{
+     RowBox[{"DCost", "[", 
+      RowBox[{"i_", ",", "p_"}], "]"}], ":=", 
+     RowBox[{"Module", "[", 
+      RowBox[{
+       RowBox[{"{", 
+        RowBox[{"u", ",", "ii"}], "}"}], ",", 
+       RowBox[{"ND", "[", 
+        RowBox[{
+         RowBox[{"Cost", "[", 
+          RowBox[{
+           RowBox[{"Table", "[", 
+            RowBox[{
+             RowBox[{"If", "[", 
+              RowBox[{
+               RowBox[{"i", "==", "ii"}], ",", "u", ",", 
+               RowBox[{"p", "[", 
+                RowBox[{"[", "ii", "]"}], "]"}]}], "]"}], ",", 
+             RowBox[{"{", 
+              RowBox[{"ii", ",", 
+               RowBox[{"Length", "[", "p", "]"}]}], "}"}]}], "]"}], ",", 
+           "lam2", ",", "x", ",", "y"}], "]"}], ",", "u", ",", 
+         RowBox[{"p", "[", 
+          RowBox[{"[", "i", "]"}], "]"}], ",", 
+         RowBox[{"Terms", "->", "10"}]}], "]"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", 
+   RowBox[{"Table", "[", 
+    RowBox[{
+     RowBox[{"DCost", "[", 
+      RowBox[{"i", ",", "p"}], "]"}], ",", 
+     RowBox[{"{", 
+      RowBox[{"i", ",", 
+       RowBox[{"Length", "[", "p", "]"}]}], "}"}]}], "]"}]}]}]], "Input",
+ CellChangeTimes->{{3.98021979703826*^9, 3.980219827377369*^9}, {
+  3.980219863431608*^9, 3.9802198819049683`*^9}, {3.980219995402974*^9, 
+  3.980219997759408*^9}, {3.980220196481331*^9, 3.980220213434679*^9}, {
+  3.980220348281417*^9, 3.980220363235091*^9}, {3.980220425600909*^9, 
+  3.9802205459722*^9}, {3.980220589175432*^9, 3.980220594692799*^9}, {
+  3.9802207059692497`*^9, 3.9802207262248*^9}, {3.9802210428429956`*^9, 
+  3.980221072838937*^9}, {3.980221280858418*^9, 3.980221289567711*^9}, {
+  3.980221422741997*^9, 
+  3.98022142313131*^9}},ExpressionUUID->"b0ded771-dccd-45f6-a46c-\
+9b2ff466801f"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+   RowBox[{"-", "2.02166414330068`*^-6"}], ",", 
+   RowBox[{"-", "3.137769860503177`*^-8"}], ",", "4.938462760785317`*^-7", ",", 
+   RowBox[{"-", "2.546636364252448`*^-6"}]}], "}"}]], "Output",
+ CellChangeTimes->{{3.980219870760887*^9, 3.9802199094790907`*^9}, 
+   3.980219962691937*^9, {3.980219998174837*^9, 3.9802200191940937`*^9}, {
+   3.980220197643147*^9, 3.980220214115034*^9}, 3.980220363622348*^9, {
+   3.98022042632763*^9, 3.9802204383576603`*^9}, 3.980220547299224*^9, {
+   3.980220712553927*^9, 3.9802207272890863`*^9}, {3.980221066313427*^9, 
+   3.980221073167602*^9}, 3.980221290310155*^9},
+ CellLabel->"Out[49]=",ExpressionUUID->"2d1c80cf-dfd3-4a34-8f1d-5ca8d4f13801"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"(*", " ", "Hessian", "*)"}], "\[IndentingNewLine]", 
+  RowBox[{
+   RowBox[{
+    RowBox[{
+     RowBox[{"HCost", "[", 
+      RowBox[{"i_", ",", "j_", ",", "p"}], "]"}], ":=", 
+     RowBox[{"Module", "[", 
+      RowBox[{
+       RowBox[{"{", 
+        RowBox[{"u", ",", "v", ",", "k", ",", "F", ",", "G"}], "}"}], ",", 
+       "\[IndentingNewLine]", 
+       RowBox[{
+        RowBox[{
+         RowBox[{"F", "[", "v_", "]"}], ":=", 
+         RowBox[{"ND", "[", 
+          RowBox[{
+           RowBox[{"Cost", "[", 
+            RowBox[{
+             RowBox[{"Table", "[", 
+              RowBox[{
+               RowBox[{"If", "[", 
+                RowBox[{
+                 RowBox[{"i", "==", "k"}], ",", "u", ",", 
+                 RowBox[{"If", "[", 
+                  RowBox[{
+                   RowBox[{"j", "==", "k"}], ",", "v", ",", 
+                   RowBox[{"p", "[", 
+                    RowBox[{"[", "k", "]"}], "]"}]}], "]"}]}], "]"}], ",", 
+               RowBox[{"{", 
+                RowBox[{"k", ",", 
+                 RowBox[{"Length", "[", "p", "]"}]}], "}"}]}], "]"}], ",", 
+             "lam2", ",", "x", ",", "y"}], "]"}], ",", "u", ",", 
+           RowBox[{"p", "[", 
+            RowBox[{"[", "i", "]"}], "]"}], ",", 
+           RowBox[{"Terms", "->", "10"}]}], "]"}]}], ";", 
+        "\[IndentingNewLine]", 
+        RowBox[{"G", ":=", 
+         RowBox[{"ND", "[", 
+          RowBox[{
+           RowBox[{"Cost", "[", 
+            RowBox[{
+             RowBox[{"Table", "[", 
+              RowBox[{
+               RowBox[{"If", "[", 
+                RowBox[{
+                 RowBox[{"i", "==", "k"}], ",", "u", ",", 
+                 RowBox[{"p", "[", 
+                  RowBox[{"[", "k", "]"}], "]"}]}], "]"}], ",", 
+               RowBox[{"{", 
+                RowBox[{"k", ",", 
+                 RowBox[{"Length", "[", "p", "]"}]}], "}"}]}], "]"}], ",", 
+             "lam2", ",", "x", ",", "y"}], "]"}], ",", 
+           RowBox[{"{", 
+            RowBox[{"u", ",", "2"}], "}"}], ",", 
+           RowBox[{"p", "[", 
+            RowBox[{"[", "i", "]"}], "]"}], ",", 
+           RowBox[{"Terms", "->", "10"}]}], "]"}]}], ";", 
+        "\[IndentingNewLine]", "\[IndentingNewLine]", 
+        RowBox[{"If", "[", 
+         RowBox[{
+          RowBox[{"i", "!=", "j"}], ",", "\[IndentingNewLine]", 
+          RowBox[{"ND", "[", 
+           RowBox[{
+            RowBox[{"F", "[", "v", "]"}], ",", "v", ",", 
+            RowBox[{"p", "[", 
+             RowBox[{"[", "j", "]"}], "]"}]}], "]"}], ",", 
+          "\[IndentingNewLine]", "G"}], "]"}]}]}], "\[IndentingNewLine]", 
+      "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", 
+   RowBox[{"H0", "=", 
+    RowBox[{"Table", "[", 
+     RowBox[{
+      RowBox[{"HCost", "[", 
+       RowBox[{"i", ",", "j", ",", "p"}], "]"}], ",", 
+      RowBox[{"{", 
+       RowBox[{"i", ",", 
+        RowBox[{"Length", "[", "p", "]"}]}], "}"}], ",", 
+      RowBox[{"{", 
+       RowBox[{"j", ",", 
+        RowBox[{"Length", "[", "p", "]"}]}], "}"}]}], "]"}]}]}]}]], "Input",
+ CellChangeTimes->{{3.980220845242152*^9, 3.9802210388102617`*^9}, {
+  3.980221083420125*^9, 3.980221210743389*^9}, {3.9802212959394293`*^9, 
+  3.980221347003302*^9}, {3.9802214368899813`*^9, 3.980221498114513*^9}, {
+  3.980221552212894*^9, 3.980221718470955*^9}, {3.980221751831498*^9, 
+  3.980221754488792*^9}, {3.980222174396463*^9, 3.9802222588377943`*^9}, {
+  3.9802223438606653`*^9, 
+  3.9802223486758127`*^9}},ExpressionUUID->"db853e79-d0e4-4a86-8783-\
+98f6bc5933c4"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+   RowBox[{"{", 
+    RowBox[{"0.7673411749127825`", ",", 
+     RowBox[{"-", "0.00004960801306647878`"}], ",", "0.889749135171745`", ",",
+      "3.3212700921235605`"}], "}"}], ",", 
+   RowBox[{"{", 
+    RowBox[{
+     RowBox[{"-", "0.0000496078733072483`"}], ",", "0.0004804061483556116`", ",", 
+     RowBox[{"-", "0.0000624478957252519`"}], ",", 
+     RowBox[{"-", "0.00022472631280071956`"}]}], "}"}], ",", 
+   RowBox[{"{", 
+    RowBox[{"0.8897491344529453`", ",", 
+     RowBox[{"-", "0.00006244763646673682`"}], ",", "1.2193302465332434`", ",",
+      "4.2116021632604665`"}], "}"}], ",", 
+   RowBox[{"{", 
+    RowBox[{"3.3212689431121136`", ",", 
+     RowBox[{"-", "0.0002247255362663118`"}], ",", "4.211598490108457`", ",", 
+     "15.106151994325286`"}], "}"}]}], "}"}]], "Output",
+ CellChangeTimes->{3.980221740709745*^9, 3.9802217766213427`*^9, 
+  3.9802222890883713`*^9},
+ CellLabel->"Out[56]=",ExpressionUUID->"9ac024a3-d751-44de-bc53-6ad3de884f93"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"MatrixForm", "[", "H0", "]"}]], "Input",
+ CellChangeTimes->{{3.980221756962586*^9, 3.980221761796362*^9}, {
+  3.9802218237372227`*^9, 3.98022182432813*^9}},
+ CellLabel->"In[57]:=",ExpressionUUID->"c25363e3-6dd7-423a-81cd-c2c937573e5b"],
+
+Cell[BoxData[
+ TagBox[
+  RowBox[{"(", "\[NoBreak]", GridBox[{
+     {"0.7673411749127825`", 
+      RowBox[{"-", "0.00004960801306647878`"}], "0.889749135171745`", 
+      "3.3212700921235605`"},
+     {
+      RowBox[{"-", "0.0000496078733072483`"}], "0.0004804061483556116`", 
+      RowBox[{"-", "0.0000624478957252519`"}], 
+      RowBox[{"-", "0.00022472631280071956`"}]},
+     {"0.8897491344529453`", 
+      RowBox[{"-", "0.00006244763646673682`"}], "1.2193302465332434`", 
+      "4.2116021632604665`"},
+     {"3.3212689431121136`", 
+      RowBox[{"-", "0.0002247255362663118`"}], "4.211598490108457`", 
+      "15.106151994325286`"}
+    },
+    GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
+    GridBoxSpacings->{"Columns" -> {
+        Offset[0.27999999999999997`], {
+         Offset[0.7]}, 
+        Offset[0.27999999999999997`]}, "Rows" -> {
+        Offset[0.2], {
+         Offset[0.4]}, 
+        Offset[0.2]}}], "\[NoBreak]", ")"}],
+  Function[BoxForm`e$, 
+   MatrixForm[BoxForm`e$]]]], "Output",
+ CellChangeTimes->{3.980221824708803*^9, 3.9802222891345053`*^9},
+ CellLabel->
+  "Out[57]//MatrixForm=",ExpressionUUID->"60a8e5e9-338f-4169-8f24-\
+577a50a6f32e"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"MatrixForm", "[", 
+  RowBox[{"Inverse", "[", "H0", "]"}], "]"}]], "Input",
+ CellChangeTimes->{{3.980222302130584*^9, 3.980222304864068*^9}},
+ CellLabel->"In[58]:=",ExpressionUUID->"3c1e775c-af2e-4c18-8894-85a24d299c82"],
+
+Cell[BoxData[
+ TagBox[
+  RowBox[{"(", "\[NoBreak]", GridBox[{
+     {"124.2107378939491`", "0.008819362051593838`", "99.68227366780131`", 
+      RowBox[{"-", "55.10069583812897`"}]},
+     {"0.008763138449467222`", "2081.5865173231978`", 
+      RowBox[{"-", "0.002461036440763724`"}], "0.029726128780833826`"},
+     {"99.68361653826865`", 
+      RowBox[{"-", "0.002420312442278083`"}], "102.15456747369905`", 
+      RowBox[{"-", "50.39738859190163`"}]},
+     {
+      RowBox[{"-", "55.10103654542312`"}], "0.029702306550748403`", 
+      RowBox[{"-", "50.39706092462403`"}], "26.23155138113158`"}
+    },
+    GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
+    GridBoxSpacings->{"Columns" -> {
+        Offset[0.27999999999999997`], {
+         Offset[0.7]}, 
+        Offset[0.27999999999999997`]}, "Rows" -> {
+        Offset[0.2], {
+         Offset[0.4]}, 
+        Offset[0.2]}}], "\[NoBreak]", ")"}],
+  Function[BoxForm`e$, 
+   MatrixForm[BoxForm`e$]]]], "Output",
+ CellChangeTimes->{3.98022230519285*^9},
+ CellLabel->
+  "Out[58]//MatrixForm=",ExpressionUUID->"635785cc-da30-4ada-b0a8-\
+311ff9635453"]
+}, Open  ]]
+},
+WindowSize->{1366.5, 683.25},
+WindowMargins->{{7.5, Automatic}, {Automatic, 24}},
+FrontEndVersion->"14.3 for Linux x86 (64-bit) (July 8, 2025)",
+StyleDefinitions->"Default.nb",
+ExpressionUUID->"9ea52935-3f71-45ad-8b9c-8ebfe3a2f1d5"
+]
+(* End of Notebook Content *)
+
+(* Internal cache information *)
+(*CellTagsOutline
+CellTagsIndex->{}
+*)
+(*CellTagsIndex
+CellTagsIndex->{}
+*)
+(*NotebookFileOutline
+Notebook[{
+Cell[554, 20, 3290, 81, 413, "Input",ExpressionUUID->"db312a24-878f-40c8-962f-6aedbdf27fe4"],
+Cell[CellGroupData[{
+Cell[3869, 105, 3458, 67, 258, "Input",ExpressionUUID->"f22ff9ef-193c-434b-91da-e9e8b3175259"],
+Cell[7330, 174, 2011, 34, 97, "Output",ExpressionUUID->"a7f74a5e-642c-4b21-9dc6-dc4532fcd2dd"],
+Cell[9344, 210, 799, 12, 34, "Output",ExpressionUUID->"6fe63fdd-9529-4ff3-95a7-2b315610cc0a"],
+Cell[10146, 224, 354, 7, 34, "Output",ExpressionUUID->"9174505c-bb9d-4cfc-a204-c6c464fab952"],
+Cell[10503, 233, 196, 3, 60, "Output",ExpressionUUID->"dfb9c680-45f0-495e-bc68-7c102a014c4a"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[10736, 241, 250, 4, 30, "Input",ExpressionUUID->"73502887-6057-4d13-881a-8932ae266804"],
+Cell[10989, 247, 193, 2, 33, "Output",ExpressionUUID->"603f6210-ea83-425d-9d61-d60dffe5c701"]
+}, Open  ]],
+Cell[11197, 252, 230, 3, 30, "Input",ExpressionUUID->"417b4445-9108-40f3-a80b-fda472fd9d54"],
+Cell[CellGroupData[{
+Cell[11452, 259, 150, 2, 29, "Input",ExpressionUUID->"3ca07280-6eb4-41d4-8ab2-dfe9bbf7b8be"],
+Cell[11605, 263, 331, 7, 34, "Output",ExpressionUUID->"029d324b-6436-44a3-8803-14c8381b5896"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[11973, 275, 1869, 46, 73, "Input",ExpressionUUID->"b0ded771-dccd-45f6-a46c-9b2ff466801f"],
+Cell[13845, 323, 723, 12, 45, "Output",ExpressionUUID->"2d1c80cf-dfd3-4a34-8f1d-5ca8d4f13801"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[14605, 340, 3534, 88, 247, "Input",ExpressionUUID->"db853e79-d0e4-4a86-8783-98f6bc5933c4"],
+Cell[18142, 430, 995, 22, 57, "Output",ExpressionUUID->"9ac024a3-d751-44de-bc53-6ad3de884f93"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[19174, 457, 259, 4, 30, "Input",ExpressionUUID->"c25363e3-6dd7-423a-81cd-c2c937573e5b"],
+Cell[19436, 463, 1181, 30, 111, "Output",ExpressionUUID->"60a8e5e9-338f-4169-8f24-577a50a6f32e"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[20654, 498, 243, 4, 30, "Input",ExpressionUUID->"3c1e775c-af2e-4c18-8894-85a24d299c82"],
+Cell[20900, 504, 1117, 27, 111, "Output",ExpressionUUID->"635785cc-da30-4ada-b0a8-311ff9635453"]
+}, Open  ]]
+}
+]
+*)
+

python/logistic/data_utils.py

@@ -5,6 +5,7 @@
 """
 
 import numpy as np
+import pandas as pd
 
 """
   Extract data fom dictionaries for specific organ
@@ -62,3 +63,50 @@ def within_bounds(vector: np.ndarray, bounds: np.ndarray) -> bool:
     
     return np.apply_along_axis(lambda x: np.all((x >= bounds[:, 0]) & (x <= bounds[:, 1])), -1, vector)
 
+
+
+def  prepare_data(x0, y0, drop_mask):
+    """
+        Prepare data for plotting and fitting by creating a DataFrame with different subsets and scales.
+
+        Parameters:
+        - x0: Original feature vector (e.g., SUV percentiles).
+        - y0: Original target vector (e.g., patient states).
+        - drop_mask: Boolean array indicating which samples to drop for the "TRIM" dataset.
+
+        Returns:
+        dataframe, scales, datasets:
+        - df_data: A pandas DataFrame containing the prepared data with columns ['X', 'Y', 'scale', 'dataset'].
+        - scales: List of scale names.
+        - datasets: List of dataset names.
+    """
+    
+    datasets = ["FULL", "TRIM"]
+    scales = ["plain", "log"]
+
+    lst = []
+    for dataset in datasets:
+        for scale in scales:
+
+            if dataset == "TRIM":
+                y_tmp = y0[~drop_mask]
+                x_tmp = x0[~drop_mask]
+            else:
+                y_tmp = y0[:]
+                x_tmp = x0[:]
+            
+            if scale == "log":
+                # remove non-positive values for log scale
+                valid_mask = (x_tmp > 0)
+                y_tmp = y_tmp[valid_mask]
+                x_tmp = x_tmp[valid_mask]
+                x_tmp = np.log(x_tmp)
+
+            df_tmp = pd.DataFrame({
+                "X" : x_tmp,
+                "Y" : y_tmp ,
+                "scale" : scale,
+                "dataset": dataset})
+            lst.append(df_tmp)
+
+    return pd.concat(lst, ignore_index=True)

python/logistic/logit_utils_gen.py

@@ -51,7 +51,7 @@
     with respect to parameters pars. We can have regularization term in cost function
     and this case we minimize cost function:
 
-        cost(pars) = nllf(pars) + lambda ||pars||^2
+        cost(pars) = nllf(pars) + lambda_0*|pars| + lambda_1*|pars|_2^2
 
     Author: Martin Horvat, January 2026
 """
@@ -150,7 +150,7 @@ class LogisticPolyRegression:
 
         self.degree = degree
         self.mono = mono
-        self.big = 200
+        self.big = 1e3
         self.small = 1e-8
         self.lam = lam
 
@@ -265,7 +265,7 @@ class LogisticPolyRegression:
             nllf           : if jac is false
             (nllf, grad)   : if jac is true
     """
-    def nllf(self, x, y, pars, jac = False):
+    def get_nllf(self, x, y, pars, jac = False):
 
         beta = self.get_beta(pars)
         X = np.column_stack([x**i for i in range(len(beta))])
@@ -310,9 +310,9 @@ class LogisticPolyRegression:
     """
         Cost function
     """
-    def cost(self, x, y, pars, jac = False):
+    def get_cost(self, x, y, pars, jac = False):
 
-        val = self.nllf(x, y, pars, jac)
+        val = self.get_nllf(x, y, pars, jac)
         
         if self.lam is not None:
             pen = self.penalty(pars, jac)
@@ -370,23 +370,23 @@ class LogisticPolyRegression:
         if method == "local":
 
             pars0 = self.get_est_pars(x, y)
-            cf = lambda pars: self.cost(x, y, pars, jac = True)
+            cf = lambda pars: self.get_cost(x, y, pars, jac = True)
             res = scipy.optimize.minimize(cf, x0 = pars0, method = 'L-BFGS-B', 
                                           jac = True, bounds = bnds, tol=1e-12)
 
         elif method == "diff_evol":
             
-            cf = lambda pars: self.cost(x, y, pars, jac = False)
+            cf = lambda pars: self.get_cost(x, y, pars, jac = False)
             res = scipy.optimize.differential_evolution(cf, bounds = bnds, 
                                                         tol = 1e-8, polish = False)
 
-            cf = lambda pars: self.cost(x, y, pars, jac = True)
+            cf = lambda pars: self.get_cost(x, y, pars, jac = True)
             res = scipy.optimize.minimize(cf, x0 = res.x, method = 'L-BFGS-B', 
                                           jac = True, bounds = bnds, tol=1e-12)
 
         elif method == "anneal":
 
-            cf = lambda pars: self.cost(x, y, pars, jac = False)
+            cf = lambda pars: self.get_cost(x, y, pars, jac = False)
             res = scipy.optimize.dual_annealing(cf, bounds = bnds)
         
         else:
@@ -428,7 +428,7 @@ class LogisticPolyRegression:
         p = self.model(x, pars)
 
         # log likelihood
-        llf = -self.nllf(x, y, pars)
+        llf = -self.get_nllf(x, y, pars)
         
         # information criteria
         k, n = len(pars), len(x)
@@ -454,33 +454,18 @@ class LogisticPolyRegression:
                 "chi2": chi2, 
                 "p-value(chi2)": p_val,  # not very useful
                 "n": n, "k": k, "dof": dof}
-
+    
     """
-        Calculation of asymptotic variance-covariance matrix of regression 
-        parameters pars
-
-           cov_{asymp}[pars] = H^{-1}
-
-        where H is hessian of nllf 
-
-           H = [d^2(nllf)/(d(pars)_a d(pars)_b ]_{a,b}
-        
-        for the logistic regression of the polynomial model:
-        
-            log(f(x)/(1 - f(x))) ~ sum_{i=0}^degree b_i(pars) x^i
+        Calculate hessian of cost function with respect to parameters pars
 
         Input:
             x: array of n floats
+            y: array of n int in {0,1}
             pars: array of r = degree+1 floats, model parameters
-            
         Return:
-            array of rxr floats; r = degree + 1
-        
-        Ref:
-          * https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients
-          * https://goodboychan.github.io/machine_learning/2020/09/14/02-Regularized-likelihood-methods.html
+            H matrix of shape (r, r)
     """
-    def cov(self, x, y, pars):
+    def get_cost_hessian(self, x, y, pars):
         
         # coefficients
         beta = self.get_beta(pars)
@@ -488,31 +473,59 @@ class LogisticPolyRegression:
         # design matrix -- add column of 1's at the beginning of your X_train matrix
         X = np.column_stack([x**i for i in range(len(beta))])
 
-        # Jacobian J = [dbeta_i/dpars_j]_{ij}
+        # Jacobian J = [dbeta_i/dpars_a]_{ia}
+        # Hessian H = [d^2 beta_i/(d(pars_a) d(pars_b))]_{i,a,b}
         J, H = self.get_jac_beta(pars, hess = True)
 
         # signs
         s = 2.0*y - 1
         
         # decision function for conditional probability Prob(Y = y| x)
-        F = s*(X @ beta)
+        V = s[:,None]*X
+        F = V @ beta
 
         # probabilities p_i = P(Y=y_i | x_i)
         p = safe_expit(F)
         q = 1 - p
 
         # calculate hessian
-        L = X @ J
-        H = (L.T*(q*p))@L - np.tensordot((s*q)@X, H, axes = ([0], [0]))
+        L = V @ J
+        H = (L.T*(q*p))@L - np.tensordot(q@V, H, axes = ([0], [0]))
 
         if self.lam is not None: 
-            Hp = H + 2*self.lam[1]*np.eye(len(pars))   # H' = H + lambda id 
-            iHp = np.linalg.inv(Hp)                    # inv(H')
+            return H + 2*self.lam[1]*np.eye(len(pars))   # H' = H + lambda id 
+       
+        return H
 
-            return iHp@H@iHp
+    """
+        Calculation of asymptotic variance-covariance matrix of regression 
+        parameters pars
+
+           cov_{asymp}[pars] = H^{-1}
+
+        where H is hessian of nllf 
+
+           H = [d^2(nllf)/(d(pars)_a d(pars)_b ]_{a,b}
+        
+        for the logistic regression of the polynomial model:
+        
+            log(f(x)/(1 - f(x))) ~ sum_{i=0}^degree b_i(pars) x^i
+
+        Input:
+            x: array of n floats
+            pars: array of r = degree+1 floats, model parameters
+            
+        Return:
+            array of rxr floats; r = degree + 1
+        
+        Ref:
+          * https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients
+          * https://goodboychan.github.io/machine_learning/2020/09/14/02-Regularized-likelihood-methods.html
+    """
+    def get_cov(self, x, y, pars):
         
         # covariance matrix C_params = H^-1
-        return np.linalg.inv(H)
+        return np.linalg.pinv(self.get_cost_hessian(x, y, pars))
 
     """
         Calculating standard errors fo model parameters for normal distribution of parameters:
@@ -529,7 +542,7 @@ class LogisticPolyRegression:
     def get_SE_pars_normal(self, cov_pars):
         
         # computing standard errors of parameters
-        return np.sqrt(np.diag(cov_pars))
+        return np.sqrt(np.clip(np.diag(cov_pars),a_min=0, a_max = None))
     
     """
         Calculating quantiles of the model parameters at given probabilities p 
@@ -549,7 +562,7 @@ class LogisticPolyRegression:
         
         # mean and standard variance parameters
         locs = mean_pars
-        scales = np.sqrt(np.diag(cov_pars))
+        scales = np.sqrt(np.clip(np.diag(cov_pars), a_min=0, a_max=None))
 
         # computing quantiles of parameters
         return locs + np.outer(scipy.stats.norm.ppf(probs), scales)
@@ -689,7 +702,7 @@ class LogisticPolyRegression:
         pars = res["pars"]
 
         # covariance matrix of parameters
-        cov = self.cov(x, y, pars)
+        cov = self.get_cov(x, y, pars)
 
         # init random generator
         rng = np.random.default_rng(seed)
@@ -844,4 +857,5 @@ class LogisticPolyRegression:
             lst.append(res_fit["pars"])
             if len(lst) == m: break
 
-        return np.array(lst)
+        return np.array(lst)
+    

python/logistic/mono_cubic1.py

@@ -3,7 +3,8 @@
         first parametrization
     We define the polynomial:
           poly(x) = sum_i beta_i * x^i
-    where the coefficients beta_i are parameterized by the vector 'pars'.
+    where the coefficients beta_i are parameterized by the vector 'pars':
+        beta_i is a function pars
 """
 
 import numpy as np

python/logistic/mono_cubic2.py

@@ -1,9 +1,13 @@
 """ 
     Function supporting monotonic cubic polynomials: 
         second parametrization
+    which looks more stable than the first on in mono_cubic1.py
+
     We define the polynomial:
           poly(x) = sum_i beta_i * x^i
-    where the coefficients beta_i are parameterized by the vector 'pars'.
+    where the coefficients beta_i are parameterized by the vector 'pars',
+    meaning:
+            beta_i is a function pars
 """
 import numpy as np