Solving and plotting equation with many variables
3 views (last 30 days)
Show older comments
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
% 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 Comments
Torsten
on 28 Feb 2024
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.
Answers (1)
Walter Roberson
on 28 Feb 2024
Edited: Walter Roberson
on 28 Feb 2024
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 Comments
Sam Chak
on 28 Feb 2024
@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
Walter Roberson
on 28 Feb 2024
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.
See Also
Categories
Find more on Quadratic Programming and Cone Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!