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

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}}
}

• Share (some of) your data. Would it be acceptable to minimize the difference between $f(x_i)$ and the middle of the $(a_i,b_i)$ interval for each $x_i$ set of values? Jun 19 '20 at 21:42
• @MarcoB I will try to pick some data. As for middles - no, if that would be preferable I could just use ordinary linear regression. The preference is to stay inside all intervals if possible, and this might mean that I need to get very close to $a_i$ or $b_i$. One could say that the penalty function is zero everywhere inside the intervals and grows rapidly as soon as I get out of them. Jun 19 '20 at 21:50
• That makes sense. Couldn't you then write an ancillary penalty function and minimize that? What form should the function $f$ have? Jun 20 '20 at 2:46
• @MarcoB Yes, I could. Do you propose to use NMinimize or NArgMin or something like that? In principle it is fine except I am not sure whether there is a global minimum, but I could try. Jun 20 '20 at 7:36

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%*)

• Great! Well, there are too many nonlinear models to choose from... But it is extremely useful to know anyway. Are you adding a constant term too? Intercept, I mean. Sorry, I should have it in the question - affine rather than linear. Jun 20 '20 at 14:18
• @მამუკაჯიბლაძე I have updated the answer. I tried a bunch of DimensionReduction methods too - AutoEncoder was the best. Jun 20 '20 at 14:30
• @მამუკაჯიბლაძე I'll have a look at adding a constant term. Jun 20 '20 at 14:50
• @მამუკაჯიბლაძე adding a constant term so we have const+point.vec did not improve it - in fact it got slightly worse and fell to 24% or lower in all combinations of using/not-using DimensionReduction and regularizing the norm. Jun 20 '20 at 15:00
• @მამუკაჯიბლაძე adding the degree of freedom might make it worse because of a local optimium, and on the other point - yes you may be right - I tried a completely harsh penalty of If[interval[[1]] <= point.vec <= interval[[2]], 0, 1] and I got 52% ! Some middle ground smooth function would be better to ensure that points that missed intervals were at least somewhat close. Jun 20 '20 at 15:33