# 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} \]

*This page was generated using Literate.jl.*