Solving and plotting equation with many variables

2024 2 26
1 个回答
更新时间：2024 2 28
11 次查看（30 天）
Ran in:

0 个投票

Ran in:
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')

13 个评论
显示 11更早的评论 隐藏 11更早的评论

Sabrina Garland
Sabrina Garland 2024-2-27
编辑：Voss 2024-2-27
@Torsten Can you please look at my code and tell me where I went wrong?
Voss
Voss 2024-2-27
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
Sabrina Garland 2024-2-27
编辑：Voss 2024-2-27
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
Voss 2024-2-27
If this is a script, you've got to clear the old integral2 variable before running the script again.
clear integral2
Sabrina Garland
Sabrina Garland 2024-2-27
I have cleared and ran it again
Torsten
Torsten 2024-2-27
编辑：Torsten 2024-2-27
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
Sabrina Garland 2024-2-28
@Torsten Rewrote it completely but still getting error
Sabrina Garland
Sabrina Garland 2024-2-28
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
Sabrina Garland 2024-2-28
@Torsten Please help. Please help anyone
Sabrina Garland
Sabrina Garland 2024-2-28
@Torsten How to solve the problem. Why is it giving this error? Please please guide me a way out
Torsten
Torsten 2024-2-28
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2024-2-28
编辑：Walter Roberson 2024-2-28
Ran in:

0 个投票

Ran in:
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')

8 个评论
显示 6更早的评论 隐藏 6更早的评论

Sabrina Garland
Sabrina Garland 2024-2-28
编辑：Voss 2024-2-28
@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
Sam Chak 2024-2-28
@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
Sabrina Garland 2024-2-28
编辑：Voss 2024-2-28
@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
Sam Chak 2024-2-28
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
Sabrina Garland 2024-2-28
编辑：Voss 2024-2-28
@Sam Chak Renamed the equation involving the word integral. Can you please review it. Please
Sam Chak
Sam Chak 2024-2-28
Ran in:
Ran in:
@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
Walter Roberson 2024-2-28
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Lognormal Distribution 的更多信息

产品

版本

R2023b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by