Integration of pdf of normal distribution

[shmarovfedor]

Just a reminder. I've already reported about this issue earlier. I was integrating a PDF of standard normal distribution on the interval [-10, 10] and the output was:

SAT with the following box:
         time_0 : [19.99999761581421, 20.00000262260438];
            P_t : [0.9999981259955913, 0.9999998047304542];
            P_0 : [0, 0];
            x_t : [10, 10];
            x_0 : [-10, -10]

When I changed the bounds of the integration to [-20, 20] the output was:

SAT with the following box:
         time_0 : [39.99999999999986, 40];
            P_t : [0, 0.0009765625];
            P_0 : [0, 0];
            x_t : [20, 20];
            x_0 : [-20, -20]

But the value of the integral should be converging to 1. Here is .smt2 file:

(set-logic QF_NRA_ODE)
(declare-fun x () Real)
(declare-fun P () Real)

(declare-fun x_0 () Real)
(declare-fun x_t () Real)

(declare-fun P_0 () Real)
(declare-fun P_t () Real)

(declare-fun time_0 () Real)

(define-ode flow_1 (
(= d/dt[x] 1.0)
(= d/dt[P] (* (/ 1.0 (^ (* 2.0 3.14159265359) 0.5)) (exp (/ (- 0.0 (^ (- x 0.0) 2.0)) 2.0)))) 
))

(assert (<= -100.0 x_0))
(assert (<= x_0 100.0))

(assert (<= -100.0 x_t))
(assert (<= x_t 100.0))

(assert (<= 0.0 P_0))
(assert (<= P_0 1.0))

(assert (<= 0.0 P_t))
(assert (<= P_t 1.0))

(assert (<= 0.0 time_0))
(assert (<= time_0 40.0))

(assert (and
                (= x_0 -20.0)
                (= P_0 0.0)
                (= x_t 20.0)
                (>= P_t 0.0)
                (= [x_t P_t] (integral 0. time_0 [x_0 P_0] flow_1))
))

(check-sat)
(exit)

Also, when I changed the standard deviation in the pdf (from 1 to 0.01) and integrated it on the interval [-10, 10], the output was:

SAT with the following box:
         time_0 : [19.99999999999999, 20.00000000000001];
            P_t : [0, 0.0009765625];
            P_0 : [0, 0];
            x_t : [10, 10];
            x_0 : [-10, -10]

[soonhokong]

Hi Fedor, thanks for the report. Yes, it's on my todo list but it's nice to have it filed here. I'll work on it this weekend.

[soonhokong]

@shmarovfedor, @pzuliani, @scungao

Problem

I'm working on a solution for the problem, but I think you may want to know what's happening here. As fedor said, when the interval is small (i.e. [-10, 10]), CAPD picks up a small time step for the ODE integration and ends up with a precise result (i.e. P_t = [0.9999981259955913, 0.9999998047304542]). Please note that the interval is tight (width <= 10^-6).

When the interval is large, CAPD starts behaving differently. For example, when an interval is [-20, 20], CAPD picks a large time-step such as [0, 1] for the ODE integration and ends up with an imprecise integration result (P_t =[0, 1]). It means that P can be any value between 0 and 1, which is obvious information. Then ICP solver branches the domain [0, 1] into a smaller region [0, 0.0009..](<= delta), and returns it.

Thoughts

I first enforced CAPD to use a small ODE time-step (i.e. 0.001) instead of letting CAPD pick one for itself. After some iterations, unfortunately, CAPD generated an exception and failed to do the ODE integration.

I'm trying to replace CAPD-3.0 in dReal with the new version CAPD-4.0 hoping that it might give us a better result. However, I think this kind of behavior is unavoidable at some level. We need to think about 1) how to detect this kind of imprecise pruning operation, and 2) a back-up plan.

[scungao]

This is a case where CAPD violates the specification of a delta-regular interval extension (see our SAT modulo ODE paper in FMCAD13).

First, let's contact the authors of CAPD about this.

In general, if we see that a pruning operator is not working (this is one scenario of failed pruning), we should probably notify the user and return the last interval assignment, before ICP branches on it and picks a random value.

You mentioned that it works for [-18,18] right? What exception is it throwing when you set time step?

[pzuliani]

Thanks for looking into this!

BTW, which condition of delta-reg extension is actually being violated?

[scungao]

In Definition 3.3, we ask the size of the interval extension to be proportional to the size of the input interval -- for the case right now, it is always too big and does not decrease w.r.t. smaller inputs.

[scungao]

What's the status of this one? @soonhokong

[ftang921]

Have you solved this problem? Thank you very much

[soonhokong]

@pzuliani, @shmarovfedor, @scungao, and @arwintang, I think the problem is on CAPD's automatic step control. I will send an email to the authors of CAPD4, and let you know when I get a feedback from them.

Please see https://gist.github.com/soonhokong/f5cb88085621f37d4a88 for details.

[soonhokong]

I just sent an email. One way to workaround the issue is to specify the ODE time-step manually. I'll re-implement this to dreal3. See https://github.com/dreal/dreal3/issues/123 to follow the progress.

[soonhokong]

@arwintang, FYI, we're working on the new release of dReal, dReal3. Please try it and let us know if you find problems.

[pzuliani]

Great, many thanks @soonhokong !

[ftang921]

@soonhokong thank you very much. I tried dreal3 with this example, but the result is unsat.

[soonhokong]

@arwintang, there was a problem in the original fomula (I edited it now). Before I edited, it was

(define-ode flow_1 (
(= d/dt[P] (* (/ 1.0 (^ (* 2.0 3.14159265359) 0.5)) (exp (/ (- 0.0 (^ (- x 0.0) 2.0)) 2.0)))) 
(= d/dt[x] 1.0)
))

Now, it's

(define-ode flow_1 (
(= d/dt[x] 1.0)
(= d/dt[P] (* (/ 1.0 (^ (* 2.0 3.14159265359) 0.5)) (exp (/ (- 0.0 (^ (- x 0.0) 2.0)) 2.0)))) 
))

The variable ordering matters in dReal3. For the updated formula, dReal3 returns delta-sat, but the result is not accurate because of the reason I mentioned in gist. I'll share the feedbacks from the CAPD team.

[ftang921]

@soonhokong thank you very much!

[soonhokong]

Got a response from CAPD team saying that they will work on the issue. They also suggested to handle Pi(3.141592...) specially. We can introduce a predefined constant "pi" and let users of dreal use it instead of using an approximated constant. When we formulate CAPD/IBEX problem, we can use predefined constants in numeric solvers. Possibly this can increase precision. I'll implement this soon.

[soonhokong]

1. For this particular example, using symbolic "pi" does not increase precision. Anyway, I think it's a good thing to have. I'll deploy the feature soon.
2. I got another email from CAPD team this morning. They said they found the source of the problem and are working on a fix.

[soonhokong]

I just added support for --ode-step option in dReal3. Using it certainly helps:

$ ./dReal p.smt2 --model --ode-step 0.1
Solution:
P_0 : [0, 1] = [0, 0]
P_t : [0, 1] = [0.9999999999999381, 1]
time_0 : [0, 40] = [39.99999999999872, 40]
x_0 : [-100, 100] = [-20, -20]
x_t : [-100, 100] = [20, 20]

delta-sat with delta = 0.00100000000000000

$ ./dReal p.smt2 --model --ode-step 0.01
Solution:
P_0 : [0, 1] = [0, 0]
P_t : [0, 1] = [0.9999999999997366, 1]
time_0 : [0, 40] = [39.99999999999048, 40]
x_0 : [-100, 100] = [-20, -20]
x_t : [-100, 100] = [20, 20]

delta-sat with delta = 0.00100000000000000

For now, this is all I can do. Let's wait for CAPD team fixing the problem so that we can handle this case using automatic time-step control.

[pzuliani]

@soonhokong excellent!