```python
import numpy as np
import sympy as sp
```

First, let's derive the inverse X matrix that maps from fluxes to "a" coefficients:


```python
X = np.vander([-1, 0, 1], 3)
np.linalg.inv(X)
```




    array([[ 0.5, -1. ,  0.5],
           [-0.5,  0. ,  0.5],
           [ 0. ,  1. ,  0. ]])



Then, let's propagate the uncertainties:


```python
np.linalg.inv(np.dot(X.T, X))
```




    array([[ 1.5,  0. , -1. ],
           [ 0. ,  0.5,  0. ],
           [-1. ,  0. ,  1. ]])



This means that a_1 and a_2 are independent with variances 0.5 and 1.5 respectively.


```python
f_minus, f_max, f_plus = sp.symbols("f_minus f_max f_plus")
f = sp.Matrix([[f_minus], [f_max], [f_plus]])
X_s = sp.Matrix(X)
A = X_s**-1 * f
a2 = A[0, 0]
a1 = A[1, 0]
a0 = A[2, 0]
print("a0 =", a0)
print("a1 =", a1)
print("a2 =", a2)
```

    a0 = 1.0*f_max
    a1 = -0.5*f_minus + 0.5*f_plus
    a2 = -1.0*f_max + 0.5*f_minus + 0.5*f_plus



```python
x0 = sp.symbols("x0")
x_max = x0 - a1 / (2*a2)
print("x_max =", x_max)
```

    x_max = x0 - (-0.5*f_minus + 0.5*f_plus)/(-2.0*f_max + f_minus + f_plus)



```python
y_max = a0 + a1 * (x_max - x0) + a2 * (x_max - x0)**2
print("y_max =", sp.simplify(y_max))
```

    y_max = 1.0*(4.0*f_max**2 - 2.0*f_max*f_minus - 2.0*f_max*f_plus + 0.25*f_minus**2 - 0.5*f_minus*f_plus + 0.25*f_plus**2)/(4.0*f_max - 2.0*f_minus - 2.0*f_plus)



```python
df_minus, df_max, df_plus = sp.symbols("df_minus df_max df_plus")
inv_Sigma = sp.Matrix([[1/df_minus**2, 0, 0], [0, 1/df_max**2, 0], [0, 0, 1/df_plus**2]])
Sigma_a = sp.simplify((X_s.T * inv_Sigma * X_s)**-1)

sig2_a1 = Sigma_a[0, 0]
sig2_a2 = Sigma_a[1, 1]
sig2_a0 = Sigma_a[2, 2]

sig_a1a2 = Sigma_a[0, 1]
sig_a0a2 = Sigma_a[0, 2]
sig_a0a1 = Sigma_a[1, 2]

print("sig2_a0 =", sig2_a0)
print("sig2_a1 =", sig2_a1)
print("sig2_a2 =", sig2_a2)

print("sig_a1a2 =", sig_a1a2)
print("sig_a0a2 =", sig_a0a2)
print("sig_a0a1 =", sig_a0a1)
```

    sig2_a0 = df_max**2
    sig2_a1 = df_max**2 + df_minus**2/4 + df_plus**2/4
    sig2_a2 = df_minus**2/4 + df_plus**2/4
    sig_a1a2 = -df_minus**2/4 + df_plus**2/4
    sig_a0a2 = -df_max**2
    sig_a0a1 = 0



```python
a0_, a1_, a2_, df = sp.symbols("a0_ a1_ a2_ df")
f = a1_ / (2*a2_)

sig2_f = sp.diff(f, a1_)*sp.diff(f, a2_)*sig_a1a2
sig2_f += sp.diff(f, a1_)**2*sig2_a1
sig2_f += sp.diff(f, a2_)**2*sig2_a2
print("x_max_var =", sp.simplify(sig2_f).evalf())

subs = [(df_minus, df), (df_plus, df), (df_max, df)]
sp.simplify(sig2_f.subs(subs))
```

    x_max_var = 0.0625*(a1_**2*(df_minus**2 + df_plus**2) + a1_*a2_*(df_minus**2 - df_plus**2) + a2_**2*(4.0*df_max**2 + df_minus**2 + df_plus**2))/a2_**4





    df**2*(a1_**2 + 3*a2_**2)/(8*a2_**4)




```python
tmp = -a1_/(2*a2_)
y_max = a0_ + a1_ * tmp + a2_ * tmp**2
print("y_max =", y_max)

sig2_y = sp.diff(y_max, a0_)*sp.diff(y_max, a1_)*sig_a0a1
sig2_y += sp.diff(y_max, a0_)*sp.diff(y_max, a2_)*sig_a0a2
sig2_y += sp.diff(y_max, a1_)*sp.diff(y_max, a2_)*sig_a1a2
sig2_y += sp.diff(y_max, a0_)**2*sig2_a0
sig2_y += sp.diff(y_max, a1_)**2*sig2_a1
sig2_y += sp.diff(y_max, a2_)**2*sig2_a2

print("sig2_y =", sp.simplify(sig2_y).evalf())
```

    y_max = a0_ - a1_**2/(4*a2_)
    sig2_y = 0.015625*(a1_**4*(df_minus**2 + df_plus**2) + 2.0*a1_**3*a2_*(df_minus**2 - df_plus**2) + 4.0*a1_**2*a2_**2*(df_minus**2 + df_plus**2) + 64.0*a2_**4*df_max**2)/a2_**4



```python
print("sig2_y =", sp.simplify(sig2_y.subs(subs)).evalf())
```

    sig2_y = df**2*(0.03125*a1_**4 + 0.125*a1_**2*a2_**2 + a2_**4)/a2_**4



```python

```

04be81bbad3ddc0a7a0447de4dc31603e02bc4c4
ipynb
Jupyter Notebook
docs/_static/notebooks/math.ipynb
iancze/bettermoments
d97bd4ae8d11a3670ca4ee97bd658c97455bcc49
MIT
2018-09-21T16:37:08.000Z
2022-03-14T08:34:35.000Z
docs/_static/notebooks/math.ipynb
iancze/bettermoments
d97bd4ae8d11a3670ca4ee97bd658c97455bcc49
MIT
2018-09-13T11:42:41.000Z
2021-08-30T10:59:45.000Z
docs/_static/notebooks/math.ipynb
iancze/bettermoments
d97bd4ae8d11a3670ca4ee97bd658c97455bcc49
MIT
2018-09-17T17:53:35.000Z
2022-01-28T21:54:51.000Z
Qwen/Qwen-72B
1. YES
2. YES