Generate a multiplicative cascade process

I am trying to simulate a multiplicative cascade process. At the first iteration i call four random numbers (between 0 and 1); these numbers are the initial probabilities. In the following step, i call again four random numbers and then multiply them by the parent probabilities as shown in the figure. I repeat the process 6 times. At each iteration, i want the four probabilities (x, y, z, and w) to satisfy the following conditions:

(1) x^2+y^2+z^2+w^2== m, where m is a constant

(2) x+y+z+w==1

I have tried to use the "solve" and "vapsolve" tools in matlab, but the problem is that i am getting negative numbers and zero values. I have also tried to set x and y so that i only have two unknowns, but it didn't help.

I will appreciate any help.

Answer 1

[darova]

Probably should think about initial m, w, z

m = 0.2;
z = rand;
w = rand;
syms x 
y = -x-z-w+1;
eq = x^2+y^2+z^2+w^2 -m;

X = solve(eq,x);
x = double(X)
Y = matlabFunction(y);
Y(x)

[Tchilabalo]

Thanks Darova! Your approach is efficient; with two unkowns i solved it as bellow. I find that for certain values of x, y, z, and w, the two conditions are satisfied. If you don't mind i have a follow up question: how to automatically return these numbers once the conditions are satisfied (see second part of the code).

z=rand 
w=rand
m=0.5
p=1-(w+z)% Intermediary step to solve for the root
t=m-(z^2+w^2)
y=(2*p+sqrt((-2*p)^2-4*(2*(-(p)^2-t))))/4% Root
x=1-(y+w+z)
x+y+w+z% verification
x^2/m+y^2/m+w^2/m+z^2/m
[x y z w]  


%Returning x y z w for which conditions are satisfied
function [x y z w]=prob(m)
while m>1
z=rand 
w=rand
p=1-(w+z)
t=m-(z^2+w^2)
y=(2*p+sqrt((-2*p)^2-4*(2*(-(p)^2-t))))/4
x=1-(y+w+z)
    if (x+y+w+z==1) && (x^2/m+y^2/m+w^2/m+z^2/m==1)
        [x y w z];
        return;
    end
end
end

[darova]

Maybe another way. Choose w and m

Then we have two equations:

Plot them with surf()

Surface intersection are values where z1==z2. Find them with contour()

Read also about fminsearch(), fmincon() (maybe it's a better solution for this)

Read about contour() (ContourMatrix at the end)

[Tchilabalo]

I have explored your second suggestion but i am having issues once i specify "w=rand". However, i didn't have the same problem with the previous approach. I am tempted to go with the first, but i am still not sure how to return the numbers that satisfy my condition. Below is what I have tried. Thanks for your time.

function A = prob(D)
z=rand; 
w=rand;
x=-0.02;
y=-0.03
while ((x<=0)&&(y<=0)&&(z<=0)&&(w<=0)&&((x+y+w+z)~=1)&&((x^2/m+y^2/m+w^2/m+z^2/m)~=1))
m=0.5^D;
p=1-(w+z);% Intermediary step to solve for the root
t=m-(z^2+w^2);
y=((2*p-sqrt((-2*p)^2-4*(2*(-(p)^2-t))))/4);% Root
x=(1-(y+w+z));
x+y+w+z% verification
x^2/m+y^2/m+w^2/m+z^2/m
end
A=[x,y,z,w]
end

[darova]

Are those numbers can be complex or negative?

while ((x<=0)&&(y<=0)&&(z<=0)&&(w<=0)&&((x+y+w+z)~=1)&&((x^2/m+y^2/m+w^2/m+z^2/m)~=1))    

Also you don't have to write all those conditions in one line:

tol = 1E-3;
b(1) = x<= 0;
% ...
% 0.999 ~= 1. Introduce some tolerance
b(5) = abs(1-x-y-z-w) < tol;
b(6) = abs(m - x*x-y*y-z*z-w*w) < tol;
% if don't want complex number
b(7) = abs(imag(x)) < tol;        % imiginary part smaller than tolerance
while any(b) == 0                 % do loop if one of condition is not met 

Remember also that not for any w or z exist such x and y that met the condition

k = 0;            % some counter of iterations
while ...
    if k > 1000
        disp('1000 iterations. X and y was not found')
        break;
    end
    k = k + 1;
end

[Tchilabalo]

Thanks very much for your help. It is working now.

[darova]

I made a mistake before

Z2(imag(Z2)~=0) = inf;         % remove complex numbers
    % should be
Z2(imag(Z2)~=0) = NaN;         % remove complex numbers