Bayesian Integrals on Toric Varieties
This notebook accompanies the paper “Bayesian Integrals on Toric Varieties” by Michael Borinsky, Anna-Laura Sattelberger, Bernd Sturmfels, and Simon Telen. The code runs in Julia v1.7.1 and uses the packages
Polymake (v0.7.1)
DynamicPolynomials (v0.3.21)
HomotopyContinuation (v2.6.3)
IterTools (v1.4.0)
HCubature (v1.5.0)
The file BITV.jl contains our main functions.
[1]:
include("BITV.jl")
polymake version 4.6
Copyright (c) 1997-2021
Ewgenij Gawrilow, Michael Joswig, and the polymake team
Technische Universität Berlin, Germany
https://polymake.org
This is free software licensed under GPL; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
[1]:
getM2 (generic function with 1 method)
Below, the examples are numbered as in the paper.
Example 4.5
We define a list of numerator and denominator factors. In this example, the numerator is \(f\) and the denominator is \(g\)
[2]:
@polyvar x[1:5]
numers = [2*x[1]^2*x[2]^2*x[3]^3*x[4]*x[5]^3 + 3*x[1]^2*x[2]*x[3]^2*x[4]^2*x[5]^4 + 5*x[1]*x[2]^2*x[3]^5*x[4]*x[5]^2];
denoms = [7*x[1]^3*x[2]^3*x[3]^2*x[5]^3+11*x[1]^3*x[2]*x[4]^2*x[5]^5+13*x[1]*x[3]^3*x[4]^3*x[5]^4+17*x[2]^2*x[3]^7*x[4]*x[5]];
We display the numerator and denominator below
[3]:
numers[1]
[3]:
[4]:
denoms[1]
[4]:
We now define the lists of exponents \(\mu\) and \(\nu\) for the numerator and denominator factors respectively. In this example, we integrate \(f^1/g^1\). We also define the matrix \(V\) of ray generators and we fix a number of samples to use. The main function is \({\tt integrate\_tropical}\). The option “verbatim = true” makes sure that some intermediate output is printed.
[5]:
μ = [1] # exponents of the factors in the numerator
ν = [1] # exponents of the factors in the denominator
V = [1 1 -1 -1 0; 0 -1 -1 1 1] # matrix of ray generators
nsamples = 100000 # number of samples
I = integrate_tropical(V,numers,denoms,x,μ,ν,nsamples; verbatim = true)
1. Sum the Newton polytopes
----------------------------------------------
2. Compute the normal fan of the sum
Σ has 6 rays and 6 maximal cones
----------------------------------------------
3. For each maximal cone, compute the winning exponent
[[-1, -1, 1, 1, 0], [-1, 0, 2, 0, -1], [1, 1, -1, -1, 0], [-2, -1, 3, 1, -1], [0, 2, 2, -2, -2], [1, 0, -2, 0, 1]]
----------------------------------------------
4. Compute the sector integral for all 6 maximal cones
<<<<<<< Itr = 37//4 >>>>>>
I_σ^tr = 1//1
Number of samples for this cone: 10811
Computed the value for sector 1: 0.47833986740530376
I_σ^tr = 2//1
Number of samples for this cone: 21622
Computed the value for sector 2: 0.5061742517431539
I_σ^tr = 3//2
Number of samples for this cone: 16216
Computed the value for sector 3: 0.36310802859909797
I_σ^tr = 1//1
Number of samples for this cone: 10811
Computed the value for sector 4: 0.43697192157606796
I_σ^tr = 1//4
Number of samples for this cone: 2703
Computed the value for sector 5: 0.07664476251139477
I_σ^tr = 7//2
Number of samples for this cone: 37838
Computed the value for sector 6: 1.0127646409048066
----------------------------------------------
5. The integral is 2.874003472739825
[5]:
2.874003472739825
We compute a reference value for the integral using numerical cubature procedures implemented in the package HCubature.jl. The optional input parameter “inttol” specifies the relative accuracy with which the sector integrals are computed. Its default value is set to \(10^{-3}\).
[6]:
I = integrate_cubature(V,numers,denoms,x,μ,ν; verbatim = false, inttol = 1e-6)
[6]:
2.872660369526546
We now use Algorithm 4.3 to draw samples from the posterior distribution. We use 100000 candidate samples.
[7]:
nsamples = 100000
samples, densities = sample_rejection(V,numers,denoms,x,μ,ν,nsamples; verbatim = false);
length(samples)
[7]:
21705
We see that about 21000 samples are accepted, and 79000 were rejected. The bound from Proposition 4.4 predicts acceptance of about \((M_1/M_2)\cdot 100000\) samples. We compute this estimate below:
[8]:
M1 = getM1(numers, denoms, μ, ν)
M2 = getM2(numers, denoms, μ, ν)
M1/M2*100000
[8]:
2916.666666666666666666666666666666666666666666666666666666666666666666666666678
We see that the bound is pessimistic. The actual expected acceptance rate is
which is close to the empirical acceptance rate.
[9]:
Itr = 37//4
(1/M2)*I/Itr*100000
[9]:
21739.051445065753667320776101503823254559491131756756756756756756756756756757
Example 5.3
Example 5.3 computes the area of our pentagon. This requires integrating the toric Hessian determinant \(\det H\) against the canonical form. We computed \(\det H\) using Maple and make use of the fact that the denominator can be written as \(q^3\).
[10]:
@polyvar x[1:5]
numers = [4*x[1]^4*x[2]^4*x[3]^2*x[4]^2*x[5]^3+x[1]^3*x[2]^5*x[3]^5*x[4]*x[5]+4*x[1]^3*x[2]^4*x[3]^4*x[4]^2*x[5]^2+9*x[1]^3*x[2]^3*x[3]^3*x[4]^3*x[5]^3+4*x[1]^3*x[2]^2*x[3]^2*x[4]^4*x[5]^4+8*x[1]^2*x[2]^3*x[3]^5*x[4]^3*x[5]^2+4*x[1]^2*x[2]^2*x[3]^4*x[4]^4*x[5]^3+x[1]*x[2]^3*x[3]^7*x[4]^3*x[5]+x[1]*x[2]*x[3]^5*x[4]^5*x[5]^3];
denoms = [x[1]^2*x[2]^2*x[5]+x[1]*x[2]^2*x[3]^2+x[1]*x[4]^2*x[5]^2+x[2]*x[3]^3*x[4]+x[3]^2*x[4]^2*x[5]];
[11]:
nn = numers[1]
[11]:
[12]:
dd = denoms[1]
[12]:
[13]:
μ = [1] # exponents of the factors in the numerator
ν = [3] # exponents of the factors in the denominator
V = [1 1 -1 -1 0; 0 -1 -1 1 1] # matrix of ray generators
nsamples = 100000 # number of samples
I = integrate_tropical(V,numers,denoms,x,μ,ν,nsamples; verbatim = false)
[13]:
2.5005096864871557
Example 6.2
Example 6.2 considers a conditional independence model (coin model). In particular, we consider rank two tensors. Here is how to define the input:
[14]:
@polyvar x[1:6] # These represent [x[0],x[1],s[0],s[1],t[0],t[1]] from the paper.
d = 2
u = [2,1,2]
numers = [binomial(d,ℓ)*(x[1]*x[3]^ℓ*x[4]^(d-ℓ)*(x[5]+x[6])^d + x[2]*x[5]^ℓ*x[6]^(d-ℓ)*(x[3]+x[4])^d) for ℓ = 0:d]
numers = [numers;prod(x)]
denoms = [x[1]+x[2];x[3]+x[4];x[5]+x[6]]
μ = [u;1]
ν = sum(u)*[1;ones(Int64,d)*d] + 2*ones(Int64,d+1)
V = [1 -1 0 0 0 0; 0 0 1 -1 0 0; 0 0 0 0 1 -1]
nsamples = 100000;
Integrating gives approximately 0.00145. To see the values \(I_{\sigma}^{\rm tr}\), choose “verbatim = true”.
[15]:
I = integrate_tropical(V,numers,denoms,x,μ,ν,nsamples; verbatim = false)
[15]:
0.00142410937323444
We check this using HCubature.jl:
[16]:
I = integrate_cubature(V,numers,denoms,x,μ,ν; verbatim = false, inttol = 1e-4)
[16]:
0.0014538991546372338
Example 6.3
This example considers the linear model and the Wachspress model associated to our pentagon.
The linear model has 5 numerator factors coming from the edges of the pentagon, and one from the toric Hessian. We have \(w = (0,1,3,1,0)\) and set \(\gamma = (5,5,5,5,5)\), but we ignore the \(\gamma\)’s at first. They will be taken into account later.
[17]:
@polyvar x[1:5]
numers = [(2*x[1]*x[2]^2*x[5]+x[2]^2*x[3]^2+x[4]^2*x[5]^2)*x[1];
(2*x[1]^2*x[2]*x[5] + 2*x[1]*x[2]*x[3]^2 + x[3]^3*x[4])*x[2];
(2*x[2]^2*x[1]+3*x[2]*x[3]*x[4]+2*x[4]^2*x[5])*x[3]^2;
(2*x[1]*x[4]*x[5]^2+x[2]*x[3]^3+2*x[3]^2*x[4]*x[5])*x[4];
(x[1]^2*x[2]^2+2*x[1]*x[4]^2*x[5]+x[3]^2*x[4]^2)*x[5];
4*x[1]^4*x[2]^4*x[3]^2*x[4]^2*x[5]^3+x[1]^3*x[2]^5*x[3]^5*x[4]*x[5]+4*x[1]^3*x[2]^4*x[3]^4*x[4]^2*x[5]^2+9*x[1]^3*x[2]^3*x[3]^3*x[4]^3*x[5]^3+4*x[1]^3*x[2]^2*x[3]^2*x[4]^4*x[5]^4+8*x[1]^2*x[2]^3*x[3]^5*x[4]^3*x[5]^2+4*x[1]^2*x[2]^2*x[3]^4*x[4]^4*x[5]^3+x[1]*x[2]^3*x[3]^7*x[4]^3*x[5]+x[1]*x[2]*x[3]^5*x[4]^5*x[5]^3]
[17]:
6-element Vector{Polynomial{true, Int64}}:
2x₁²x₂²x₅ + x₁x₂²x₃² + x₁x₄²x₅²
2x₁²x₂²x₅ + 2x₁x₂²x₃² + x₂x₃³x₄
2x₁x₂²x₃² + 3x₂x₃³x₄ + 2x₃²x₄²x₅
2x₁x₄²x₅² + x₂x₃³x₄ + 2x₃²x₄²x₅
x₁²x₂²x₅ + 2x₁x₄²x₅² + x₃²x₄²x₅
4x₁⁴x₂⁴x₃²x₄²x₅³ + x₁³x₂⁵x₃⁵x₄x₅ + 4x₁³x₂⁴x₃⁴x₄²x₅² + 9x₁³x₂³x₃³x₄³x₅³ + 4x₁³x₂²x₃²x₄⁴x₅⁴ + 8x₁²x₂³x₃⁵x₄³x₅² + 4x₁²x₂²x₃⁴x₄⁴x₅³ + x₁x₂³x₃⁷x₄³x₅ + x₁x₂x₃⁵x₄⁵x₅³
Each of the \(p_i\) has a pole along \(q = 0\), where \(q\) is
[18]:
denoms = [x[1]^2*x[2]^2*x[5]+x[1]*x[2]^2*x[3]^2+x[1]*x[4]^2*x[5]^2+x[2]*x[3]^3*x[4]+x[3]^2*x[4]^2*x[5]]
[18]:
1-element Vector{Polynomial{true, Int64}}:
x₁²x₂²x₅ + x₁x₂²x₃² + x₁x₄²x₅² + x₂x₃³x₄ + x₃²x₄²x₅
Here are the data used in Example 6.3.
[19]:
μ = [20,16,10,15,23,1]
ν = [sum(μ)+2];
V = [1 1 -1 -1 0; 0 -1 -1 1 1]
I = integrate_cubature(V,numers,denoms,x,μ,ν; verbatim = false, inttol = 1e-4)
[19]:
1.2478259257919542
To compute \(I_u\), we should scale by \(5^{-84} \cdot (2/5)\). The first factor comes from \(\gamma\), the second divides by the area of the pentagon.
[20]:
I_u = I*5^(-84)*(2/5)
[20]:
9.654585472476343e-60
We now turn to the Wachspress model. This involves the adjoint polynomial
\(A(y)\) pulls back to the degree zero rational function with numerator
[21]:
@polyvar x[1:5]
numer_A = 8*x[1]^4*x[2]^4*x[5]^2+12*x[1]^3*x[2]^4*x[3]^2*x[5]+20*x[1]^3*x[2]^2*x[4]^2*x[5]^3+4*x[1]^2*x[2]^4*x[3]^4+12*x[1]^2*x[2]^3*x[3]^3*x[4]*x[5]+32*x[1]^2*x[2]^2*x[3]^2*x[4]^2*x[5]^2+8*x[1]^2*x[4]^4*x[5]^4+8*x[1]*x[2]^3*x[3]^5*x[4]+12*x[1]*x[2]^2*x[3]^4*x[4]^2*x[5]+12*x[1]*x[2]*x[3]^3*x[4]^3*x[5]^2+12*x[1]*x[3]^2*x[4]^4*x[5]^3+3*x[2]^2*x[3]^6*x[4]^2+8*x[2]*x[3]^5*x[4]^3*x[5]+4*x[3]^4*x[4]^4*x[5]^2
[21]:
and denominator \(q^2\).
[22]:
numers = [(2*x[1]*x[2]^2*x[5]+x[2]^2*x[3]^2+x[4]^2*x[5]^2)*x[1];
(2*x[1]^2*x[2]*x[5] + 2*x[1]*x[2]*x[3]^2 + x[3]^3*x[4])*x[2];
(2*x[2]^2*x[1]+3*x[2]*x[3]*x[4]+2*x[4]^2*x[5])*x[3]^2;
(2*x[1]*x[4]*x[5]^2+x[2]*x[3]^3+2*x[3]^2*x[4]*x[5])*x[4];
(x[1]^2*x[2]^2+2*x[1]*x[4]^2*x[5]+x[3]^2*x[4]^2)*x[5];
4*x[1]^4*x[2]^4*x[3]^2*x[4]^2*x[5]^3+x[1]^3*x[2]^5*x[3]^5*x[4]*x[5]+4*x[1]^3*x[2]^4*x[3]^4*x[4]^2*x[5]^2+9*x[1]^3*x[2]^3*x[3]^3*x[4]^3*x[5]^3+4*x[1]^3*x[2]^2*x[3]^2*x[4]^4*x[5]^4+8*x[1]^2*x[2]^3*x[3]^5*x[4]^3*x[5]^2+4*x[1]^2*x[2]^2*x[3]^4*x[4]^4*x[5]^3+x[1]*x[2]^3*x[3]^7*x[4]^3*x[5]+x[1]*x[2]*x[3]^5*x[4]^5*x[5]^3;
x[1]^2*x[2]^2*x[5]+x[1]*x[2]^2*x[3]^2+x[1]*x[4]^2*x[5]^2+x[2]*x[3]^3*x[4]+x[3]^2*x[4]^2*x[5]]
[22]:
7-element Vector{Polynomial{true, Int64}}:
2x₁²x₂²x₅ + x₁x₂²x₃² + x₁x₄²x₅²
2x₁²x₂²x₅ + 2x₁x₂²x₃² + x₂x₃³x₄
2x₁x₂²x₃² + 3x₂x₃³x₄ + 2x₃²x₄²x₅
2x₁x₄²x₅² + x₂x₃³x₄ + 2x₃²x₄²x₅
x₁²x₂²x₅ + 2x₁x₄²x₅² + x₃²x₄²x₅
4x₁⁴x₂⁴x₃²x₄²x₅³ + x₁³x₂⁵x₃⁵x₄x₅ + 4x₁³x₂⁴x₃⁴x₄²x₅² + 9x₁³x₂³x₃³x₄³x₅³ + 4x₁³x₂²x₃²x₄⁴x₅⁴ + 8x₁²x₂³x₃⁵x₄³x₅² + 4x₁²x₂²x₃⁴x₄⁴x₅³ + x₁x₂³x₃⁷x₄³x₅ + x₁x₂x₃⁵x₄⁵x₅³
x₁²x₂²x₅ + x₁x₂²x₃² + x₁x₄²x₅² + x₂x₃³x₄ + x₃²x₄²x₅
[23]:
denoms = [numer_A]
[23]:
1-element Vector{Polynomial{true, Int64}}:
8x₁⁴x₂⁴x₅² + 12x₁³x₂⁴x₃²x₅ + 20x₁³x₂²x₄²x₅³ + 4x₁²x₂⁴x₃⁴ + 12x₁²x₂³x₃³x₄x₅ + 32x₁²x₂²x₃²x₄²x₅² + 8x₁²x₄⁴x₅⁴ + 8x₁x₂³x₃⁵x₄ + 12x₁x₂²x₃⁴x₄²x₅ + 12x₁x₂x₃³x₄³x₅² + 12x₁x₃²x₄⁴x₅³ + 3x₂²x₃⁶x₄² + 8x₂x₃⁵x₄³x₅ + 4x₃⁴x₄⁴x₅²
We now insert the data from Example 6.3
[24]:
u = [20,16,10,15,23]
μ = [u...,1,2*84-sum(u)-3]
ν = [84];
V = [1 1 -1 -1 0; 0 -1 -1 1 1];
We compute the integral and scale by \(2^{13} \cdot (2/5)\) to obtain \(I_u\):
[25]:
R = integrate_cubature(V,numers,denoms,x,μ,ν; verbatim = false, inttol = 1e-4)
I_u = R*2^(13)*2/5
[25]:
1.2185066463399776e-66
Example 6.4
We now use our code for Bayesian model decision. As in Example 6.4, we compute a Bayesian factor from two different toric models.
[26]:
@polyvar x[1:5]
c1 = [2,3,5,7,11,13]
numers1 = [c1[1]*x[2]*x[3]^3*x[4]; c1[2]*x[1]*x[2]^2*x[3]^2; c1[3]*x[3]^2*x[4]^2*x[5];
c1[4]*prod(x); c1[5]*x[1]^2*x[2]^2*x[5]; c1[6]*x[1]*x[4]^2*x[5]^2;
4*x[1]^4*x[2]^4*x[3]^2*x[4]^2*x[5]^3+x[1]^3*x[2]^5*x[3]^5*x[4]*x[5]+4*x[1]^3*x[2]^4*x[3]^4*x[4]^2*x[5]^2+9*x[1]^3*x[2]^3*x[3]^3*x[4]^3*x[5]^3+4*x[1]^3*x[2]^2*x[3]^2*x[4]^4*x[5]^4+8*x[1]^2*x[2]^3*x[3]^5*x[4]^3*x[5]^2+4*x[1]^2*x[2]^2*x[3]^4*x[4]^4*x[5]^3+x[1]*x[2]^3*x[3]^7*x[4]^3*x[5]+x[1]*x[2]*x[3]^5*x[4]^5*x[5]^3];
denoms1 = [sum(numers1[1:6]);x[1]^2*x[2]^2*x[5]+x[1]*x[2]^2*x[3]^2+x[1]*x[4]^2*x[5]^2+x[2]*x[3]^3*x[4]+x[3]^2*x[4]^2*x[5]]
c2 = [32,16,8,4,2,1]
numers2 = [c2[1]*x[2]*x[3]^3*x[4]; c2[2]*x[1]*x[2]^2*x[3]^2; c2[3]*x[3]^2*x[4]^2*x[5];
c2[4]*prod(x); c2[5]*x[1]^2*x[2]^2*x[5]; c2[6]*x[1]*x[4]^2*x[5]^2;
4*x[1]^4*x[2]^4*x[3]^2*x[4]^2*x[5]^3+x[1]^3*x[2]^5*x[3]^5*x[4]*x[5]+4*x[1]^3*x[2]^4*x[3]^4*x[4]^2*x[5]^2+9*x[1]^3*x[2]^3*x[3]^3*x[4]^3*x[5]^3+4*x[1]^3*x[2]^2*x[3]^2*x[4]^4*x[5]^4+8*x[1]^2*x[2]^3*x[3]^5*x[4]^3*x[5]^2+4*x[1]^2*x[2]^2*x[3]^4*x[4]^4*x[5]^3+x[1]*x[2]^3*x[3]^7*x[4]^3*x[5]+x[1]*x[2]*x[3]^5*x[4]^5*x[5]^3];
denoms2 = [sum(numers2[1:6]);x[1]^2*x[2]^2*x[5]+x[1]*x[2]^2*x[3]^2+x[1]*x[4]^2*x[5]^2+x[2]*x[3]^3*x[4]+x[3]^2*x[4]^2*x[5]];
Here are the data used in Example 6.4:
[ ]:
u = [1,2,4,8,16,32]
μ = [u...,1];
ν = [sum(u),3];
V = [1 1 -1 -1 0; 0 -1 -1 1 1];
nsamples = 100000
I1 = (2/5)*integrate_cubature(V,numers1,denoms1,x,μ,ν; verbatim=false, inttol = 1e-5)
I2 = (2/5)*integrate_cubature(V,numers2,denoms2,x,μ,ν; verbatim=false, inttol = 1e-5)
[I1 I2]
The Bayes factor is the ratio of these marginal likelihood integrals.
[ ]:
I1/I2