I have a joint probability density f(x,y,z) and I wish to find the conditional distribution X|Y=y,Z=z, which is equivalent to treating x as data and y and z as parameters (constants).
For example, if I have X|Y=y,Z=z being the pdf of a N(1-2y,3z^2+2), the function would be:
syms x y z
f(y,z) = 1/sqrt(2*pi*(3*z^2+2)) * exp(-1/(2*(3*z^2+2)) * (x-(1-2*y))^2);
I would like to compare it to the following:
syms mu s L a b
Normal(mu,s) = (1/sqrt(2*pi*s^2)) * exp(-1/(2*s^2) * (x-mu)^2);
Exponential(L) = L * exp(-L*x);
Gamma(a,b) = (b^a / gamma(a)) * x^(a-1)*exp(-b*x);
Beta(a,b) = (1/beta(a,b)) * x^(a-1)*(1-x)^(b-1);
Question
How do I make a program whichDistribution that would be able to print which of these four, f is equivalent to (up to proportionality) with respect to the variable x, and what are the parameters? E.g. f and x as above, the distribution is Normal, mu=1-2*y, s=3*z^2+2.
NB: there would not always be a unique solution, since some distributions are are equivalent (e.g. Gamma(1,L)==Exponential(L))
Desired outputs
syms x y z
f = 1/sqrt(2*pi*(3*z^2+2)) * exp(-1/(2*(3*z^2+2)) * (x-(1-2*y))^2)
whichDistribution(f,x) %Conditional X|Y,Z
% Normal(1-2*y,3*z^2+2)
syms x y
f = y^(1/2)*exp(-(x^2)/2 - y/2 * (1+(4-x)^2+(6-x)^2)) % this is not a pdf because it is missing a constant of proportionality, but it should still work
whichDistribution(f,x) %Conditional X|Y
% Normal(10*y/(2*y+1), 1/(2*y+1))
whichDistribution(f,y) %Conditional Y|X
% Gamma(3/2, x^2 - 10*x + 53/2)
f = exp(-x) %also missing a constant of proportionality
whichDistribution(f,x)
% Exponential(1)
f = 1/(2*pi)*exp(-(x^2)/2 - (y^2)/2)
whichDistribution(f,x)
% Normal(0,1)
whichDistribution(f,y)
% Normal(0,1)
What I have tried so far:
- Using
solve():
q = solve(f(y,z) == Normal(mu,s), mu, s)
Which gives wrong results, since parameters can't depend on x:
>> q.mu
ans =
(z1^2*(log((2^(1/2)*exp(x^2/(2*z1^2) - (x + 2*y - 1)^2/(6*z^2 + 4)))/(2*pi^(1/2)*(3*z^2 + 2)^(1/2))) + pi*k*2i))/x
>> q.s
ans =
z1
- Attempting to simplify
f(y,z)up to proportionality (inxvariable) using apropto()function that I wrote:
>> propto(f(y,z),x)
ans =
exp(-(x*(x + 4*y - 2))/(2*(3*z^2 + 2)))
>> propto(Normal(mu,s),x)
ans =
exp((x*(2*mu - x))/(2*s^2))
This is almost on the money, since it is easy to spot that s^2=3*z^2 + 2 and 2*mu=-(4*y - 2), but I don't know how to deduce this programmatically.
In case it is useful: propto(f,x) attempts to simplify f by dividing f by children of f which don't involve x, and then output whichever form has the least number of children. Here is the routine:
function out = propto(f,x)
oldf = f;
newf = propto2(f,x);
while (~strcmp(char(oldf),char(newf))) % if the form of f changed, do propto2 again. When propto2(f) == f, stop
oldf = newf;
newf = propto2(oldf,x);
end
out = newf;
end
function out = propto2(f,x)
t1 = children(expand(f)); % expanded f
i1 = ~has([t1{:}],x);
out1 = simplify(f/prod([t1{i1}])); % divides expanded f by terms that do not involve x
t2 = children(f); % unexpanded f
i2 = ~has([t2{:}],x);
out2 = simplify(f/prod([t2{i2}])); % divides f by terms that do not involve x
A = [f, symlength(f); out1, symlength(out1); out2, symlength(out2)];
A = sortrows(A,2); % outputs whichever form has the fewest number of children
out = A(1,1);
end
function L = symlength(f)
% counts the number of children of f by repeatingly applying children() to itself
t = children(f);
t = [t{:}];
L = length(t);
if (L == 1)
return
end
oldt = f;
while(~strcmp(char(oldt),char(t)))
oldt = t;
t = children(t);
t = [t{:}];
t = [t{:}];
end
L = length(t);
end
edit: added desired outputs
edit2: clarified the desired function