I have a non-linear system of equations in Matlab and I would like to find the set of its solutions. The system is represented below. 
Imagine for a moment that the Matlab function mvncdf can be used together with solve. My system of equations looks like this
clear %SOME VALUES USED BELOW U1true=1; U2true=-1; U3true=-2; U4true=1; rhotrue=-0.5; mutrue=[0 0]; sigmatrue=[2*(1-rhotrue) 1-rhotrue; 1-rhotrue 2*(1-rhotrue)]; %SET VARIABLES syms U1 U2 U3 U4 rho %SOLVE SYSTEM WITH 6 EQUATIONS S=solve(mvncdf([ U1 -U2+U1], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U1true -U2true+U1true], mutrue, sigmatrue) ,... mvncdf([ U2 U2-U1], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U2true U2true-U1true], mutrue, sigmatrue),... mvncdf([-U1 -U2], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([-U1true -U2true], mutrue, sigmatrue),... mvncdf([ U3 -U4+U3], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U3true -U4true+U3true], mutrue, sigmatrue) ,... mvncdf([ U4 U4-U3], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U4true U4true-U3true], mutrue, sigmatrue),... mvncdf([-U3 -U4], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([-U3true -U4true], mutrue, sigmatrue)); %DISPLAY SOLUTIONS M=[S.U1,S.U2,S.U3,S.U4, S.rho]; This approach clearly does not work for many reasons, most importantly because mvncdf cannot be used together with solve, as discussed here for example.
As an alternative method to solve my system, I thought about using slicesample as follows
This code implements my idea.
%MAIN clear rng default Utrue=[1;-1;-2;1]; rhotrue=-0.5; number = 100000; initials=[Utrue;rhotrue]; %[U1;U2;U3;U4;rho] rnd = slicesample(initials,number,'pdf',@dens); %ANCILLARY FUNCTIONS TO SAVE function density=dens(param) T=10^(-6); if param(end)>=1 || param(end)<=-1 %rho outside (-1,1) density=0; %force density to be equal to zero else Utrue=[1;-1;-2;1]; rhotrue=-0.5; mutrue=[0 0]; sigmatrue=[2*(1-rhotrue) 1-rhotrue; 1-rhotrue 2*(1-rhotrue)]; truth=mvncdf([ Utrue(1) -Utrue(2)+Utrue(1);... Utrue(2) Utrue(2)-Utrue(1);... -Utrue(1) -Utrue(2) ;... Utrue(3) -Utrue(4)+Utrue(3);... Utrue(4) Utrue(4)-Utrue(3);... -Utrue(3) -Utrue(4)], mutrue, sigmatrue); density=exp(-criterion(param,truth)/T); end end function Q=criterion(param, truth) U=param(1:end-1); rho=param(end); mu=[0,0]; sigma=[2*(1-rho) 1-rho; 1-rho 2*(1-rho)]; candidates=mvncdf([ U(1) -U(2)+U(1);... U(2) U(2)-U(1);... -U(1) -U(2) ;... U(3) -U(4)+U(3);... U(4) U(4)-U(3);... -U(3) -U(4)], mu, sigma); Q=(candidates(1)-truth(1))^2+(candidates(2)-truth(2))^2+(candidates(3)-truth(3))^2+... (candidates(4)-truth(4))^2+(candidates(5)-truth(5))^2+(candidates(6)-truth(6))^2; %this is zero at a solution of the system end Question: using the approach that I propose, as I increase the variable number I get more and more solutions for different values of rho, which is in line with comment above that the system has one and only one solution with respect to U for every value of rho. However, to span the entire (-1,1) interval for rho it seems to me that I would need to set the variable number immensely high. Do you have (i) suggestions that can help me to improve the code above to solve this issue, or (ii) better approaches to provide the set of solutions of my system of equations?
0 comments:
Post a Comment