LinearModelFit with the dependent variable given in intervals instead of exact values

Question

I have data of the form $\{x_{i1},...,x_{ik},\{a_i,b_i\}\}$, $i=1,...,N$, and I want to find a function $f(x_1,...,x_k)$ (say, linear) such that violation of $a_i\le f(x_{i1},...,x_{ik})\le b_i$ for all $i$ is as small as possible. Does LinearModelFit or FindFit or some similar command have this functionality? If not, what would be an optimal code to deal with this?

As promised, here is a small sample of the kind of data I've got.

{
{36.8,-0.34,98.7,23.8,-0.71,92.6,33.7,-0.42,96.3,0,15.42,-0.05,-0.67,-0.42,-0.16,{2,24}},
{33.7,-0.42,98.3,11,-1.22,87.4,16.2,-0.99,91.3,0,14.64,-0.68,-0.75,-0.9,-0.76,{2,24}},
{6.1,-1.55,93.6,10.7,-1.24,85.9,36.5,-0.35,97,0,15.27,0.02,-1.57,-1.02,-0.07,{0,1}},
{19.7,-0.85,96.6,25.8,-0.65,92.9,44.8,-0.13,98.8,0,15.06,-0.14,-1.02,-0.7,-0.1,{2,24}},
{0.2,-2.93,89.5,0.1,-3.19,67.6,8.1,-1.4,86.7,0,14.84,-1.71,-3.21,-3.13,-1.87,{0,0}},
{37.7,-0.31,98.8,5,-1.65,81.8,4.6,-1.68,84.5,0,14.5,-1.47,-0.23,-1.17,-1.46,{0,0}},
{13.3,-1.11,95.3,7.3,-1.46,78.4,17.2,-0.94,89,0,12.2,-2.15,-1.43,-2.28,-1.57,{0,0}},
{7.5,-1.44,94.5,2,-2.05,77.5,7.3,-1.45,86.6,0,14.85,-0.98,-1.46,-1.56,-1.21,{0,0}},
{30.3,-0.52,97.7,22.1,-0.77,87.8,35.4,-0.38,95.3,0,13.17,-1.24,-0.71,-1.2,-0.93,{0,0}},
{8.2,-1.39,94.4,9.9,-1.29,86.1,29.8,-0.53,95.3,0,14.64,-0.43,-1.53,-1.24,-0.43,{0,0}},
{14.1,-1.08,95.9,17.7,-0.93,90,32.3,-0.46,96.1,0,16.19,0.32,-1.06,-0.4,0.13,{0,0}}
{15.7,-1.01,95.8,9.4,-1.32,85,18.6,-0.89,92.2,0,14.41,-0.67,-1.03,-1.09,-0.72,{1,2}},
{30.9,-0.5,97.9,30.1,-0.52,94.1,42.5,-0.19,98.4,0,15.45,0.09,-0.6,-0.32,0.09,{1,2}},
{23.1,-0.74,97.2,14.3,-1.07,89.9,28,-0.58,94.8,0,15.93,-0.02,-1.26,-0.79,-0.19,{1,2}},
{6.3,-1.53,93.6,13.4,-1.11,87.4,43.9,-0.15,98.7,0,15.44,0.18,-1.53,-0.87,0.1,{0,1}},
{13.2,-1.12,95.6,15.9,-1,89.4,40.1,-0.25,97.8,0,15.31,0.01,-1.33,-0.82,-0.07,{0,1}},
{34.3,-0.41,98.6,19.2,-0.87,91.3,23.5,-0.72,94.2,0,16.27,-0.18,-0.34,-0.34,-0.25,{1,2}},
{24.7,-0.69,97.3,24.2,-0.7,92.4,40.4,-0.24,97.9,0,15.92,0.23,-1.05,-0.5,0.08,{0,1}},
{0.1,-3.18,87.4,0.7,-2.46,74,22.9,-0.74,93.3,0,15.3,-0.02,-3.33,-2.3,-0.35,{1,2}},
{16.3,-0.98,96.1,9.4,-1.32,86.2,23.5,-0.72,93.6,0,14.98,-0.34,-1.27,-1.02,-0.48,{0,1}},
{14.7,-1.05,96,2.7,-1.93,79.6,2.8,-1.9,85.3,0,14.17,-1.4,-1.21,-1.69,-1.48,{0,1}},
{25.1,-0.67,97.3,20.4,-0.83,91.1,33.8,-0.42,96.3,0,15.7,0.03,-1.05,-0.62,-0.12,{0,1}}
}

Answer 1

Here the cost is zero if inside an interval. If outside, the cost increases with the square of the distance to the nearest edge of the interval. I minimize the total cost with NMinimize:

vectors = Most /@ data;
intervals = Last /@ data;
vars = Array[x, Last@Dimensions@vectors];
penalty[point_, vec_, interval_] :=
 If[interval[[1]] <= point.vec <= interval[[2]], 0,
  Min[(point.vec - #)^2 & /@ interval]]
cost = Total[MapThread[penalty[vars, #1, #2] &, {vectors, intervals}]];
result = vars /. Last[NMinimize[cost, vars]]
membership = MapThread[IntervalMemberQ[Interval[#2], result.#1] &, {vectors, intervals}]
N[Count[membership, True]/Length[membership]]
(* result 28.5% *)

Unfortunately, under the dot-product / linear model only 28.5% of $f(\mathbf{x}_i)$ fell into the intervals. If I replace the complicated penalty with a simple square distance to midpoint...

penalty[point_, vec_, interval_] := (point.vec - (interval[[2]] - interval[[1]])/2)^2

... I still get 28.5%. Perhaps a nonlinear model is more appropriate.

Using the original penalty, if I create a DimensionReduction function of your vectors using the Method->"AutoEncoder" I can achieve 38.1% by simply adding these lines after the first line vectors = Most /@ data;

drf = DimensionReduction[vectors, Method -> "AutoEncoder"];
vectors = drf[vectors];

To use this on an unseen vector v you would apply result.drf[v]. Finally, and this is probably going into over-fitting territory as we haven't done any cross-validation, I've also constrained the $L_\infty$ norm of $\mathbf{f}$ and switched to the midpoint penalty. This takes us all the way to 42.8%. Putting it all together:

vectors = Most /@ data;
intervals = Last /@ data;
drf = DimensionReduction[vectors, Method -> "AutoEncoder"];
vectors = drf[vectors];
vars = Array[x, Last@Dimensions@vectors];
penalty[point_, vec_, interval_] := (point.vec - (interval[[2]] - interval[[1]])/2)^2
cost = Total[MapThread[penalty[vars, #1, #2] &, {vectors, intervals}]];
result = vars /. Last[NMinimize[{cost, Norm[vars, Infinity] < 2}, vars]]
membership = MapThread[IntervalMemberQ[Interval[#2], result.#1] &, {vectors, intervals}]
N[Count[membership, True]/Length[membership]]
(*result 42.8%*)