# Monte Carlo Uncertainty Analysis

Monte Carlo *Uncertainty Analysis* is a technique used to investigate how the inherent uncertainty associated with *measured input parameters* affects the output of a mathematical model.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import scipy.integrate as spi
import scipy.optimize as spo

## Overview of Uncertainty

### Sources of Uncertainty
Input parameters in models are often *measured quantities* (e.g., mass, time, distance). Even with a perfectly executed measurement, there is always an *uncertainty* associated with the result, often stemming from the limitations or inaccuracy of the tools used.

For example, a volumetric flask might have a volume measurement of $100 \pm 0.08$ mL, meaning the measurement is truly a range of values.

### The Role of Random Error
When many measurements are taken, the distribution of results often approaches a *normal curve* due to *random errors*. Random errors cause both over- and underestimates of the true value, which tend to cancel out on average. Monte Carlo leverages this natural distribution to model uncertainty.

## Monte Carlo Uncertainty Procedure
The general strategy is to incorporate the measured uncertainty $c$ into the model by simulating it as **random noise**:

1. **Define the Random Variable:** For a measured quantity $x$ with uncertainty $\pm c$, define a random variable $\hat{x} = x + \varepsilon$.
2. **Define the Error $\varepsilon$:** $\varepsilon$ is a *normally distributed* random variable with a *mean of $0$* and a *standard deviation of $c$* (the uncertainty level).
* Choosing a **zero mean** assumes the system behaves, on average, like the measured value (i.e., no systematic errors).
* Choosing the **standard deviation** equal to $c$ allows the variation to be correctly injected into the parameter's value.
3. **Simulate and Solve:** Solve the model multiple times using the randomly generated $\hat{x}$ value for the input parameter in each simulation.
4. **Visualize:** Plot the distribution of the final results to assess the output uncertainty.

## Example: Racecar Velocity
The velocity formula $v = d/t$ is calculated using uncertain measurements for distance $d$ and time $t$.

| Quantity | Value | Uncertainty (±c) |
| :--- | :---: | :---: |
| Distance ($d$) | $54 \text{ m}$ | $1 \text{ m}$ |
| Time ($t$) | $2 \text{ s}$ | $0.1 \text{ s}$ |

The calculated velocity without uncertainty is $v = 54/2 = 27 \text{ m/s}$.

### Monte Carlo Simulation
We generate $N=10,000$ samples for $d$ and $t$ from Normal distributions centered at their measured values, with a standard deviation equal to the uncertainty.

# Apply monte carlo to see how uncertainty affects the velocity 
d = 54; t = 2
d_err = 1; t_err = 0.1 
N = 10000 

# Generate Monte Carlo samples using N(mean=value, scale=error)
d_mc = d + np.random.normal(loc=0.0, scale=d_err,size=N)
t_mc = t + np.random.normal(loc=0.0, scale=t_err,size=N)

v_mc = d_mc/t_mc

plt.hist(v_mc,bins = 20)
plt.axvline(x = v, color = 'r') #vertical line for exact value
plt.xlabel('velocity')
plt.ylabel('frequency')
plt.show()


#### Analysis
The resulting distribution of velocities is centered at $27 \text{ m/s}$. The spread of the distribution quantifies the uncertainty in the final velocity calculation: a measured velocity could range from approximately $23 \text{ m/s}$ to $33 \text{ m/s}$.

## Example: Energy Balance Model
We apply Monte Carlo uncertainty analysis to the *Global Energy Balance Model* to determine how uncertainties in measured physical constants affect the calculated temperature solution $T(t)$.

### Model and Uncertain Parameters

#### The Energy Balance Model
The deterministic differential equation for Earth's temperature ($T$) is:

$$
k \frac{dT}{dt} = \pi R^2 Q(1-\alpha(T)) - 4 \pi R^2 \sigma \varepsilon T^4
$$

| Variable | Definition |
| :---: | :--- |
| $k$ | Heat capacity of the Earth (J/K) |
| $Q$ | Solar energy reaching Earth's surface (s$^{-1}$ m$^{-2}$) |
| $R$ | Radius of the Earth (m) |
| $\sigma$ | Stefan-Boltzmann constant (W/(m$^2$K$^4$)) |
| $\alpha(T)$ | Albedo (unitless, dependent on T or constant) |
| $\varepsilon$ | Emissivity (unitless, constant $0.622$) |

#### Uncertain Constant Values
The measured constants are defined with a central value and an associated uncertainty ($\pm c$):

| Constant | Value | Uncertainty (±c) |
| :---: | :---: | :---: |
| $k$ | $1.0 \times 10^{23} \text{ J/K}$ | $\pm 0.22 \times 10^{9} \text{ J/K}$ |
| $R$ | $6.3781 \times 10^6 \text{ m}$ | $\pm 10 \text{ m}$ |
| $Q$ | $1365 \text{ W}/(\text{s}\cdot\text{m}^2)$ | $\pm 27.3 \text{ W}/(\text{s}\cdot\text{m}^2)$ |
| $\sigma$ | $5.6704 \times 10^{-8} \text{ W}/(\text{m}^2\text{K}^4)$ | $\pm 0.000021 \times 10^{-8} \text{ W}/(\text{m}^2\text{K}^4)$ |

#### Albedo Functions
The albedo can be modeled as constant

$$
\alpha_1(T) = 0.3
$$

or temperature-dependent

$$
\alpha_2(T)=
\begin{cases}
0.7, & T \leq 247 \text{ K} \\
3.52296-0.011429T, & 247 < T < 282 \text{ K} \\
0.3, & T \geq 282 \text{ K}
\end{cases}
$$

### Monte Carlo Implementation: Constant Albedo $\alpha_1$
For the first analysis, we use the simple constant albedo $\alpha_1(T)=0.3$ and inject *Normal random noise* (mean 0, standard deviation = uncertainty) into the four uncertain constants ($k, R, Q, \sigma$). We then solve the ODE for each set of noisy parameters to get a distribution of solutions.

# --- Deterministic Constants ---
Q_det = 1365 
R_det = 6.3781e6 
k_det = 1e23 
sigma_det = 5.670374419e-8 
eps = 0.622
alpha = 0.3

# --- Uncertainty Values ---
k_err = 0.22e9
R_err = 10
Q_err = 27.3
sigma_err = 0.000021e-8

# --- Monte Carlo Setup ---
N = 1000 
t = np.linspace(0,1e9,100)
T0 = 267
T_mc_results = []

# --- Monte Carlo Simulation ---
for _ in range(N):
    # 1. Define Noisy Parameters (Value + N(0, error))
    Q_mc = Q_det + np.random.normal(loc=0, scale=Q_err)
    R_mc = R_det + np.random.normal(loc=0, scale=R_err)
    k_mc = k_det + np.random.normal(loc=0, scale=k_err)
    sigma_mc = sigma_det + np.random.normal(loc=0, scale=sigma_err)
    
    # 2. Define ODE function with Noisy Parameters
    f = lambda T, t: (np.pi * R_mc**2 * Q_mc * (1 - alpha) - 
                      4 * np.pi * R_mc**2 * eps * sigma_mc * T**4) / k_mc
    
    # 3. Solve the ODE
    T_mc = spi.odeint(f, T0, t)
    T_mc_results.append(T_mc[:, 0])

# --- Deterministic Solution (for comparison) ---
f_det = lambda T,t: (np.pi * R_det**2 * Q_det * (1 - alpha) - 
                     4 * np.pi * R_det**2 * eps * sigma_det * T**4) / k_det
T_determin = spi.odeint(f_det, T0, t)

# --- Plotting Results ---
plt.figure(figsize=(10, 6))
for T_series in T_mc_results:
    plt.plot(t, T_series, color='gray', alpha=0.2)
plt.plot(t, T_determin, color='red', linewidth=2, label='Deterministic Solution')
plt.title('Energy Balance Model Solutions with Parameter Uncertainty ($\ \\alpha_1(T)=0.3$)')
plt.xlabel('Time (s)')
plt.ylabel('Temperature (K)')
plt.legend()
plt.grid(True)
plt.show()


### Monte Carlo Implementation: Temperature-Dependent Albedo $\alpha_2$
In this more complex analysis, we use the *temperature-dependent albedo* $\alpha_2(T)$ and focus the uncertainty on the *initial temperature* $T_0$, simulating $T_0 = 267 \pm 1 \text{ K}$. This helps determine how a small initial uncertainty affects the long-term solution, especially around the critical ice/water boundaries defined by $\alpha_2(T)$.

# --- Deterministic Constants ---
Q = 1365 
R = 6.3781e6 
k = 1e23 
sigma = 5.670374419e-8 
eps = 0.622

# --- Albedo 2 Function ---
def alpha2(T):
    if T <= 247:
        return 0.7
    elif 247 < T < 282:
        return 3.52296 - 0.011429 * T
    else:
        return 0.3

# Vectorize the Albedo function for use in ODE
alpha_vec = np.vectorize(alpha2)

# --- Monte Carlo Setup ---
N = 1000 
t = np.linspace(0, 1e9, 100)
T0_det = 267
T0_err = 1 # Uncertainty of +/- 1 K
T_mc_results2 = []

# --- Monte Carlo Simulation ---
for _ in range(N):
    # 1. Define Noisy Initial Temperature
    T0_mc = T0_det + np.random.normal(loc=0, scale=T0_err)
    
    # 2. Define ODE function with alpha2
    f1 = lambda T, t: (np.pi * R**2 * Q * (1 - alpha_vec(T)) - 
                       4 * np.pi * R**2 * eps * sigma * T**4) / k
    
    # 3. Solve the ODE
    T_mc = spi.odeint(f1, T0_mc, t)
    T_mc_results2.append(T_mc[:, 0])

# --- Plotting Results ---
plt.figure(figsize=(10, 6))
for T_series in T_mc_results2:
    plt.plot(t, T_series, color='gray', alpha=0.2)
plt.title('Energy Balance Model Solutions with Initial Temperature Uncertainty ($\ \\alpha_2(T)$)')
plt.xlabel('Time (s)')
plt.ylabel('Temperature (K)')
plt.grid(True)
plt.show()


### Analysis
By comparing the results from the $\alpha_1(T)$ and $\alpha_2(T)$ models, we can see that models with a temperature-dependent albedo are generally more sensitive to initial conditions. A small initial error (like the $\pm 1 \text{ K}$) can lead to a wider range of final equilibrium temperatures because the system can be pushed past a tipping point (e.g., above or below the $282 \text{ K}$ or $247 \text{ K}$ thresholds in $\alpha_2(T)$). This highlights the necessity of using Monte Carlo to assess sensitivity in complex, nonlinear models.