# Function to marginalize a table over given parameters

Let us say I have generated a table tab with entries {p1,..,pd,chi2[p1,..,pd]} of dimension {ngrid[1],..,ngrid[d],d+1} where chi2 is a function of d parameters. More precisely, tab is generated using the following code:

tab = Table[
{p1, p2, p3, p4, p5} = {pTab[1][[i1]], pTab[2][[i2]], pTab[3][[i3]],pTab[4][[i4]], pTab[5][[i5]]};
{p1, p2, p3, p4, p5, chi2[p1, p2, p3, p4, p5]}
, {i1, 1, ngrid[1]}, {i2,1,ngrid[2]}, {i3,1,ngrid[3]}, {i4,1,ngrid[4]}, {i5,1,ngrid[5]}]


where pTab[i] are the tables giving the pi parameters.

Now, I want to marginalize the table over some parameters. More precisely I want to marginalize the corresponding likelihood ($\exp(-\chi^2/2)$). This can be done over all parameters but p2 in the following way:

newtab = tab[[1, All, 1, 1, 1, {2, 6}]];
Do[
newtab[[j,2]]=-2Log[Total[Exp[-tab[[All,j,All,All,All,6]]/2],4]];
,{j, 1,ngrid[2]}];


If one wants to marginalize over all parameters but p1 and p3:

newtab = tab[[All, 1, All, 1, 1, {1, 3, 6}]]
Do[
newtab[[j, k, 3]] = -2Log[Total[Exp[-tab[[j,All,k,All,All,6]]/2],4]];
,{j, 1, ngrid[1]}, {k, 1, ngrid[3]}];


and so on.

It would be great to have a function that can do the tasks above and which works for any number d of parameters and which can accept any combination of parameters to marginalize over. Something like marginalize[{p2,p4,p5},tab] where tab is my table and {p2,p4,p5} are the parameters to marginalize over. Do you know how to code such a function?

-
I'm not sure whether you're asking about the math or about how to code it in Mathematica – belisarius has settled Jul 7 '13 at 15:02
Maybe you could provide us with a simple example: say two ro three dimensions and a clear explanation of what the "marginalization" is supposed to do. For example, it is not clear whether this is much more than just normalizing the column/row of the matrix. – bill s Jul 7 '13 at 15:20
marginalization is supposed to do exactly what shown in the examples. The problem is that I don't know how to code a function that can do flexibly those jobs. – Valerio Jul 7 '13 at 16:18
But Valero, we don't know what your function does because you haven't defined all the relevant data: what is tab? What is ngrid? If you provide an example where it is clear what you are after, you will be more likely to get help. – bill s Jul 7 '13 at 20:41
@bill s I added some explanations, is it now clear? – Valerio Jul 18 '13 at 12:31

I propose this function

marginalize[pp_, tab_] :=
Module[{newtab, id, d = Dimensions[tab][[-1]] - 1},
id = ConstantArray[1, d];
id[[pp]] = All;
newtab = Transpose[tab[[Sequence @@ id, {Sequence @@ pp, -1}]], Ordering[pp]];
With[{nn = Sequence @@ Table[Unique["n"], {Length[pp]}]},
With[{lim = Sequence @@ Transpose@{{nn}, Dimensions[tab][[pp]]},
mm = (id = ConstantArray[All, d]; id[[pp]] = {nn}; Sequence @@ id)},
Do[newtab[[nn, -1]] = -2 Log[Total[Exp[-tab[[mm, -1]]/2], Infinity]];, lim];
]];
newtab
];


It does not depend on any external parameters like the number of dimensions, ngrid, etc.

I think tab can be generated by

d = 5;
ClearAll[pTab];
Do[pTab[i] = Range[-1, 1, 0.1], {i, d}];
chi2 = Total[#^2] &; (* I don't know your chi2 and pTab *)

tab = Module[{T, n}, With[{F = chi2}, Replace[
Hold[Array[Append[T, F[T]] &, Length[pTab[#]] & /@ Range[d]]] /.
T -> Table[Hold[pTab[n][[Slot[n]]]] /. n -> i, {i, d}],
Hold[x_] :> x, {0, Infinity}, Heads -> True]]];

tab // DeveloperPackedArrayQ


True

It takes the number of dimensions d as parameter. It's quite complicated but it generate packed array, which is much faster.

-
thanks! only one thing: the parameters pp have to be sorted. Or, would it be possible to modify your function such that if i give pp={3,2}, then i get a table which has the third parameter in the first dimension and the second parameter in the second dimension? Otherwise i can do this later using: newtab = newtab[[All, All, {2, 1, 3}]]; newtab = Transpose[newtab, {2, 1}] – Valerio Sep 12 '13 at 10:01
@Valerio I add transposition to my code. I think it now supports unsorted pp`. Please check it. – ybeltukov Sep 12 '13 at 16:09
It's working alright, thanks! – Valerio Sep 13 '13 at 10:15