How do I add summation to a symbolic equation in order to integrate across set bounds?

The goal is to integrate from a set number to infinity (or some larger number) for function GF-PDF, which is a lognormal distribution with a geometric mean GF (GFg) and spread (sigma). Each GF-PDF is specifically collected for/interpolated to a dry size (Dp) and the other distribution parameters (GFg, sigma) are fit to the data through a data inversion algorithm. For one mode in the distribution, the equation is:

The issue is that sometime there are multiple modes within one GF-PDF and so the equation becomes:

...where for each mode (k), there is a number fraction (f0,k), geometric mean GF (GFg,k), and spread (sigmak) that described each mode and the sum of all the modes is the distribution. The f0, GFg, and sigma are different for all the modes making up the GF-PDF.

The integration I'm trying to do is:

...where I can integrate above a specific value for GFc (that was calculated for a chosen Sc at a chosen Dp) and get the area/fraction of particles above that GFc since the GF-PDF is at unity. I know how to use trapz(), but integration as a whole has been confusing me. From what I have been reading, I need to input the GF-PDF as a symbolic function, but I'm uncertain on how to add the summation to the symbolic equation in the cases where 2 or more modes are present. I could create the single mode equation as a symbolic function, but I'm again not sure how to integrate across the summation of the equations.

Below are functions I created to create the GF-PDFs, but I know that they aren't symbolic. I add some sample parameters and the range of GF too.

binsp = (exp(1/60)*0.7)-0.7; % logspace for GF range
gf_rng = 0.90:binsp:2.5; % range of GF in log-space for the GF-PDFs

% Example values for a bimodal (k = 2) GF-PDF:
g0s = [1.1,1.5]; % geometric mean GF
sig0s = [1.04,1.05]; % spread in GF
f0s = [0.4,0.6]; % number fraction of mode

% GF-PDF for GF range
gfpdf = GFPDFf(gf_rng,g0s,sig0s,f0s); % uses functions at bottom

% I could make a symbolic function for just the one lognormal mode, I'm not
% sure how to sum them both in it though...
gf_rng = sym('gf_rng');
g0s = sym('g0s');
sig0s = sym('sig0s');
f0s = sym('f0s');
GFPDF_1m = (f0s./sqrt(2.*pi.*(log(sig0s)).^2)).*exp((-0.5).*((log(gf_rng./g0s).^2)./((log(sig0s)).^2)))

GFPDF_1m = 

% I'm uncertain where to go from here if I want to integrate from GFc to
% infinity for the bimodal curve...
GFc = 1.3; % random GFc value


% Functions for GF-PDF creations
% Calculating the GFPDF: 
function npdf = GFPDFf(gf_rng,g0s,sig0s,f0s)
if width(g0s)>1
    for i = 1:width(g0s)
        pdf(i,:) = lnPDF(gf_rng,g0s(1,i),sig0s(1,i),f0s(1,i));
    end
    npdf = sum(pdf);
else
    npdf = lnPDF(gf_rng,g0s,sig0s,f0s);
end
end

% Plotting Lognormal Modes:
function pdf = lnPDF(gf_rng,g0s,sig0s,f0s)
pdf = (f0s./sqrt(2.*pi.*(log(sig0s)).^2)).*exp((-0.5).*((log(gf_rng./g0s).^2)./((log(sig0s)).^2)));
end

Answer 1

[Torsten]

Use logncdf:

https://uk.mathworks.com/help/stats/logncdf.html

Use

to compute f_act with the help of "logncdf" and use the normal "sum" command to sum over k.

No symbolic computations are needed.

[Paul]

Because a) GF-PDF does not have a 1/GF term out front, and b) the variable of integration in f_act is ln(GF), wouldn't substitution of u = ln(GF) indicate that normcdf is the appropriate function? Assuming that the equations as written are correct.

[Torsten]

d(ln GF) = 1/GF dGF

[Paul]

True. If using logncdf does the lower limit of integration need to change from GF_c to exp(GF_c)?

Maybe I'm just confused about the notation. The integrand in f_act is a function of GF, but the differential is d(log(GF)). Are the limits of integration in GF or log(GF) ?

[Torsten]

I think it's written as a Riemann-Stieltjes integral. For differentiable h, it holds that

https://de.wikipedia.org/wiki/Stieltjesintegral

@Mara

For the case given, I think the code would be:

% Example values for a bimodal (k = 2) GF-PDF:
g0s = [1.1,1.5]; % geometric mean GF
sig0s = [1.04,1.05]; % spread in GF
f0s = [0.4,0.6]; % number fraction of mode

% I'm uncertain where to go from here if I want to integrate from GFc to
% infinity for the bimodal curve...
GFc = 1.3; % random GFc value

s = f0s(1)*(1-logncdf(GFc,log(g0s(1)),log(sig0s(1)))) + f0s(2)*(1-logncdf(GFc,log(g0s(2)),log(sig0s(2))))

s = 0.5990

[Paul]

logncdf is vectorized and has the 'upper' flag, which might improve accuracy in some cases, so we can save some keystrokes

g0s = [1.1,1.5]; % geometric mean GF
sig0s = [1.04,1.05]; % spread in GF
f0s = [0.4,0.6]; % number fraction of mode
% I'm uncertain where to go from here if I want to integrate from GFc to
% infinity for the bimodal curve...
GFc = 1.3; % random GFc value
s = f0s(1)*(1-logncdf(GFc,log(g0s(1)),log(sig0s(1)))) + f0s(2)*(1-logncdf(GFc,log(g0s(2)),log(sig0s(2))))

s = 0.5990

s = sum(f0s.*logncdf(GFc,log(g0s),log(sig0s),'upper'))

s = 0.5990

[Mara]

Thank you for all of your answers, they have all be very helpful. As you pointed out there is an error I think in my differential equation. The method/paper I'm adapting this integration from had used an inversion program that creates a piece-wise linear function for the GF-PDF while the inversion program I used creates multi-modal lognormal distributions. And so I thought that all the integrals needed to be changed because of the type of function being in linear versus log space. When reading about the lognormal distribution, I saw that they used dln(x) in the integral when talking about unity (equation below) and so I figured that if I were going to integrate with GF, that I needed to be integrating by dln(GF).

Whenever I plotted it, I would use just the normal GFs and not the log form of them, but when I used trapz to get area, I thought because the GF is technically ln(GF) in the function, I need to use ln(GF) when integrating. I could be completely wrong about this though.

[Torsten]

I must admit that I don't understand what you want to tell us.

How does your question actually change after your comment ?

Why do you define

binsp = (exp(1/60)*0.7)-0.7; % logspace for GF range
gf_rng = 0.90:binsp:2.5; % range of GF in log-space for the GF-PDFs

?

Do you want to truncate GF-PDF and restrict it to the interval [0.9,2.5] ? Then use

format long
g0s = [1.1,1.5]; % geometric mean GF
sig0s = [1.04,1.05]; % spread in GF
f0s = [0.4,0.6]; % number fraction of mode
GFc = 1.3; % random GFc value
pdf = @(x)f0s(1)*lognpdf(x,log(g0s(1)),log(sig0s(1)))+f0s(2)*lognpdf(x,log(g0s(2)),log(sig0s(2)));
normalizer = integral(pdf,0.9,2.5);
pdfn = @(x)(0.9<=x).*(2.5>=x).*pdf(x)/normalizer;
s = integral(pdfn,GFc,2.5)

s = 
   0.598996979312686

[Mara]

Sorry, so the question doesn't change, I more wanted to explain why I wrote out the equation that way since it seems that there was discussion of it in the replies. The gf_rng is the range of GF I use, which in your answer would be the x, and defined it like above to make sure it was evenly spaced in log-space, though it's actually incorrectly done now that I'm looking at it again. Either way, thank you for your time and help on answering my question.

[Torsten]

If you know that the GF are lognormally distributed with mu = g0s and sigma = sig0s, it doesn't matter which range of the distribution you use. So my first answer (without the normalization) should be correct.

Or are the g0s and sig0s determined for a lognormal distribution restricted to the interval [0.9:2.5] ? Then I had to rethink my answer again.

[Mara]

The GF are lognormally distribution with g0s=mu and sigma = sig0s. The range is for plotting the pdf and are technically restricted to that range just due to measured values from all experimental data (the modes realistically won't appear above 2.5 or below 0.9). I believe your first answer (the non-restricted one) works as the mu and sigma won't be values that will place the mode outside this range if the inversion works correctly and the data is good quality.

Answer 2

[Umar]

Hi @Mara,

I think you're actually overthinking this a bit. The numerical approach you've already started is the right way to go. Using symsum and int for this problem won't really help you because the sum of multiple lognormal distributions doesn't have a closed-form solution, so MATLAB would just return an unevaluated symbolic expression anyway.

Since you already have actual numerical values for your modes, the numerical integration approach is both more practical and more accurate. Here's how to complete what you started.

The key thing you're missing is normalization. Your total integral should be 1.0, but without normalizing after you sum the modes, it won't be. Here's the complete working code:

clear all;
clc;
% Parameters
GFc = 1.3;
binsp = (exp(1/60)*0.7)-0.7;
gf_rng = 0.90:binsp:2.5;
% Bimodal parameters
g0s = [1.1, 1.5];
sig0s = [1.04, 1.05];
f0s = [0.4, 0.6];
% Calculate PDF - the summation happens here
gfpdf = zeros(size(gf_rng));
for i = 1:length(g0s)
    % Calculate lognormal PDF for mode i
    pdf_mode = (f0s(i)./sqrt(2.*pi.*(log(sig0s(i))).^2)) .* ...
        exp((-0.5).* ...
    ((log(gf_rng./g0s(i)).^2)./((log(sig0s(i))).^2)));
    
    % Add to total PDF (this is your summation over k)
    gfpdf = gfpdf + pdf_mode;
end
% Normalize the PDF - this is critical
gfpdf_normalized = gfpdf / trapz(gf_rng, gfpdf);
% Integration from GFc to end of range
idx = gf_rng >= GFc;
f_act = trapz(gf_rng(idx), gfpdf_normalized(idx))

f_act = 0.6700

fprintf('Fraction above GFc = %.4f\n', f_act);

Fraction above GFc = 0.6700

% Verification - should be 1.0
total_integral = trapz(gf_rng, gfpdf_normalized);
fprintf('Total integral (should be ~1): %.4f\n', total_integral);

Total integral (should be ~1): 1.0000

% Plot
figure;
plot(gf_rng, gfpdf_normalized, 'b-', 'LineWidth', 2);
hold on;
xline(GFc, 'r--', 'LineWidth', 2);
fill([gf_rng(idx), fliplr(gf_rng(idx))], ...
    [gfpdf_normalized(idx), zeros(size(gfpdf_normalized(idx)))], ...
    'r', 'FaceAlpha', 0.3);
xlabel('GF');
ylabel('PDF');
title('Bimodal GF-PDF (Normalized)');
legend('GF-PDF', 'GFc cutoff', sprintf('Area = %.4f', f_act));
grid on;

*Results*: please see attached.

To answer your specific question about where to go from here: the summation is already happening in your loop where you do gfpdf = gfpdf + pdf_mode. That's literally the sum over k modes. You don't need symsum for that. The integration from GFc to infinity is just using trapz on the portion of your range above GFc, which the code above shows with the idx indexing.

For your bimodal case with GFc = 1.3, you should get a fraction around 0.67, which makes sense because most of your second mode centered at 1.5 is above the cutoff, plus a bit of the tail from the first mode at 1.1.

This approach scales easily too. If you end up with 3 or 4 modes later, just add more values to your g0s, sig0s, and f0s arrays and the loop handles it automatically. The functions you already created work fine, you just needed the normalization step and to use trapz for the integration part.

I know you mentioned integration as a whole has been confusing, but I think that confusion was coming from trying to make it symbolic. Once you realize you can just use trapz directly on your numerical arrays (no symbolic function needed), it becomes straightforward. The idx indexing just selects the part of your gf_rng array that's above GFc, and trapz integrates that portion.

Hope this clears things up.

[Umar]

Sorry, i tried to format it but it end up like this.

[Mara]

Thank you for your answer! I do have a quick question about the normalization. I usually normalize using trapz, but since it's lognormal modes created from a lognormal function, I've used ln(GF_rng) as my x values instead of just GF_rng. Is there a right and wrong way to calculate the area using trapz for lognormal functions? In previous classes, we would take the area under the curve for particle size distributions that were lognormal using trapz(log(D),dN/dlogD) and so I didn't want to use the wrong x values for the function, if that makes sense.

[Umar]

First, I would like to say thank you @Torsten for formatting my posted comments,really appreciate and helpful. Okay @Mara, this is actually a great question about the normalization! This is actually a really important distinction that comes up with lognormal distributions. The answer is that the current approach using trapz(gf_rng, gfpdf) with linear x-values is correct for your case. Here's why. Your lognormal PDF formula is giving you probability density per unit GF, not per unit log(GF). When you calculate the PDF at different GF values, those numbers tell you the probability density at those actual GF points. So when you integrate to find areas or normalize, you need to integrate over the actual GF values, not their logarithms. In your previous class with particle size distributions, you used trapz(log(D), dN/dlogD) because that data was specifically formatted as "per logarithmic interval" - the dlogD literally means per log interval. That's a different format than what you have now. Your current PDF isn't normalized per log interval, it's normalized per regular interval. Here's a simple way to think about it: if you plot your PDF on a graph with regular GF values on the x-axis (not log-scaled), and you want to find the area under that curve, you'd measure the width of each trapezoid using the actual GF spacing. That's what trapz(gf_rng, gfpdf) does. If you were to use trapz(log(gf_rng), gfpdf), you'd be mixing things up - you'd be treating your PDF as if it were per log interval when it's actually per regular interval, and your normalization would come out wrong. You can verify this by checking that trapz(gf_rng, gfpdf_normalized) equals 1.0, which it should for a properly normalized PDF. That confirms you're integrating correctly. So stick with what the answer provided - trapz(gf_rng, gfpdf) is the right way to go for your lognormal PDF. Hope this helps!