How do I apply iteration method to an integral?

Hello everyone,

I need to code an equation. That equation includes a new parameter that is not in my data, and therefore the parameter needs to be estimated/calculated before I can proceed.

There is a formula for the parameter, and the formula includes an integral, with the lowest limit being exactly that parameter I am trying to estimate. Therefore, that parameter is both on the left and on the right side of the equation. I am trying to figure out how to do it.

I think I understand the general idea. I need to use iterative method: make an initial guess, write a loop with the formula and it will converge to some value, which will be the final answer (please let me know if I am wrong)

I can’t figure out the steps I need to do. I’ve tried multiple things, but I couldn’t get it to run.

Here’s the formula:

• cdf is stored in stp.F;
• alpha is stored in mod.alpha;
• vector with q values : mod.vecq.

Could you please help me or give me some pointers? Thank you in advance.

Answer 1

[Alan Stevens]

Look up help on function fzero.

[Angie]

Thank you very much! It looks promising and shoul do the trick. I found many examples online, but none of them have integrals, so I'm currently trying to figure out how to apply it to my case. I will give an update when I've succeeded (or not:)

[Angie]

Hi Alan. Thank you very much for your pointer. It did indeed help me. I still had to take quite a while to figure out how to code the integral properly and then finish the rest of my task (it was just a small but very crucial part of it), but (I think) I've succeeded.

For those wondering, here's the code (could have slight errors since my actual integral was more sophisticated than the simplified version in this topic):

file ql_eq.m:

function y = ql_eq(ql, mod, stp, Fx)
qmax = max(mod.vecq);
% Interpolation of the CDF F(x) over the range of pvec (if needed)
    Fx = @(x) interp1(mod.vecq, stp.F, x, 'linear', 'extrap');  % Interpolate CDF values
fun = @(x) (1 - Fx(x)) ./ (mod.alpha + (1 - Fx(x)));
integral_val = integral(fun, ql, qmax, 'ArrayValued', true);
y = integral_val - pinf;
end

file ql_sol.m:

ql_guess = 1; % set your guess
options = optimset('Display', 'iter');  % Show iterations
ql_solution = fzero(@(x) ql_eq(x, mod, stp), ql_guess, options);
disp(ql_solution);
ql=find(mod.vecq >= ql_solution, 1, 'first');

[Alan Stevens]

Good!

(We can't actually run your code as we don't have access to mod, stp and pinf.)

[Angie]

Yes of course, I understand. But if at some point in the future someone has their own values and is looking for a solution to a similar problem hopefully the code can help them to move in the right direction. Thank you again, your pointer saved me a lot of headache :)