Hi there,
I am trying to plot a probability distribution (PDF) which can only be calculated using integration. Essentially I’m trying to plot the function:
P(y) = \int_{-\infty}^{\infty} P(y|x) P(x) dx
Where P(y|x) is the PDF of y for a given value of x, and P(x) is the PDF of x. I’ve attached the full functional form of all variables as a screenshot. But essentially, both P(y|x) and P(x) are gaussian distributions, however I’m having difficulty defining explicitly what each distribution is a function of.
I have no difficulty plotting the integrated form of just P(x) (which results in an erf(x) cumulative distribution function), but it seems to be running into issues when I combine multiple functions together.
Here is the code that I’m using, along with the set of equations"
function calculate_simple_integral(Data, Obs_Data, slope, intercept, σ_b)
# Extract data and parameters
ΔT_data = Data.ΔT # x data points
N = length(ΔT_data) # Number of x data points
ΔT_Mean = mean(ΔT_data) # Mean of the x data points
σ_ΔT = std(ΔT_data) # `Standard deviation data points
ΔT_Obs_Mean = Obs_Data[1] # Mean of Observational x data points
σ_ΔT_Obs = Obs_Data[2] # Standard deviation of Observational x data points
@syms x
@syms y
# Define f(x)
f(x) = slope * x + intercept
# Define Prediction Error of f(x), \sigma_f(x)
σ_f(x) = σ_b * σ_ΔT * sqrt(1 + N + ((x - ΔT_Mean)^2 / (σ_ΔT^2)))
# Define P(x)
P(x) = (1 / sqrt(2 * π * σ_ΔT_Obs^2)) * exp(-((x - ΔT_Obs_Mean)^2) / (2 * σ_ΔT_Obs^2))
# Define P{y | x}
P_given_y_x(y, x) = (1 / sqrt(2 * π * σ_f(x)^2)) * exp(-((y - f(x))^2) / (2 * σ_f(x)^2))
# Define the integrand: P{y | x} * P(x)
integrand(x, y) = P_given_y_x(y, x) * P(x)
# Calculate the integral: \int_{-\infty}^{\infty} P{y | x} * P(x) dx
integral = integrate(integrand(x, y), x)
return integral
end
Many thanks for taking the time to read through,
Pazzy Boardman