# Bounds in Probability

Adapted from: SOSTOOLS' SOSDEMO8 (See Section 4.8 of SOSTOOLS User's Manual)

The probability adds up to one.

μ0 = 1
1

The mean is one.

μ1  = 1
1

The standard deviation is 1/2.

σ = 1/2
0.5

The second moment E(x^2) is:

μ2 = σ^2 + μ1^2
1.25

We define the moments as follows:

using DynamicPolynomials
@polyvar x
monos = [1, x, x^2]
using SumOfSquares
μ = measure([μ0, μ1, μ2], monos)
Measure{Float64, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}([1.0, 1.0, 1.25], DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[1, x, x²])

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver);

We create a polynomial with the monomials in monos and JuMP decision variables as coefficients as follows:

@variable(model, poly, Poly(monos))
$$$({\_}_{1}) + ({\_}_{2})x + ({\_}_{3})x^{2}$$$

Nonnegative on the support:

K = @set 0 <= x && x <= 5
con_ref = @constraint(model, poly >= 0, domain = K)
$$$({\_}_{1}) + ({\_}_{2})x + ({\_}_{3})x^{2} \text{ is SOS}$$$

Greater than one on the event:

@constraint(model, poly >= 1, domain = (@set 4 <= x && x <= 5))
$$$({\_}_{1} - 1) + ({\_}_{2})x + ({\_}_{3})x^{2} \text{ is SOS}$$$

The bound (we use LinearAlgebra for the ⋅ syntax for the scalar product):

using LinearAlgebra
@objective(model, Min, poly ⋅ μ)

${_}{1} + {\}{2} + 1.25 {\}_{3}$

We verify that we found a feasible solution:

optimize!(model)
primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1

The objective value is 1/37:

objective_value(model)
0.027027027945194515

The solution is (12x-11)^2 / 37^2:

value(poly) * 37^2
$$$120.9999573610103 - 263.99994311062636x + 143.99998960526992x^{2}$$$