Solving and plotting equation with many variables

I have two random variables a and b. a is drawn from uniform distribution (m, 1+m) where, 0<m<=1. b is drawn from uniform distribution (0,1). Conditional expected value is calculated as - These expected values are used to input into following equation to solve for x as a function of k and m- Solution of x is then substituted into following integration which is optimized for k (after solving the integration, it is differentiated with respect to k and an expression of k is derived) Since a and b has been removed through integration, solution of k will just be a function of m. I want Matlab to pick only those solution of k which is strictly between (1/2,1) [there is only unique such k] I want to plot this k against m (I'll get unique k for each m) I also want to solve second integration - I'll also get another k on optimizing it (again k picked should be between 1/2,1). I want to plot these both k against m in same plot so that I can get comparative graph.

% Define the range of m values 
m = linspace(0.01, 1, 20);
% Initialize arrays to store the values of k 
k1 = zeros(size(m)); 
k2 = zeros(size(m)); 
k3 = zeros(size(m));
% Calculate the values of k for each m 
for i = 1:length(m) 
    E_ab = m(i) * (integral2(@(a, b) a.*unifpdf(a).*unifpdf(b), 0, m(i), m(i), 1+m(i)) /...
                   integral2(@(a, b) 1.*unifpdf(a).*unifpdf(b), 0, m(i), m(i), 1+m(i))) + ...
                   (1 - m(i)) * (integral2(@(a, b) a.*unifpdf(a).*unifpdf(b), m(i), 1, @(b)b, 1+m(i)) /...
                   integral2(@(a, b) 1.*unifpdf(a).*unifpdf(b), m(i), 1, @(b)b, 1+m(i))); 
    E_ba = (1 - m(i)) * (integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), m(i), 1, 0, @(a)a) /...
           integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), m(i), 1, 0, @(a)a)) + ...
           m(i) * (integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), 1, 1+m(i), 0, 1) /...
           integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), 1, 1+m(i), 0, 1)); 
    E_aab = integral2(@(a, b) a.*unifpdf(a).*unifpdf(b), 0,1, m(i), @(b)b) /...
            integral2(@(a, b) 1.*unifpdf(a).*unifpdf(b), 0, 1, m(i), @(b)b); 
    E_bba = integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), m(i), 1, @(a)a, 1) / ...
            integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), m(i), 1, @(a)a, 1); 
    k = sym('k');
    % Solve for x
    syms x
    eqn     = sqrt(k) * E_ab * ( - x) + sqrt(1-k) * E_ba  == sqrt(k) * E_bba  + sqrt(1-k) * E_aab * ( - x);
    sol     = solve(eqn, x);
    a       = unifrnd(0,1);
    b       = unifrnd(0,1);
    r1      = matlabFunction((sqrt(k) * a .* a + sqrt(1-k) * b .* b)); 
    r2      = matlabFunction((sqrt(k) * b .* b + sqrt(1-k) * a .* a));
    r3      = matlabFunction((sqrt(k) * E_ab .* a + sqrt(1-k) * E_ba .* b));
    r4      = matlabFunction((sqrt(k) * E_bba .* b + sqrt(1-k) * E_aab .* a));

    intgeqn4 = integral2(r1, 0, m(i), m(i), 1+m(i));
    intgeqn5 = integral2(r1, m(i), 1, @(b)b, 1+m(i));
    intgeqn6 = integral2(r2, m(i), 1, m(i), @(b)b);

    % Substitute x into the integration equation
    intgeqn7 = integral2(r3, 0, m(i), m(i), 1+m(i));
    intgeqn8 = integral2(r3, m(i), 1, @(b)b-x, 1+m(i));
    intgeqn9 = integral2(r4, m(i), 1, m(i), @(b)b-x);

    % Solve for k
    syms k
    eqn2    = intgeqn7 + intgeqn8 + intgeqn9;
    eqn3    = intgeqn4 + intgeqn5 + intgeqn6;
    sol2    = solve(eqn2, k);
    sol3    = solve(eqn3, k);

    % Store the values of k
    k1(i)   = double(sol2(1));
    k2(i)   = double(sol2(2));
    k3(i)   = double(sol3(1));
end

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

Error using symengine>@(k)sqrt(-k+1.0).*3.146244985517477e-2+sqrt(k).*1.334869879283682e-1
Too many input arguments.

Error in integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

% Plot the values of k against m 
plot(m, k1, 'b', 'LineWidth', 2) 
hold on 
plot(m, k2, 'r', 'LineWidth', 2) 
hold on 
plot(m, k3, 'g', 'LineWidth', 2) 
xlabel('m') 
ylabel('k') 
legend('k1', 'k2', 'k3') 
title('Values of k against m')

[Sabrina Garland]

@Torsten Please help

[Sabrina Garland]

@Torsten Can you please look at my code and tell me where I went wrong?

[Voss]

Here's one thing that's wrong:

integral2 = integral2(E_ba, m(i), 1+m(i), m(i), 1+m(i));

You create a variable with the same name as a function you want to subsequently call. Subsequent calls to integral2 will reference the variable, not the function, and you'll get an error about attempting to index a variable using a function handle.

[Sabrina Garland]

Thanks @Voss. But even changing that, I am getting the error

Unable to use a value of type function_handle as an index.

What am I doing wrong?

[Voss]

If this is a script, you've got to clear the old integral2 variable before running the script again.

clear integral2

[Sabrina Garland]

I have cleared and ran it again

[Voss]

OK

[Torsten]

The name "integral3" again conflicts with the MATLAB built-in function "integral3". And where do you define "integral3" and "integral4" in your code ?

[Sabrina Garland]

@Torsten Rewrote it completely but still getting error

[Sabrina Garland]

Getting error - Error using integral2Calc>integral2t/tensor Integrand output size does not match the input size.

Error in integral2Calc>integral2t (line 55) [Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9) [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105) Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

[Sabrina Garland]

@Torsten Please help. Please help anyone

[Sabrina Garland]

@Torsten How to solve the problem. Why is it giving this error? Please please guide me a way out

[Torsten]

The integrals

integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), 1, 1+m(i), 0, 1)
integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), 1, 1+m(i), 0, 1)

in the computation of E_ba will always be 0 because the mass of unifpdf is concentrated on [0 1] and you evaluate the function from 1 to 1+m in b-direction. Thus you get a 0/0 division and the result for E_ba is NaN.

Answer 1

[Walter Roberson]

The solution of solve(eqn, x) is in four parts, with different conditions for each part -- conditions depending on the values of a and b relative to various constants. But it doesn't matter, because you never use the solution after you compute it.

The solution for solve(eqn2, k) is empty, at least over reals. You are adding three values that cannot be negative, and you are solving for 0. The only potential solution is over complex values.

I have to question the inequalities for E_ab and so on, whether they should be a>=b and similar.

% Define the range of m values
 m = linspace(0.01, 1, 30);
 % Initialize arrays to store the 
 %values of k
 k1 = zeros(size(m));
 k2 = zeros(size(m));
 % Calculate the values of k 
% for each m
 syms a b k x
 for i = 1:length(m)
     % Define the integrands for the conditional expected values
     %E_ab = @(a, b) (a > b) .* a;
     %E_ba = @(a, b) (a < b) .* b;
     %E_aab = @(a, b) (b > m(i) & b < 1+m(i)) .* a;
     %E_bba = @(a, b) (a > m(i) & a < 1) .* b;
     E_ab(a,b) = piecewise(a > b, a, 0);
     E_ba(a,b) = piecewise(a < b, b, 0);
     E_aab(a,b) = piecewise(b > m(i) & b < 1+m(i), a, 0);
     E_bba(a,b) = piecewise(a > m(i) & a < 1, b, 0);
     
     % Calculate the integrals for E(a|a>b) and E(b|a>b)
     integral1_ = int(int(E_ab, a, m(i), 1+m(i)), b, m(i), 1+m(i));
     integral2_ = int(int(E_ba, a, m(i), 1+m(i)), b, m(i), 1+m(i));
     % Calculate the integrals for E(a|ab) and E(b|a>b)
     integral5_ = int(int(E_aab, a, m(i), 1), b, 0, 1);
     integral6_ = int(int(E_bba, a, m(i), 1), b, 0, 1);
     % Calculate the conditional expected values
     E_a_ab = (integral1_ + integral5_) / (integral1_ + integral2_ + integral5_ + integral6_);
     E_b_ab = (integral2_ + integral6_) / (integral1_ + integral2_ + integral5_ + integral6_);
     % Define the missing variables

     integral3_ = int(int(sqrt(k) * E_ab(a, b) .* a + sqrt(1-k) * E_ba(a, b) .* b, a, m(i), 1+m(i)), b, m(i), 1+m(i));
     integral4_ = int(int(E_ab(a, b) .* a + sqrt(1-k) * E_ba(a, b) .* b, a, m(i), 1+m(i)), b, m(i), 1+m(i));
     % Solve for x

     %the solution is not used anywhere...
     eqn = sqrt(k) * E_a_ab * (a - x) + sqrt(1-k) * E_b_ab * a == sqrt(k) * E_bba * a + sqrt(1-k) * E_aab * (a - x);
     sola = solve(eqn, x, 'returnconditions', true);
     sol = sola.x;
     
     % Substitute x into the integration equation
     integral7_ = int(int(sqrt(k) * E_ab(a, b) .* a + sqrt(1-k) * E_ba(a, b) .* b, a, m(i), 1+m(i)), b, m(i), 1+m(i));
     integral8_ = int(int(E_ab(a, b) .* a + sqrt(1-k) * E_ba(a, b) .* b, a, m(i), 1+m(i)), b, m(i), 1+m(i));     
     integral9_ = int(int(sqrt(k) * E_ba(a, b) .* b + sqrt(1-k) * E_ab(a, b) .* a, a, 0, 1), b, m(i), 1);
     % Solve for k
     eqn2 = integral7_ + integral8_ + integral9_;
     sol2 = solve(eqn2, k);
     if isempty(sol2)
         sol2 = nan;
     end
     if length(sol2) < 2
         sol2(2) = nan;
     end
     % Store the values of k
     k1(i) = double(sol2(1));
     k2(i) = double(sol2(2));
end
 % Plot the values of k against m
 plot(m, k1, 'b', 'LineWidth', 2)
 hold on
 plot(m, k2, '--r', 'LineWidth', 2)
 xlabel('m')
 ylabel('k')
 legend('k1', 'k2')
 title('Values of k against m')

[Sabrina Garland]

@Walter Roberson Thank you for taking time to write and for looking at my problem. The code in my comments was rewritten because the x was not being used. I have edited the code in main question an hour ago. Please can you look at it. When I am running it I am getting an error that integrand input size doesn't match output size. I don't know what that means. Can you please look at my code in the main question. Please

[Sam Chak]

@Sabrina Garland: The code in my comments was rewritten because the x was not being used. I have edited the code in main question an hour ago.

Please, for the sake of clarity, update the code once again, assigning distinct names to the numbered integral equations to avoid conflicts with the built-in integral function in MATLAB. Additionally, kindly log the update and present the integral equations in LaTeX format (preferred) or in image form (if necessary).

Update 1:

inteqn1 = ...
inteqn2 = ...

[Sabrina Garland]

@Sam Chak Thank you for taking time to review my query. I already have updated my code and attached all the relevant equations too. Please look at the main problem. Can you please help me with the code, I am getting error of too many input arguments for the equation 'integral4'. Can you please please help

[Sam Chak]

Hi @Sabrina Garland, if you have fully updated the code based on my suggestions, then the numbered integral equation 'integral4' should no longer exist. Please review my previous comment once again to ensure that you haven't missed anything.

[Sabrina Garland]

@Sam Chak Renamed the equation involving the word integral. Can you please review it. Please

[Sabrina Garland]

@Torsten Please help. @Sam Chak @Walter Roberson please help

[Sam Chak]

@Sabrina Garland, I randomly selected a value between 0.01 and 1 for m. As you can observe, E_ba produces NaN, and E_aab yields Inf.

Additionally, I noticed that the previously displayed images have been removed for unknown reasons. I haven't examined your code (there may be errors) as it is rather inconvenient to click on the images one by one.

% Define the range of m values 
m = 0.5;
% Initialize arrays to store the values of k 
k1 = zeros(size(m)); 
k2 = zeros(size(m)); 
k3 = zeros(size(m));
% Calculate the values of k for each m 
for i = 1:length(m) 
    E_ab = m(i) * (integral2(@(a, b) a.*unifpdf(a).*unifpdf(b), 0, m(i), m(i), 1+m(i)) /...
                   integral2(@(a, b) 1.*unifpdf(a).*unifpdf(b), 0, m(i), m(i), 1+m(i))) + ...
                   (1 - m(i)) * (integral2(@(a, b) a.*unifpdf(a).*unifpdf(b), m(i), 1, @(b)b, 1+m(i)) /...
                   integral2(@(a, b) 1.*unifpdf(a).*unifpdf(b), m(i), 1, @(b)b, 1+m(i)))
    E_ba = (1 - m(i)) * (integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), m(i), 1, 0, @(a)a) /...
           integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), m(i), 1, 0, @(a)a)) + ...
           m(i) * (integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), 1, 1+m(i), 0, 1) /...
           integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), 1, 1+m(i), 0, 1))
    E_aab = integral2(@(a, b) a.*unifpdf(a).*unifpdf(b), 0,1, m(i), @(b)b) /...
            integral2(@(a, b) 1.*unifpdf(a).*unifpdf(b), 0, 1, m(i), @(b)b)
    E_bba = integral2(@(b, a) b.*unifpdf(a).*unifpdf(b), m(i), 1, @(a)a, 1) / ...
            integral2(@(b, a) 1.*unifpdf(a).*unifpdf(b), m(i), 1, @(a)a, 1)
    k = sym('k');
    % Solve for x
    syms x
    eqn     = sqrt(k) * E_ab * ( - x) + sqrt(1-k) * E_ba  == sqrt(k) * E_bba  + sqrt(1-k) * E_aab * ( - x);
    sol     = solve(eqn, x);
    a       = unifrnd(0,1);
    b       = unifrnd(0,1);
    r1      = matlabFunction((sqrt(k) * a .* a + sqrt(1-k) * b .* b)); 
    r2      = matlabFunction((sqrt(k) * b .* b + sqrt(1-k) * a .* a));
    r3      = matlabFunction((sqrt(k) * E_ab .* a + sqrt(1-k) * E_ba .* b));
    r4      = matlabFunction((sqrt(k) * E_bba .* b + sqrt(1-k) * E_aab .* a));

    intgeqn4 = integral2(r1, 0, m(i), m(i), 1+m(i));
    intgeqn5 = integral2(r1, m(i), 1, @(b)b, 1+m(i));
    intgeqn6 = integral2(r2, m(i), 1, m(i), @(b)b);

    % Substitute x into the integration equation
    intgeqn7 = integral2(r3, 0, m(i), m(i), 1+m(i));
    intgeqn8 = integral2(r3, m(i), 1, @(b)b-x, 1+m(i));
    intgeqn9 = integral2(r4, m(i), 1, m(i), @(b)b-x);

    % Solve for k
    syms k
    eqn2    = intgeqn7 + intgeqn8 + intgeqn9;
    eqn3    = intgeqn4 + intgeqn5 + intgeqn6;
    sol2    = solve(eqn2, k);
    sol3    = solve(eqn3, k);

    % Store the values of k
    k1(i)   = double(sol2(1));
    k2(i)   = double(sol2(2));
    k3(i)   = double(sol3(1));
end

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

E_ab = 0.4583

E_ba = NaN

E_aab = Inf

E_bba = 0.6667

Error using symengine>@(k)sqrt(-k+1.0).*3.954039498510441e-2+sqrt(k).*2.101315801518972e-1
Too many input arguments.

Error in integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

[Walter Roberson]

    syms x
    eqn     = sqrt(k) * E_ab * ( - x) + sqrt(1-k) * E_ba  == sqrt(k) * E_bba  + sqrt(1-k) * E_aab * ( - x);
    sol     = solve(eqn, x);

You do not use the output sol anywhere, so there is not point in computing it.

    r1      = matlabFunction((sqrt(k) * a .* a + sqrt(1-k) * b .* b)); 

a and b are numeric. k is syms. So r1 is a function of one variable, k.

    intgeqn4 = integral2(r1, 0, m(i), m(i), 1+m(i));

You are asking to integrate r1 (a function of k) over two (unnamed) variables.