Generate a multiplicative cascade process

2019 5 31
1 个回答
回答已采纳
更新时间：2019 6 3
4 次查看（30 天）

0 个投票

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.

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

 采纳的回答

darova
darova 2019-5-31

0 个投票

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)

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

Tchilabalo
Tchilabalo 2019-5-31
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
darova 2019-6-1
Maybe another way. Choose w and m
Then we have two equations:
Plot them with surf()
img.png
Surface intersection are values where z1==z2. Find them with contour()
img1.png
Read also about fminsearch(), fmincon() (maybe it's a better solution for this)
Read about contour() (ContourMatrix at the end)
Tchilabalo
Tchilabalo 2019-6-1
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
darova 2019-6-1
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
Tchilabalo 2019-6-1
Thanks very much for your help. It is working now.
darova
darova 2019-6-3
I made a mistake before
Z2(imag(Z2)~=0) = inf; % remove complex numbers
% should be
Z2(imag(Z2)~=0) = NaN; % remove complex numbers

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Random Number Generation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by