Solving and plotting equation with many variables
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Torsten
el 28 de Feb. de 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.
Respuestas (1)
Walter Roberson
el 28 de Feb. de 2024
Editada: Walter Roberson
el 28 de Feb. de 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 comentarios
Sam Chak
el 28 de Feb. de 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
el 28 de Feb. de 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.
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!