estimation of values using NMaximize over symbolic array entries

I'm trying to generate instances of a dataset that satisfy specific conditions, based on the following equation:

Each of four non-zero terms (term 1, A; term 3, B; term 4, C; term 5, D) is a scalar. Based on previous help in the forum, I've managed to put together the following script to generate instances for which each of the four terms predominate:

Y=Array[y,{3,3}];
groupmeans=Mean/@Y;
mean=Simplify[Mean[groupmeans]];

a=Simplify[Total[(groupmeans-mean)^4]];
b=Simplify[6  (groupmeans-mean)^2 .Total[(Y-groupmeans)^2,{2}]];
c=Simplify[4 (groupmeans-mean) .Total[(Y-groupmeans)^3,{2}]];
d=Simplify[Total[Total[(Y-Mean[Mean[Y]])^4]]];
vars=Flatten[Y];
skb=NMaximize[{Boole[a>b]+Boole[a>c]+Boole[a>d]},vars,MaxIterations->10^5]
skv=NMaximize[{Boole[b>a]+Boole[b>c]+Boole[b>d]},vars,MaxIterations->10^6]
sks=NMaximize[Boole[c>a]+Boole[c>b]+Boole[c>d],vars,MaxIterations->10^5]
skw=NMaximize[Boole[d>a]+Boole[d>b]+Boole[d>c],vars,MaxIterations->10^5]


As have understood the input, the above generates a 3 x 3 ($$j$$ x $$i$$) array of symbolic input values, and symbolically assigns a calculation for the group means and grand mean. These symbolic placeholders are then passed to the respective terms of the above equation, a-d, and the conditions in NMaximize should, in principle, produce 9 values that satisfy the corresponding conditions.

When I run the above script, I can generate sets of values for skv and sks. In the above example, the numerical output for skv has terms that are enormous (10^30000). The sks values are more reasonable. However, when I plug them back into the equation, the conditions do not hold.

I am assuming the script I have is either configured incorrectly in the NMaximize calls, or in an assignment of a-d terms that is not correct, given the equation. Is there a way to fix the script to generate the desired output? Ideally, the values should only be positive numbers (0, $$\infty$$).

I have Mathematica 9.

• Comments are not for extended discussion; this conversation has been moved to chat. – Kuba Nov 3 '18 at 19:04
• The question was cardinally edited. In particular, FindInstance was replaced by NMaximize and the title was changed twice. This is not a good practice, Please, don't delete my comment dropbox.com/s/fgg98uzkem5af7s/screen03.11.18.docx?dl=0 – user64494 Nov 5 '18 at 19:51
• @user64494 that is correct but it does not matter anymore, does it? It can only confuse future readers. You are free to join the linked chat to convince us that this comment should stay, you didn't try so far, despite my question. – Kuba Nov 5 '18 at 21:18

This tries to translate your image at the very top of your original post directly into Mathematica code.

k = 3; ni = 3;
yBariDot[i_] := Sum[y[i, j], {j, 1, ni}]/ni;
yBarDotDot[] := Sum[y[i, j], {i, 1, k}, {j, 1, ni}]/(k*ni);

SKtot = Sum[(yBariDot[i] - yBarDotDot[])^4, {i, 1, k}, {j, 1, ni}] +
4*Sum[(yBariDot[i] - yBarDotDot[])^3*(y[i, j] - yBariDot[i]), {i, 1, k}, {j, 1, ni}] +
6*Sum[(yBariDot[i] - yBarDotDot[])^2*(y[i, j] - yBariDot[i])^2, {i, 1, k}, {j, 1, ni}] +
4*Sum[(yBariDot[i] - yBarDotDot[])*(y[i, j] - yBariDot[i])^3, {i, 1, k}, {j, 1, ni}] +
Sum[(y[i, j] - yBariDot[i])^4, {i, 1, k}, {j, 1, ni}];

(* and now the four terms from SKtot which you want to have satisfy inequalities *)
a = Sum[(yBariDot[i] - yBarDotDot[])^4, {i, 1, k}, {j, 1, ni}];
b = 6*Sum[(yBariDot[i] - yBarDotDot[])^2*(y[i, j] - yBariDot[i])^2, {i, 1, k}, {j, 1, ni}];
c = 4*Sum[(yBariDot[i] - yBarDotDot[])*(y[i, j] - yBariDot[i])^3, {i, 1, k}, {j, 1, ni}];
d = Sum[(y[i, j] - yBariDot[i])^4, {i, 1, k}, {j, 1, ni}];

vars = Flatten[Array[y, {k, ni}]];
conditions = Map[# > 0 &, vars];
skb = NMaximize[{Boole[a > b] + Boole[a > c] + Boole[a > d], conditions}, vars]
skv = NMaximize[{Boole[b > a] + Boole[b > c] + Boole[b > d], conditions}, vars]
sks = NMaximize[{Boole[c > a] + Boole[c > b] + Boole[c > d], conditions}, vars]
skw = NMaximize[{Boole[d > a] + Boole[d > b] + Boole[d > c], conditions}, vars]


Notice this now incorporates your added condition that all y[i,j]>0.

When I run that I get this output

{1., {y[1, 1] -> 1.42388, y[1, 2] -> 1.75502, y[1, 3] -> 2.04228,
y[2, 1] -> 1.09188, y[2, 2] -> 0.620853, y[2, 3] -> 1.6464,
y[3, 1] -> 1.48234, y[3, 2] -> 1.87998, y[3, 3] -> 1.48247}}

{3., {y[1, 1] -> 2.04438, y[1, 2] -> 2.27041, y[1, 3] -> 2.42588,
y[2, 1] -> 1.10609, y[2, 2] -> 0.269069, y[2, 3] -> 2.3213,
y[3, 1] -> 1.47011, y[3, 2] -> 2.54533, y[3, 3] -> 1.71664}}

{2., {y[1, 1] -> 1.67519, y[1, 2] -> 0.478676, y[1, 3] -> 1.60231,
y[2, 1] -> 1.32746, y[2, 2] -> 1.37434, y[2, 3] -> 1.7539,
y[3, 1] -> 1.86218, y[3, 2] -> 1.00843, y[3, 3] -> 1.68511}}

{3., {y[1, 1] -> 5.18063, y[1, 2] -> 1.90695, y[1, 3] -> 4.17558,
y[2, 1] -> 3.71691, y[2, 2] -> 0.596049, y[2, 3] -> 4.6259,
y[3, 1] -> 3.90677, y[3, 2] -> 2.79298, y[3, 3] -> 5.49723}}


Now I use that output to test each result for validity.

In[16]:= {Boole[a > b] + Boole[a > c] + Boole[a > d],
conditions} /. {y[1, 1] -> 1.4238765790809782,
y[1, 2] -> 1.75502338307695, y[1, 3] -> 2.042280852749116,
y[2, 1] -> 1.0918808384496819, y[2, 2] -> 0.6208527474852398,
y[2, 3] -> 1.6463955261954661, y[3, 1] -> 1.482339964456537,
y[3, 2] -> 1.8799832736930584, y[3, 3] -> 1.4824728315771387}

Out[16]= {1, {True, True, True, True, True, True, True, True, True}}

In[17]:= {Boole[b > a] + Boole[b > c] + Boole[b > d],
conditions} /. {y[1, 1] -> 2.0443840549247256,
y[1, 2] -> 2.2704092153250945, y[1, 3] -> 2.4258771603051663,
y[2, 1] -> 1.1060948904114927, y[2, 2] -> 0.26906917232850847,
y[2, 3] -> 2.3213006883182383, y[3, 1] -> 1.4701091327440556,
y[3, 2] -> 2.545326429416577, y[3, 3] -> 1.7166395917271662}

Out[17]= {3, {True, True, True, True, True, True, True, True, True}}

In[18]:= {Boole[c > a] + Boole[c > b] + Boole[c > d],
conditions} /. {y[1, 1] -> 1.6751902804210805,
y[1, 2] -> 0.47867608781606996, y[1, 3] -> 1.6023121928338604,
y[2, 1] -> 1.3274636095241925, y[2, 2] -> 1.3743360783660326,
y[2, 3] -> 1.7539049108235778, y[3, 1] -> 1.8621812898930021,
y[3, 2] -> 1.008433224445793, y[3, 3] -> 1.685113152710088}

Out[18]= {2, {True, True, True, True, True, True, True, True, True}}

In[19]:= {Boole[d > a] + Boole[d > b] + Boole[d > c],
conditions} /. {y[1, 1] -> 5.180634908092935,
y[1, 2] -> 1.9069494885810279, y[1, 3] -> 4.1755803675439225,
y[2, 1] -> 3.7169070293128708, y[2, 2] -> 0.5960494074966672,
y[2, 3] -> 4.6258987168248735, y[3, 1] -> 3.9067720907993273,
y[3, 2] -> 2.7929797597944135, y[3, 3] -> 5.497225073531151}

Out[19]= {3, {True, True, True, True, True, True, True, True, True}}


That does not use hundreds of thousands or millions of iterations where NMaximize might crawl off towards infinity, perhaps because there is no real solution and it is just groping for floating point roundoff errors. But it does seem to correctly find values for your variables which satisfy some of your conditions.

If you have a little time and could look at my translation from your graphic to my Mathematica code and see if there are possibly any errors in that and then explain those to me then I would be very appreciative. It might be interesting to compare your values for a,b,c,d and these new values for a,b,c,d when a solution is found using either one of the methods. If those are the same then that might provide a little more confidence in both, if they differ then that might provide a reason to go look for misunderstandings.

Please test this without mercy and see if you can expose any of my mistakes. Thank you

• Comments are not for extended discussion; this conversation has been moved to chat. – Kuba Nov 5 '18 at 21:01