# Do any of the Mathematica CopulaDistributions (or others) fit well this sampled bivariate copula with uniform marginals over [0,1]?

I'm currently developing a data set that consists of two $$50 \times 50$$ matrices, which I designate as q1 and Q1. I strongly believe (bordering on formal proof [cf. Corollary 1 in marginalinvariance]) that the ratio of Q1 to q1 constitutes a sample from a (symmetric) bivariate copula (alternatively termed "permutons or doubly-stochastic measures" WikiCopula) $$f(x,y)$$ with uniform marginals over [0,1] for which $$x$$ and $$y$$ ("Bloch radii" of quantum bit [qubit] systems--also "quadratic Casimir invariants" CasimirInvariants) are "repulsive" in nature RepulsiveBehavior, that is the 45-degree line $$x=y$$ has relatively low values.

Unfortunately, the sampling process to generate the data sets yields more lower values of $$x$$ and $$y$$ $$\in [0,1]$$ than higher values--so some trimming of q1 and Q1 might be appropriate in analyses.

So, my question is can the ratio of Q1 to q1 be well fitted by any of the Mathematica CopulaDistribution functions--or any others (outside of those in the Mathematica inventory)? A quick examination of the Mathematica Help discussion of CopulaDistribution, a multinormal copula seems most plausible of them. (What would its parameters be?)

The developing Q1 and q1 data sets as well as a plot--are displayed in BivariateCopulaRepulsive, which I hope is usable by interested parties. I intend to update q1 and Q1 as their entries further increase in size.

So, the $$\{i,j\}$$ cells of q1 and Q1 should be considered to correspond to $$x=\frac{2 i-1}{100}$$ and $$y=\frac{2 j-1}{100}$$.

Two calculations now inserted near the end of the linked Mathematica program show the near uniformity ($$\approx 0.122$$ over [0,1]) of the two marginal distributions. (It appears from certain published analyses, that the theoretical value is $$\frac{4}{33} \approx 0.1212$$.)

• perhaps, you might consider Frank copula.
– kglr
Dec 19 '21 at 23:17
• If I understand what you want, a copula distribution would just be a 2-dimensional function without any probabilistic interpretation? Why not fit a 2-dimensional surface to q1 and Q1 and then take the ratio of the smoothed curves (each of which would not have any zero values to mess up ratios)?
– JimB
Dec 20 '21 at 3:34
• Thanks for the nice thought, JimB--will pursue. For the moment, doing an update on ongoing data-generating analysis. Dec 20 '21 at 11:13
• JimB: I'm kind of rusty Mathematica-wise (been "incapacitated"--in a manner of speaking--for a year, long story [not Covid-related, though]). How might I fit these 2-D surfaces (do I need to choose some functional form)? Specific command suggestions would be appreciated. Dec 20 '21 at 11:56
• Thanks for the Frank copula suggestion--kgir! Off-hand from the Copula Kernels help page, the Multinormal (Binormal) with uniform marginals over [0,1] seems most promising (also conceptually--since Gaussian distributions are so universal in application). I will have to try to fit the parameters of the Binormal--suggestions as to how to proceed? Dec 20 '21 at 14:03

If only the ratios of the scaled probability density functions (with the scaling being the ratio of the two sample sizes), then I don't see why finding a reasonable joint density would matter (either using standard bivariate densities or something using a copula).

Here is a description as to how one would use a probability density function to fit a smooth curve to each marginal distribution of q1 and Q1 so that the ratio of those marginal pdf's can be better estimated. (Note that in standard statistical terms the ratio of two pdf's even if relatively constant across the domain is not a "uniform marginal". I don't know where that term comes from and I hope it's not been in common use in the physics community.)

So creating the "row" marginals from the binned bivariate samples as you have in your code the counts for each bin can displayed as follows:

nq1 = Total[q1] // Total;  (* Number of q1 samples *)
nQ1 = Total[Q1] // Total;  (* Number of Q1 samples *)

(* Row totals *)
p1 = Table[{(2 i - 1)/100, Total[q1[[i, All]]]}, {i, 1, 50}];
P1 = Table[{(2 i - 1)/100, Total[Q1[[i, All]]]}, {i, 1, 50}];

GraphicsRow[{ListPlot[p1, PlotLabel -> "q1 row marginal"],
ListPlot[P1, PlotLabel -> "Q1 row marginal"]}]

Both marginal distributions look like Gaussian distributions but I've found that a generalized gamma distribution fits a bit better. (And I'm sure there are either better or maybe a theoretically-based distribution.) But the point is to get a non-zero density so that ratios of the pdf's can be constructed.

So we fit a generalized gamma distribution and show the fits:

(* Fit a generalized gamma distribution *)
p1Fit = FindDistributionParameters[WeightedData[p1[[All, 1]], p1[[All, 2]]],
(* {α -> 0.793649, β -> 0.362512, γ -> 3.37196} *)

P1Fit = FindDistributionParameters[WeightedData[P1[[All, 1]], P1[[All, 2]]],
(* {α -> 0.801702, β -> 0.360379, γ -> 3.34451} *)

(* Plot marginals and generalized gamma distribution fit *)
GraphicsRow[{Show[ListPlot[p1, PlotLabel -> "q1 row marginal"],
Plot[(nq1*0.02) PDF[GammaDistribution[α, β, γ, 0] /. p1Fit, x], {x, 0, 1}]],
Show[ListPlot[P1, PlotLabel -> "Q1 row marginal"],
Plot[(nQ1*0.02) PDF[GammaDistribution[α, β, γ, 0] /. P1Fit, x], {x, 0, 1}]]}]

The fits look reasonable so now we take the ratio of the pdf's multiplied by the ratio of the sample sizes with two different PlotRange's.

(* Plot the ratio of the row counts and the marginal pdfs scaled by the total counts *)
ratio = p1;
ratio[[All, 2]] = P1[[All, 2]]/p1[[All, 2]];
Show[ListPlot[ratio, PlotRange -> {{0, 1}, {0, 0.2}}],
Plot[(nQ1*0.02) PDF[GammaDistribution[α, β, γ, 0] /. P1Fit, x]/
((nq1*0.02) PDF[GammaDistribution[α, β, γ, 0] /. p1Fit, x]),
{x, 0, 1}, PlotStyle -> Red]]

Show[ListPlot[ratio],
Plot[(nQ1*0.02) PDF[GammaDistribution[α, β, γ, 0] /. P1Fit, x]/
((nq1*0.02) PDF[GammaDistribution[α, β, γ, 0] /. p1Fit, x]),
{x, 0, 1}, PlotStyle -> Red]]

While I have not placed confidence bands for the ratio, I suspect that the ratio is only relatively constant as opposed to being strictly constant. The apparent lack-of-fit on the area to the right of 0.55 on the horizontal axes is likely due to the minimal number of samples in that area.

The average value of the ratio appears to just be the ratio of the sample sizes:

$$1532135/12378293 \approx 0.123776$$

The same approach can be taken with the column totals.

This is just an extended comment enhancing my lack of understanding of what you're asking and stating.

From downloading the definitions of q1 and Q2 from yesterday the marginal distributions all look very Gaussian rather than Uniform.

nq1 = Total[q1] // Total;
q1x = Table[{(2 i - 1)/100, Total[q1[[All, i]]]/nq1}, {i, Length[q1]}];
q1y = Table[{(2 i - 1)/100, Total[q1[[i, All]]]/nq1}, {i, Length[q1]}];
ListPlot[{q1c, q1r}, PlotStyle -> {{LightGray, PointSize[0.025]}, {Red, PointSize[0.01]}}]

nQ1 = Total[Q1] // Total;
Q1x = Table[{(2 i - 1)/100, Total[Q1[[All, i]]]/nQ1}, {i, Length[Q1]}];
Q1y = Table[{(2 i - 1)/100, Total[Q1[[i, All]]]/nQ1}, {i, Length[Q1]}];
ListPlot[{Q1c, Q1r}, PlotStyle -> {{LightGray, PointSize[0.025]}, {Red, PointSize[0.01]}}]

The marginal distributions are essentially identical for the 2 dimensions in q1 (and the same for Q1).

This suggests that a bivariate normal might be a good approximation for both q1 and Q1. Then one can take the ratio of the two bivariate probability density functions.

q1Data = Flatten[Table[{(2 i - 1)/100, (2 j - 1)/100, q1[[i, j]]/nq1}, {i, 50}, {j, 50}], 1];
Q1Data = Flatten[Table[{(2 i - 1)/100, (2 j - 1)/100, Q1[[i, j]]/nQ1}, {i, 50}, {j, 50}], 1];

q1xMean = q1Data[[All, 1]] . q1Data[[All, 3]];
q1yMean = q1Data[[All, 2]] . q1Data[[All, 3]];
q1xSD = Sqrt[(q1Data[[All, 1]]^2) . q1Data[[All, 3]] - q1xMean^2];
q1ySD = Sqrt[(q1Data[[All, 2]]^2) . q1Data[[All, 3]] - q1yMean^2];
q1Cov = (q1Data[[All, 1]] q1Data[[All, 2]]) . q1Data[[All, 3]] - q1xMean q1yMean;
q1Cor = q1Cov/(q1xSD q1ySD);

Q1xMean = Q1Data[[All, 1]] . Q1Data[[All, 3]];
Q1yMean = Q1Data[[All, 2]] . Q1Data[[All, 3]];
Q1xSD = Sqrt[(Q1Data[[All, 1]]^2) . Q1Data[[All, 3]] - Q1xMean^2];
Q1ySD = Sqrt[(Q1Data[[All, 2]]^2) . Q1Data[[All, 3]] - Q1yMean^2];
Q1Cov = (Q1Data[[All, 1]] Q1Data[[All, 2]]) . Q1Data[[All, 3]] - Q1xMean Q1yMean;
Q1Cor = Q1Cov/(Q1xSD Q1ySD);

qpdf = PDF[BinormalDistribution[{q1xMean, q1yMean}, {q1xSD, q1ySD}, q1Cor], {x, y}];
Qpdf = PDF[BinormalDistribution[{Q1xMean, Q1yMean}, {Q1xSD, Q1ySD}, Q1Cor], {x, y}];

ContourPlot[qpdf/Qpdf, {x, 0, 1}, {y, 0, 1}, PlotRange -> All,
Contours -> {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 2, 3, 4, 5, 6},
ContourLabels -> True, ContourShading -> False]

• I'm sorry JimB if we're not connecting or if I've been messing up the q1 and Q! in my updated postings of wolframcloud.com/obj/slater/Published/… But at the end of it all the non-degenerate marginal probabilities are shown as $\approx 0.12$ clearly exhibiting uniformity not Gaussianity. I get the marginals--as indicated -in the program by summing over either the rows and columns of Q1 and similarly for q1 and taking the ratio of these 50-length arrays. Dec 20 '21 at 23:19
• You might have noticed I can be very thick when it comes to understanding physics. I'll go to the updated postings and report back.
– JimB
Dec 20 '21 at 23:32
• I ran your code which produces p1 and P1 (the row sums). I would call those lists the row "marginals". What you are calling a "marginal" is the ratio of those two lists of counts. ListPlot[p1] and ListPlot[P1] produces the Gaussian looking distributions which I displayed above. I used the "density" (counts/total counts) for each list as opposed to just using the counts. If what you have are random samples from two distributions, then a ratio of the densities makes sense (to me). But you are using a ratio of the counts.
– JimB
Dec 20 '21 at 23:46
• The above approach (bivariate normal) could be used where instead of having the volume under the surfaces both be 1.0 (which is what the density does for you) the volume could be made to be the total count for q1 and Q1. I don't understand why the counts would be used but that's your call. The point is that you'd get two smooth surfaces (or two smooth curves if the row or column totals are used) that never go to zero and then the ratios could be taken.
– JimB
Dec 20 '21 at 23:49

In sec. 6.2.1 of the Wikipedia article, "Copula (probability theory)" WikiCopula, the Ali-Mikhail-Haq copula distribution is given by $$$$C_{\theta}=\frac{u v}{1-\theta (1-u)(1-v)}.$$$$

We start our analysis here with a $$50 \times 50$$ data matrix we denote S, that is the ratio of the $$50 \times 50$$ matrix of (3,420,375) counts Q1--which pass a certain "separability" test--to the total $$50 \times 50$$ matrix of (27,615,417) counts, which we designate as q1, $$$$S=\frac{Q1}{q1}.$$$$ (The Q1 and q1 matrices are given in BivariateCopulaRecursive.)

Since our domain of interest is the unit square, $$[0,1] \times [0,1]$$, we transform the two variables to $$u=\frac{2 i -1}{100}$$ and $$v=\frac{2 j-1}{100}$$, with $$i$$ and $$j$$ running from 1 to 50.

So, the single parameter we need to estimate is $$\theta$$.

Of the 2,500 cells in $$S$$, the number greater than zero is 1,114 (the other cells which we presently ignore are either 0 or indeterminate). The command

w = 0; Do[If[NumberQ[S[[i, j]]] == True && S[[i, j]] > 0, w = w + 1],{i, 1, 50}, {j, 1, 50}];w

gives us for this number of ratios to be incorporated in the analysis, w = 1114.

Now, we do a least-squares minimization

data = Array[0 &, w]; expr = data; diff = 0; w = 0;Do[If[NumberQ[S[[i, j]]] == True && S[[i, j]] != 0, w = w + 1; data[[w]] = S[[i, j]];expr[[w]] = Simplify[AMH]; diff = diff + (data[[w]] - expr[[w]])^2], {i, 1, 50}, {j,1, 50}]; w

to estimate the single AMH parameter $$t$$, on the sum (diff) of the squares of the differences of the 1,111 entries of data and the AMH formula for the entries.

The result of the subsequent minimization

c = NMinimize[diff, t]

is

{14.3606, {t -> -1.90858}}.

Further,

T = t /. c[[2]]; BMH = Simplify[AMH /. t -> T]; R = Array[0 &, {50, 50}]; Do[R[[i, j]] = BMH, {i, 1, 50}, {j, 1, 50}]; v = 50; ListPlot3D[R, DataRange -> {{1/50, v/50}, {1/50, v/50}}]

gives us the $$50 \times 50$$ matrix of predicted AMH values (see the second-to-end ListPlot3D in WikiCopulas). The second ListPlot3D is the difference between this fitted AMH model and the (raw, strictly data-based) first ListPlot3D of Q1/q1 given in BivariateCopulaRepulsive that showed apparent repulsion along the $$x=y$$ line.

Within this framework, we plan to estimate the parameters of the other copula models given in the wiki article, and see if we can achieve least-squares fits superior to the AMH one of 14.3606.

Clayton (15.9567), Frank (14.9727), Gumbel (14.424) and independence (14.7482) models also now included in cloud notebook--but still slightly inferior fit to AMH (14.3606).

However, Joe copula (13.4324) with $$t= 0.493663$$ gives best fit to this point--but still far from terrific (see plot now available).

So, can our analyses also be performed instead of with the use of the wiki function $$C_{\theta}$$, the Mathematica command CopulaDistribution[{"AMH" ]? Presently, I don't see how. They only seem to work for very specific settings.

• Maybe CopulaDistribution[{"AMH", \[Alpha]}, {GammaDistribution[a1, b1, g1, 0], GammaDistribution[a2, b2, g2, 0]}]. I haven't got FindDistributionParameters to work using WeightedData to estimate the parameters (yet).
– JimB
Dec 23 '21 at 5:28
• JimB: Thanks as always for your skilled input. However, I think you are unnecessarily complicating the (conceptually simple?) problem at hand (strongly motivated by the quant-phys "invariance" phenomena observed and of major interest--I'll post some papers in the subsequent comment)--which is to find/define a probability distribution over the unit square that has UNIFORM marginals. Have you looked at the Copula Kernels (11) section under Scope in the CopulaDistribution Help page? The great emphasis is on having UniformDistribution as the desired marginals. Dec 23 '21 at 12:59
• Indicated references from prior comment--arxiv.org/abs/1408.3666--note the phrase "independent of η" in the extended abstract. Also, arxiv.org/abs/1610.01410--with Invariance as the first word of the title. Dec 23 '21 at 13:06
• I'll take a look at these early next week. In the meantime I'll still claim the bivariate distributions of q1 and Q1 do not have uniform distributions. A copula by definition has uniform marginal distributions but is just a means to an end which is to induce some dependence between two random variables having known (or at least estimated) marginal distributions. And it still sounds like it is the ratio of the pdf's of q1 and Q1 that is of interest rather than the ratio of q1 to Q1.
– JimB
Dec 24 '21 at 6:09
• No, no, no--it IS the ratio of Q1 to q1 that is of interest--that was the subject (in a slightly different setting) investigated in arxiv.org/abs/1506.08739. The extended abstract has "Moreover, we study the forms that the separability probabilities take as joint (bivariate) functions of the radii (rA,rB) of the Bloch balls of both single-qubit subsystems. Here, a form of Bloch radii repulsion for separable two-qubit systems emerges in our several analyses". rA and rB are the two variables. The marginal uniiformity/invariance was a major discovery, and this seeks to expand upon that. Dec 24 '21 at 18:25