Uncertainty Propagation

Published

April 9, 2024

Here we will see a typical example of uncertainty propagation through a nonlinear function. Assume that we want to propagate the uncertainty in \(x\) through the function $ y = f(x)$ where \(f\) is nonlinear. To begin with we will assume that \(x\) is normally distributed with mean \(\mu_x\) and standard deviation \(\sigma_x\). Remember that we have already argued out that for most of the real life variables, the normal distribution is a good approximation of the true distribution because of the central limit theorem.

As a simple example, let \(f(x) = \sin(x)\).

The first step is to sample the input variable \(x\) from the normal distribution as follows. We need to first decide on the number of sampling points \(N\). This would typically depend on the complexity of the function \(f\) and the desired accuracy of the output. For this example, we will choose \(N = 1000000\).

Code 
# Define the number of sampling points
N = 1000000

# Input the uncertainty in x
μx = 0
σx = pi/16

# Sample the input variable
x̃ = μx .+ σx * randn(N)

Let us check the mean and standard deviation of the sampled input variable.

Code 
μx̃ = mean(x̃)
σx̃ = std(x̃)

Clearly the calculated mean 0.00022385939277158012 is close to the input mean 0 and the standard deviation 0.19635724242675154 is close to the input standard deviation 0.19634954084936207. May be not close enough. If we increase the number of sampling points, the mean and standard deviation predictions can be improved.

Next we will evaluate the function \(f(x)\) at the sampled input points.

Code 
# Define the sine function
function f(x)
    return sin(x)
end

# Monte Carlo Simulations
ỹ = f.(x̃)

# Define the range of the mean and standard deviation of the perturbed output variable for each simulation
μy = mean(ỹ)
σy = std(ỹ)

# Plot x̃ histogram and ỹ histogram in two subplots
p1 = histogram(x̃, label="x̃", alpha=0.5)
title!("μx = $(round(μx, digits=3)), σx = $(round(σx, digits=3))")

p2 = histogram(ỹ, label="ỹ", alpha=0.5)
title!("μy = $(round(μy, digits=3)), σy = $(round(σy, digits=3))")

plot(p1, p2, layout = (2, 1))
# Define the derivative of the sine function
function df(x)
    return cos(x)
end

# Investigating the convergence rate of Monte Carlo simulations
# Define the number of simulations
M = 14

μx = 0.0
σx = pi/16

# Initialize arrays to store the mean and standard deviation of the perturbed output variable
μy = zeros(M)
σy = zeros(M)

# Perform Monte Carlo simulations
for i in 1:M
    x̃ = μx .+ σx * randn(2^i)
    ỹ = f.(x̃)
    μy[i] = mean(ỹ)
    σy[i] = std(ỹ)
end

# Plot the convergence of the mean and standard deviation of the perturbed output variable

p1 = plot(1:M, μy, label="Mean", xlabel="Number of Simulations", ylabel="Value", title="Convergence of Mean")
p2 = plot(1:M, σy, label="Standard Deviation", xlabel="Number of Simulations", ylabel="Value", title="Convergence of Standard Deviation")

plot(p1, p2, layout = (2, 1))


# Investigating the convergence rate of Monte Carlo simulations
# about μx = 0.0

# Define the number of simulations
M = 6

μx = 0.0
σx = pi/16

# Initialize arrays to store the mean and standard deviation of the perturbed output variable
μy = zeros(M)
σy = zeros(M)

# Perform Monte Carlo simulations
for i in 1:M
    x̃ = μx .+ σx * randn(10^i)
    ỹ = f.(x̃)
    μy[i] = mean(ỹ)
    σy[i] = std(ỹ)
end

# Plot the convergence of the mean and standard deviation of the perturbed output variable

p1 = scatter(1:M, μy, label="Mean", xlabel="Number of Simulations", ylabel="Value", title="Convergence of Mean")
p2 = scatter(1:M, σy, label="Standard Deviation", xlabel="Number of Simulations", ylabel="Value", title="Convergence of Standard Deviation")

plot(p1, p2, layout = (2, 1))

# Plot the error in the mean and standard deviation of the perturbed output variable in log-log scale
# We also plot 1/sqrt(N) for reference as a bold black line

error_mean = abs.(μy .- sin(μx))
error_std = abs.(σy .- abs(cos(μx) * σx))

p1 = scatter(1:M, error_mean, label="Mean Error", xlabel="Number of Simulations", ylabel="Error", title="Error in Mean")
plot!(1:M, 1 ./ sqrt.(10 .^ (1:M)), label="1/sqrt(N)", color=:black, linewidth=2)

p2 = scatter(1:M, error_std, label="Standard Deviation Error", xlabel="Number of Simulations", ylabel="Error", title="Error in Standard Deviation")
plot!(1:M, 1 ./ sqrt.(10 .^ (1:M)), label="1/sqrt(N)", color=:black, linewidth=2)

plot(p1, p2, layout = (2, 1), yscale=:log10)