Tools/RatioDraw.m at master · Buczman/Tools · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
function Y = RatioDraw(m1,m2,s1,s2,r,z1,z2)

% Draw empirical ratio distribution of X1, X2, s1, s2, r, on [z1 z2]

% f = @(w) ((1-r.^2).*s1.*s2.*exp(-(m1-w.*m2).^2 ./ (2.*(w.^2.*s2.^2-2.*w.*r.*s2.*s1+s1.^2)))) .* ((((2.*pi).^0.5 .* (((m1.*s2.*(w.*s2-r.*s1)+m2.*s1.*(s1-w.*r.*s2)).^2)./((1-r.^2).*(w.^2.*s2.^2-2.*w.*r.*s2.*s1+s1.^2)) ).^0.5) .* (erf (( (m1.*s2.*(w.*s2-r.*s1)+m2.*s1.*(s1-w.*r.*s2)).^2./((1-r.^2).*(w.^2.*s2.^2-2.*w.*r.*s2.*s1+s1.^2)) ).^0.5./(2.^0.5.*s1.*s2)))) ./ (s1.*s2) + (2.*exp(-(m1.*s2.*(w.*s2-r.*s1)+m2.*s1.*(s1-w.*r.*s2)).^2./(2.*(1-r.^2).*s1.^2.*s2.^2.*(w.^2.*s2.^2-2.*w.*r.*s2.*s1+s1.^2))))) ./ (2.*pi.*(w.^2.*s2.^2-2.*w.*r.*s2.*s1+s1.^2));
% f = @(w) ...
% ((m1./s1-r.*m2./s2).*w./s1 + (m2./s2-r.*m1./s1).*1./s2).*(exp(1./(2.*(1-r.^2)).*(((m1./s1-r.*m2./s2).*w./s1 + (m2./s2-r.*m1./s1).*1./s2).^2./(sqrt(w.^2./s1.^2-2.*r.*w./(s1.*s2)+1./(s2).^2)).^2-(m1.^2./s1.^2-2.*r.*m1.*m2./(s1.*s2)+m2.^2./s2.^2))))./(sqrt(2.*pi).*s1.*s2.*(sqrt(w.^2./s1.^2-2.*r.*w./(s1.*s2)+1./(s2).^2)).^3) .*...
% (normcdf(((m1./s1-r.*m2./s2).*w./s1 + (m2./s2-r.*m1./s1).*1./s2)./(sqrt(1-r.^2).*(sqrt(w.^2./s1.^2-2.*r.*w./(s1.*s2)+1./(s2).^2))),0,1) - normcdf(-((m1./s1-r.*m2./s2).*w./s1 + (m2./s2-r.*m1./s1).*1./s2)./(sqrt(1-r.^2).*(sqrt(w.^2./s1.^2-2.*r.*w./(s1.*s2)+1./(s2).^2))),0,1)) + ...
% sqrt(1-r.^2)./(pi.*s1.*s2.*(sqrt(w.^2./s1.^2-2.*r.*w./(s1.*s2)+1./(s2).^2)).^2).*...
% exp(-(m1.^2./s1.^2 - 2.*r.*m1.*m2./(s1.*s2)+ m2.^2./s2.^2)./(2.*(1-r.^2)));

f = @(w) ...
    (m1.*w./s1.^2 - r.*(m1+m2.*w)./(s1.*s2) + m2./s2.^2).*(exp(((m1.*w./s1.^2 - r.*(m1+m2.*w)./(s1.*s2) + m2./s2.^2).^2 - (m1.^2./s1.^2 - 2.*r.*m1.*m2./(s1.*s2) + m2.^2./s2.^2).*(sqrt(w.^2./s1.^2 - 2.*r.*w./(s1.*s2) + 1./s2.^2)).^2) ./ (2.*(1-r.^2).*(sqrt(w.^2./s1.^2 - 2.*r.*w./(s1.*s2) + 1./s2.^2)).^2))) ./ (sqrt(2.*pi).*s1.*s2.*(sqrt(w.^2./s1.^2 - 2.*r.*w./(s1.*s2) + 1./s2.^2)).^3).*...
    (normcdf((m1.*w./s1.^2 - r.*(m1+m2.*w)./(s1.*s2) + m2./s2.^2)./(sqrt(1-r.^2).*(sqrt(w.^2./s1.^2 - 2.*r.*w./(s1.*s2) + 1./s2.^2))),0,1) - normcdf(-(m1.*w./s1.^2 - r.*(m1+m2.*w)./(s1.*s2) + m2./s2.^2)./(sqrt(1-r.^2).*(sqrt(w.^2./s1.^2 - 2.*r.*w./(s1.*s2) + 1./s2.^2))),0,1))+...
    sqrt(1-r.^2)./(pi.*s1.*s2.*(sqrt(w.^2./s1.^2 - 2.*r.*w./(s1.*s2) + 1./s2.^2)).^2).*...
    exp(-(m1.^2./s1.^2 - 2.*r.*m1.*m2./(s1.*s2) + m2.^2./s2.^2)./(2.*(1-r.^2)));

% figure
% ezplot(f ,[-100,100]);

set(gcf,'defaultlinelinewidth',2);


if ~exist('z1','var')

    z1 = (m1-s1)/(m2-s2);
    z2 = (m1-s1)/(m2+s2);
    z3 = (m1+s1)/(m2-s2);
    z4 = (m1+s1)/(m2+s2);
    zz = [z1,z2,z3,z4];
    zmax = max(abs(zz));

    %     z1 = m1/m2 - zmax;
    %     z2 = m1/m2 + zmax;

    fplot(f, [m1/m2-zmax,m1/m2+zmax]);

else

    % z1 = RatioQuantile(m1,m2,s1,s2,r,0.025);
    % z2 = RatioQuantile(m1,m2,s1,s2,r,0.975);

    fplot(f, [z1,z2]);

end

%
%
% f = @(x) quad(@(w) myFunc(w),t0,x)-cen;
%
%
% options=optimset('Display','off');
% Y = fsolve(f,m1/m2,options);