{VERSION 6 0 "IBM INTEL NT" "6.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 1 10 208 0 0 1 0 0 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 0 128 0 0 1 2 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 0 0 0 179 1 0 0 0 0 0 0 0 0 0 1 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "MS Sans Serif" 1 10 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 3 0 3 0 2 2 0 1 }{PSTYLE "Heading 1" -1 3 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 } 1 1 0 0 8 4 3 0 3 0 2 2 0 1 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE " " -1 -1 "Times" 1 10 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 3 0 0 0 0 3 0 3 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "Computing the cumulative b ivariate normal distribution using recursions" }}{PARA 0 "" 0 "" {TEXT -1 19 "Axel Vogt, Jan 2005" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 5 "about" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "Numerical procedures t o compute the standard bivariate normal distribution" }}{PARA 0 "" 0 " " {TEXT -1 75 "N2(x,y,rho) through recursions. Exactness should be of \+ the size of the used" }}{PARA 0 "" 0 "" {TEXT -1 71 "onedimensional cu mulative normal distribution, hence set needed Digits." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 "For 'small' correlat ion a bivariate application of the univariate series due" }}{PARA 0 " " 0 "" {TEXT -1 74 "to Marsaglia is used and for 'large' correlation a n approach of Vasicek is" }}{PARA 0 "" 0 "" {TEXT -1 63 "applied (i ha d problems with his recursion, so i build my own)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 75 "For that the standard BVN is transformed to the univariate case (cf Vasicek" }}{PARA 0 "" 0 "" {TEXT -1 78 "or Abramowitz & Stegun). That integrals have series for w hich the coefficients" }}{PARA 0 "" 0 "" {TEXT -1 71 "can be computed \+ by recursion (based on the incomplete Gamma function or" }}{PARA 0 "" 0 "" {TEXT -1 49 "looking at the integrals by partial integration)." } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 22 "The main functions are" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 74 "reduce_N2(x,y,rho) which reduces (symbolical) a bivariate normal integral," }}{PARA 0 "" 0 "" {TEXT -1 80 "(cf Abramowitz & Ste gun, \247 26.3 [especially Eq 26.3.20], p 936/937) to the cases" }} {PARA 0 "" 0 "" {TEXT -1 78 "N2(x,0,rho_new) where one variable is 0 ( by changing the correlation) and then" }}{PARA 0 "" 0 "" {TEXT -1 76 " expresses it as a one-dimensional integral N2_simplified(x,0,rho_new). It is" }}{PARA 0 "" 0 "" {TEXT -1 73 "enclosed for completeness only \+ and not needed, but may serve for testing." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 78 "N2_as_sum(x,y,rho,nSteps::posin t) which does the actual numerical computation:" }}{PARA 0 "" 0 "" {TEXT -1 85 "it first reduces to the case just described (and usually \+ splits in two integrals) and" }}{PARA 0 "" 0 "" {TEXT -1 80 "then it c alls N2_reduced_as_sum (having one variable less). This procedure care s" }}{PARA 0 "" 0 "" {TEXT -1 79 "for signs and calls N2_series. The l ater decides whether N2_Marsaglia_series or" }}{PARA 0 "" 0 "" {TEXT -1 61 "N2_Vasicek_series will be used, so these are the actual core." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 "So use \+ N2_as_sum(1.0, 2.0, 0.8, 200) to compute the bivariate normal " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 3 " " } {XPPEDIT 18 0 "int(int(exp(-1/2*(xi^2-2*rho*xi*eta+eta^2)/(1-rho^2))/s qrt(1-rho^2)/Pi/2,eta = -infinity .. y),xi = -infinity .. x);" "6#-%$i ntG6$-F$6$**-%$expG6#,$**\"\"\"F.\"\"#!\"\",(*$%#xiGF/F.**F/F.%$rhoGF. F3F.%$etaGF.F0*$F6F/F.F.,&F.F.*$F5F/F0F0F0F.-%%sqrtG6#,&F.F.*$F5F/F0F0 %#PiGF0F/F0/F6;,$%)infinityGF0%\"yG/F3;,$FCF0%\"xG" }{TEXT -1 2 " " } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 58 "for x= 1 .0, y=2.0 and rho=0.8 with at most 200 recursions." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 80 "nSteps is the maximum num ber of iterations for the recursions to be used whithin" }}{PARA 0 "" 0 "" {TEXT -1 82 "the core solutions. Since both report when they quit (according to exactness used)" }}{PARA 0 "" 0 "" {TEXT -1 82 "that nu mber gards against 'endless' loops. In the above example it will stop \+ after" }}{PARA 0 "" 0 "" {TEXT -1 37 "around 2*20 steps for 12 digits \+ with " }{XPPEDIT 18 0 ".839454198052;" "6#-%&FloatG6$\"-_!)>a%R)!#7" } {TEXT -1 11 " as result." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 85 "Proofs are not included. And i have not (yet?) don e the error estimations for cutting" }}{PARA 0 "" 0 "" {TEXT -1 75 "se ries to sums (depending on the used recursion steps, Vasicek has it) . .. " }}{PARA 0 "" 0 "" {TEXT -1 88 "\nAnd for the 'Marsaglia series' o ne has probably not to reduce to the 'simplified' case," }}{PARA 0 "" 0 "" {TEXT -1 57 "but can recurse directly (for correlation not too la rge)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 " At maple-new http://groups.yahoo.com/group/maple-new/message/212 i ask ed for" }}{PARA 0 "" 0 "" {TEXT -1 81 "help having trouble with a cras h. The two case are discussed as extreme examples," }}{PARA 0 "" 0 "" {TEXT -1 76 "and have been my motivation to write this sheet from snip pets laying around." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 47 "Finally the limiting case |rho| = 1 is treated." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 93 "As i do n ot have Fortran and high precision integrators ... i hope i am not hea vily wrong ..." }}{PARA 0 "" 0 "" {TEXT -1 35 "any comments and hints \+ are welcome." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 11 "References:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 190 "George Marsaglia, Evalua ting the Normal Distribution, Journal of Statistical Software (2004), \nhttp://www.jstatsoft.org/counter.php?id=100&url=v11/i04/v11i04.pdf&c t=1 \nfor the univariate case" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 43 "Oldrich Alfons Vasicek (1998), Moody's KM V," }}{PARA 0 "" 0 "" {TEXT -1 101 "http://www.moodyskmv.com/research/ whitepaper/A_Series_Expansion_for_the_Bivariate_Normal_Integral.pdf" } }{PARA 0 "" 0 "" {TEXT -1 77 "as series in 1-rho^2 (which is like Drez ner & Wesolowsky, compare Genz below)" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 23 "Abramowitz & Stegun ..." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 8 "Related:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 57 "Alan Miller 'fo r classical' Fortran Code refering to Genz" }}{PARA 0 "" 0 "" {TEXT -1 36 "http://users.bigpond.net.au/amiller/" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 61 "Alan Genz, Numerical Computatio n of Rectangular Bivariate ..." }}{PARA 0 "" 0 "" {TEXT -1 49 "http:// www.sci.wsu.edu/math/faculty/genz/homepage" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "This code works, it is fast and good for ~ 14 - 15 decimal places." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "restart; kernelopts(versi on);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%VMaple~9.52,~IBM~INTEL~NT,~De c~17~2004~Build~ID~175529G" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 122 "#signum(0); signum(0.); signum(-0.0); evalf(%);\n#eval(_Envsignum 0); anames('environment');\n#_Envsignum0 := 0:\nDigits:= 14;" }}} {SECT 1 {PARA 3 "" 0 "" {TEXT -1 25 "Notations and conventions" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "pdfN:= x -> exp(-x^2/2)/sqrt (2*Pi); \nN:= x -> (1+erf(x/sqrt(2)))/2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 225 "pdfN2:= (x,y,rho) -> 1/sqrt(1-rho^2)/(2*Pi)*exp(-(x^ 2-2*rho*x*y+y^2)/(2*(1-rho^2)));\n``;\nN2:= (x,y,rho) -> \n 1/sqrt(1- rho^2)/2/Pi*Int(Int(exp(-1/2*(xi^2-2*rho*xi*eta+eta^2)/((1-rho^2))), \+ eta=-infinity..y),xi=-infinity..x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 141 "N2_simplified:= (x,y,rho) -> Int(pdfN(xi)*N((y-rho*x i)/sqrt(1-rho^2)),xi=-infinity..x);\n#``;\n#N2_reduced:= (x,rho) -> N2 _simplified(x,0,rho);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 10 "Note that " }{XPPEDIT 18 0 "N2(x,y,rho) = N2_simplified(x,y,rho);" "6#/-%#N2G6%%\"xG%\"yG%$rhoG-%.N2_simplified G6%F'F(F)" }{TEXT -1 14 " (reference?) " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "assume(-1 " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 4 "code" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 874 "reduce_N2:=proc(x,y,rho)\nl ocal r,rho_x,rho_y,result;\n\ndescription \"reduce cdfN2 to N2_simplif ied(x,0,rho_new) by A&S 26.3.20\";\n\n#if 1 " 0 "" {MPLTEXT 1 0 1102 "N2_as_sum:=proc(x,y,rho,nSteps::posint)\n# decompos e to get the reduced case y=0 by changing rho\nlocal r,rho_x,rho_y,res ult;\n\ndescription \"reduce cdfN2 to N2_simplified(x,0,rho_new) by A& S 26.3.20\";\n\nif 1 < abs(rho) then\n# return FAIL;\nelif abs(rho) = 1 then\n# return evalf( N2_absRhoEqualOne(x,y,rho) );\nelse\nend if; \n\nif signum(rho)=0 then\n return evalf(N(x)*N(y));\nelif (signum(x) =0 and signum(y)=0) then\n return evalf(1/4*(Pi+2*arctan(rho/(1-rho^2 )^(1/2)))/Pi);\nelse\nend if;\n\nif signum(x)=0 then\n return evalf(N 2_reduced_as_sum(y,rho,nSteps));\nelif signum(y)=0 then\n return eval f(N2_reduced_as_sum(x,rho,nSteps));\nelse\nend if;\n\nr:= ( 1/sqrt(x^2 -2*rho*x*y+y^2));\nrho_x:=( (rho*x-y)*r*signum(x) );\nrho_y:=( (rho*y- x)*r*signum(y) );\n\nif signum(x*y) = -1 then\n result:=N2_reduced_as _sum(x,rho_x,nSteps)+\n N2_reduced_as_sum(y,rho_y,nSteps) - 1/2;\n# print(rho_x, rho_y);\nelif signum(x*y) = 1 then\n result:=N2_reduced _as_sum(x,rho_x,nSteps)+\n N2_reduced_as_sum(y,rho_y,nSteps);\n#pri nt(rho_x, rho_y);\nelse\n result:=reduce_N2(x,y,rho);\nend if;\nretur n evalf(result);\nend proc; maplemint(%); ``;\n" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 976 "N2_reduced_as_sum:=proc(x,r,nSteps::posint)\n # reduce to positive variables (needed only for Vasicek)\n# assumed re al variables\n\ndescription \"prepare to call series computation with \+ positive arguments\";\n \nif type(x,numeric) and type(r,numeric) then \n # ok \nelse\n # exit with a formal integral\n return N2_simplifi ed(x,0,r);\nend if;\n\nif 1<=abs(r) then\n return N2_simplified(x,0,r );\nend if;\n\n# trivial cases first, signum(0)=0 assumed\nif signum(x ) = 0 then\n return evalf( 1/4+1/2*1/Pi*arctan(r/(1-r^2)^(1/2)) );\ne nd if;\nif signum(r) = 0 then\n return evalf( N(x)/2 );\nend if;\n\ni f (0 " 0 "" {MPLTEXT 1 0 386 "N2_series:=proc(x,r,nSteps: :posint)\n# now everthing is a float and strictly positive\n\ndescript ion \"decide which series is used (with positive arguments)\";\n\nif a bs(r) <= evalf(1/sqrt(2)) then\n return N2_Marsaglia_series(x,-r,nSte ps); # sign changes !\nelif abs(r) < 1.0 then\n return N2_Vasicek_ser ies(x,r,nSteps);\nelse\n return N2_simplified(x,0,r);\nend if;\nend p roc; maplemint(%); ``;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 871 "N2_Vasicek_series:=proc(x,rho,nSteps::posint)\nlocal beta0,A,B,k, sumA;\n\ndescription \"cdfN2(x,0,rho) through int(pdfN2(x,0,r,r=rho..1 ) by recursion using nSteps\";\n\nDigits:=Digits+3;\n\nbeta0 := (1-rho ^2)^(1/2)*signum(x)*exp(1/2/(rho+1)/(rho-1)*x^2);\n\nsumA:=0.0;\nA:= e valf(- (x*2*sqrt(2*Pi)*N(-abs(x)/(1-rho^2)^(1/2)) -2*beta0)); # A(0)\n #A:= evalf(- ((x)*2*sqrt(2*Pi)*N(-(x)/(1-rho^2)^(1/2)) -2*beta0)); # A (0)\nsumA:=evalf(sumA+A);\n\nB:= evalf(- 2*beta0*(1-rho^2)); # B(1)\nA := evalf(-1/6*x^2*A-1/6*B); # A(1)\nsumA:=evalf(sumA+A);\n\nfor k from 2 to nSteps do\n B:= evalf(-1/2*(-1+rho^2)*(2*k-1)/(k-1)*B); #\n A: = evalf((x^2*A*(-2*k+1)-B)*1/2/k/(1+2*k)); #\n sumA:=evalf(sumA+A);\n if sumA = sumA - A then break; end if; # sum does not change anymore \nend do;\n\n#return evalf(sumA/(4*Pi));\nprint('iterations_Vasicek'=k );\nreturn 0.5 - evalf(sumA/(4*Pi));\nend proc; maplemint(%); ``;\n" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 751 "N2_Marsaglia_series:=proc (x::numeric,rho::numeric,nSteps::posint)\nlocal a::numeric, i::integer ,A::numeric,B::numeric,sumA::numeric;\n\ndescription \"int(pdfN(xi)*cd fN(rho*xi),xi=-infinity..x) by recursion using nSteps\";\n\nDigits:=Di gits+3;\n\na:=evalf(rho/sqrt(1-rho^2)); # sign changes !\n\nA:=0; B:=0 ; sumA:=0;\n\nA:=evalf(-1/2*1/(1+a^2)*exp(-1/2*(1+a^2)*x^2)/Pi*a); # A (1)\nB:=evalf( (a*x)^2 * A); # B(1)\n\nsumA:=evalf(A);\nfor i from 1 t o nSteps-1 do\n#print(A);\n A:=evalf( (2*i*a^2/(1+a^2)*A + B)/(1+2*i) ); # A(n+1)\n B:=evalf( (a*x)^2 * B /(1+2*i)); # B(n+1)\n sumA:=eval f(sumA + A);\n if sumA = sumA - A then break; end if; # sum does not \+ change anymore\nend do; \n\nprint('iterations_Marsaglia'=i);\nreturn s umA + 0.5*evalf(N(x));\nend proc;: maplemint(%); ``;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 8 "Examples" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 52 "Set y=0 for first tests and call t he 'core' directly" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 201 "xTst :=-0.20; rhoTst:=0.5; nTst:=30;\n``;\n'N2_Marsaglia_series(xTst,-rhoTs t,nTst)': '%'=%; # sign changes !\n'N2_reduced(xTst,rhoTst)': '%'= eva lf(%);\n'eval(N2(xTst,0,rhoTst),infinity=80)': '%'= evalf(%);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "Do the same, but with high correla tion" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 184 "xTst:=-3.20; rhoTs t:=0.9; nTst:=30;\n``;\n'N2_Marsaglia_series(xTst,-rhoTst,nTst)': '%'= %;\n'N2_reduced(xTst,rhoTst)': '%'= evalf(%);\n'eval(N2(xTst,0,rhoTst) ,infinity=80)': '%'= evalf(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "So increase number of steps and compare the two methods" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 187 "xTst:=-3.20; yTst:=0; rhoTst:=0.9; nTst:=300;\n``;\n'N2_Marsaglia_series(xTst,-rhoTst,nTst)': '%'=%;\n'N 2_Vasicek_series(xTst,rhoTst,nTst)': '%'=%;\n'N2_as_sum(xTst,yTst,rhoT st,nTst)': '%'=%;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 16 "which shows \+ that" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 91 " - Marsaglia will converge after more iterations (196, but have patient s for very high rho)," }}{PARA 0 "" 0 "" {TEXT -1 81 "- it is neccessa ry to care for signs (if calling directly the 'Vasicek' routine)," }} {PARA 0 "" 0 "" {TEXT -1 26 "- and combining is better\n" }}{PARA 0 " " 0 "" {TEXT -1 63 "And just that is, was the procedures around the tw o methods do." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 139 "xTst:=-1. 20; yTst:=1.7; rhoTst:=0.9; nTst:=30;\n``;\n'N2_as_sum(xTst,yTst,rhoTs t,nTst)': '%'=%;\n'reduce_N2(xTst,yTst,rhoTst)': '%'= evalf(%);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "Now take extreme variables at an a verage correlation:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 194 "xTs t:=0.001; yTst:=5; rhoTst:=0.5; nTst:=30;\n``;\n'N2_as_sum(xTst,yTst,r hoTst,nTst)': '%'=%;\n'reduce_N2(xTst,yTst,rhoTst)': '%'= evalf(%);\n' eval(N2(xTst,yTst,rhoTst),infinity=80)': '%'= evalf(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 46 "This is somewhat strange: greater 1/2 or \+ not? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 80 " Actually it is above 1/2 i think, but Maple needs too much time to com pute that." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 79 "And it may be worth to think whether it is really a good to red uce the integral" }}{PARA 0 "" 0 "" {TEXT -1 24 "allways in that form \+ ..." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "Run through (some) combina tions [and generating by random would be better]:" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 424 "remDigits:=Digits: Digits:=14:\nyTst:=2; nT st:=40;\nfor xTst in \{1.0\} do #\{0.01, 0.1, 1.0, 3.0\} do # \nfor rT st in \{0.5\} do #\{0.7, 0.8, 0.9, 0.95, 0.99\} do # \n\nfor sx in \{- 1,+1\} do\n for sr in \{-1,+1\} do\n print(sx,sr);\n `error`=N2 _as_sum(sx*xTst,yTst,sr*rhoTst,nTst)-evalf(reduce_N2(sx*xTst,yTst,sr*r hoTst));\n print(%);\n end do; # sr \nend do; # sx\nsx:='sx': sr:= 'sr':\n\nend do; # rTst\nend do; # xTst\nDigits:=remDigits:" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 106 "For extreme correlation use high \+ numbers of steps and look, what happens close\nbelow the the needed st eps:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 133 "remDigits:=Digits: Digits:=100:\n'N2_as_sum(1,2,0.000000001,1000)': '%'=%;\n'N2_as_sum(1 ,2,0.000000001,100)': '%'=%;\nDigits:=remDigits:" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 78 "Compare reported steps and use some more steps (of course it does not change):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 133 "remDigits:=Digits: Digits:=100:\n'N2_as_sum(1,2,0.999999999,100 0)': '%'=%;\n'N2_as_sum(1,2,0.999999999,100)': '%'=%;\nDigits:=remDigi ts:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 15 "ex treme example" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 167 "Now attack my qu estions posted at maple-new on 29 Dec 2004, \nhttp://groups.yahoo.com/ group/maple-new/message/212\nwhere i got answers from Robert Israel an d Chris 'CW'." }}{PARA 0 "" 0 "" {TEXT -1 75 "\nThat was to compute In t(pdfN(x)*N((b-rho*x)/sqrt(1-rho^2)),x=-infinity..a)" }}{PARA 0 "" 0 " " {TEXT -1 71 "for a=b=2 and rho = 0.999999999 or rho= - 0.999999999 ( 9 times the 9)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 74 "Using the notation above this is N2_simplified(a,b,rho) a nd restricting to" }}{PARA 0 "" 0 "" {TEXT -1 76 "14 digits Maple comp utes it (via the NAG library?) but has difficulties with" }}{PARA 0 " " 0 "" {TEXT -1 17 "higher exactness." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 35 "Chris CW reports (with 128 digits):" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 ".97724890478596692786 1224346192468443840763758605052838842081675512696805867075671095742423 95711932276563251547307155689221054584;" "6#-%&FloatG6$\"[s%ea5A*obrIZ :DjlFK>r&RUUd4rc2ne!op7bn\"3U)QG00'ePwSQWoC>YVAhy#p'fy/*[s(*!$G\"" } {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 66 "and N2_as_sum would be t he same (128 digits and after 2*15 steps) " }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 ".97724890478596692786122434619246844384076375 8605052838842081675512696805867075671095742423957119322765632515473071 55689221054584;" "6#-%&FloatG6$\"[s%ea5A*obrIZ:DjlFK>r&RUUd4rc2ne!op7b n\"3U)QG00'ePwSQWoC>YVAhy#p'fy/*[s(*!$G\"" }{TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 71 "Chris CW (with 96 digits, not sure that i understood his answer right):" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{XPPEDIT 18 0 ".954494886386402429208688647347934 284261823176219001417827270552206855420921731296582001085575261;" "6#- %&FloatG6$\"[qh_d&3,?e'HJ<#4Ubo?_0FFyT,!>iu/R(**zmK?8V mf_vWc]7LpmDZV*f&eTO5O(*\\a*!#'*" }{TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "In the following i want t o know, which one is 'correct'." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 104 "'N2_simplified(a,b,rho)'=\n 'Int(pdfN(x)*N((b-rho*x)/sqrt(1-rho^2)),x=-infinity..a)';\nsubs(x=xi,% ): is(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 59 "The above numerical \+ procedures can compute it for both rho:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 150 "remDigits:=Digits: Digits:=100:\n'N2_as_sum(2,2,+0.9 99999999,100)'; evalf(%,96);\n``;\n'N2_as_sum(2,2,-0.999999999,100)'; \+ evalf(%,96);\nDigits:=remDigits:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "The second result differs from the value reported by Chris, which \+ was" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 99 "0.954494886386402429 2086886473479342842618231762190014178272705522068554209217312965820010 85575261;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 47 "For the rest let us take a closer look at this." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 150 "First \+ compare, what a direct use of numerical integration for the standard ( double integral)\nand reduced version would give for that data at 14 d igits" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 138 "rTst:=-0.9;\n'N2_ simplified(2,2,rTst)': '%'=evalf(%);\n'reduce_N2(2,2,rTst)': '%'=evalf (%);\n'N2(2,2,rTst)': '%'= evalf(subs(infinity=40,%));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "rTst:=-0.999999999;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 108 "'N2_simplified(2,2,rTst)': '%'=eva lf(%);\n'reduce_N2(2,2,rTst)': '%'=evalf(%);\n'N2(2,2,rTst)': '%'= eva lf(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 90 "So for usual exactness \+ that is ok (of course), if one does not work with double integrals," } }{PARA 0 "" 0 "" {TEXT -1 70 "but with the reduced form (but may be no t without troubles, cf above)." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "Estimate, where to cut off integration 'certainly':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 181 "'reduce_N2(2,2,rho)': '%'= simplify(%); \n\nremDigits:=Digits: Digits:=100:\nmyProblem:='reduce_N2(2,2,rTst)'; :\nDigits:=14:\n`myProblem at 14 digits`=evalf(myProblem);\nDigits:=re mDigits:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 238 "remDigits:=Dig its: Digits:=80:\n#'integrand(myProblem/2)'; %;\n#plot(%,xi=-3..2);\n' fsolve(integrand(myProblem/2)=1e-1000, xi)': '%'=%;\nrhs(%):\ncutoff:= max(70,%);\nmyProblem_cutted_off:='(subs(infinity=cutoff,myProblem))'; \nDigits:=remDigits:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "and chec k, what the integral would be with the cut off used (at 45 digits):" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 115 "remDigits:=Digits: Digits :=35:\n'N2_as_sum(2,2,rTst,100)'; %;\n``;\nmyProblem_cutted_off; evalf (%);\nDigits:=remDigits:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "Hence for 35 digits it is 'ok' and for more i lost patients with Maple." }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 "B ut then i wanted to know it ..." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "Let us see, why Maple has troubles with the chosen parameter:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 115 "'rTst': '%'=%;\n'eval(integrand(N2_simplified(2,2,rho)), rho=rTst )';\nplot(%,xi=-3..2);\nplot(%%,xi=-2-1e-3..-2+1e-3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "As i asserted (and use, but did not prove explicitely): N2(2,2,rho) = reduce_N2(2,2,rho)" }}{PARA 0 "" 0 "" {TEXT -1 59 "and i want to use the later one, calling the integrand A_ r:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "''reduce_N2(2,2,rho)/ 2'': %= simplify(eval(%));\n rhs(%):\n A_r:= integrand(%);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "but this itself is a N2_simplified " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "''N2_simplified(2,0,-be ta)'': %= simplify(eval(%));\n rhs(%):\n A_s:= integrand(%);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "and calling the integrand A_s one \+ determines the beta so they become equal:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 133 "#A_r/A_s; simplify(%);\n(xi*(1-rho)^(1/2)/(2*rho+2 )^(1/2))=(1/2*beta*xi/(1-beta^2)^(1/2)*2^(1/2));\n%/xi: [solve(%,beta) ];\nbeta1:=%[1];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "check, that w ith " }{XPPEDIT 18 0 "beta1;" "6#%&beta1G" }{TEXT -1 29 " the desired \+ indentity holds:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "'reduce _N2(2,2,rho)'-'2*N2_simplified(2,0,-beta1)'=0;\nsimplify(%): is(%);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "A_t:='eval(A_s,beta=beta1) '; simplify(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 66 "Let us see, wh ere Maple may have troubles for the translated task:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "'eval(A_r,[rho=rTst])';\nplot(%,xi=-3..2) ;\nplot(%%,xi=-1e-4..1e-4);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "So i split integration at 0" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "2*(N2_simplified(0,0,-beta) + Int(1/2*exp(-1/2*xi^2)*2^(1/2)/Pi^( 1/2)*(1/2+1/2*erf(1/2*beta*xi/(1-beta^2)^(1/2)*2^(1/2))),xi = 0 .. 2)) ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 19 "The first integral " } {XPPEDIT 18 0 "N2_simplified(0,0,-beta);" "6#-%.N2_simplifiedG6%\"\"!F &,$%%betaG!\"\"" }{TEXT -1 10 " writes as" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 67 "reduce_N2(0,0,-beta);\neval(%,beta=beta1): \nlowerP art:=simplify(%);;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 52 "while the s econd integral now is evaluated by Maple:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 154 "'eval(Int(1/2*exp(-1/2*xi^2)*2^(1/2)/Pi^(1/2)*(1/2 +1/2*erf(1/2*beta*xi/(1-beta^2)^(1/2)*2^(1/2))),xi = 0 .. 2), beta=bet a1)': \nremainingPart:=simplify(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "hoping, that Maple can do it easily - which is not quite the ca se, " }}{PARA 0 "" 0 "" {TEXT -1 26 "but splitting again helps:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 192 "remDigits:=Digits: Digits:= 96:\n'eval(1/4*2^(1/2)/Pi^(1/2)*Int(exp(-1/2*xi^2)*(1+erf(1/2*(-2*rho+ 2)^(1/2)*xi/(rho+1)^(1/2))),xi = 0 .. 1e-2),rho=rTst)';\nremaining1:=e valf(%);\nDigits:=remDigits:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 89 "c alling Maple for the below with that much digits crashes with an error in \"LIBGMP-3.DLL\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 230 "re mDigits:=Digits: Digits:=96:\n'eval(1/4*2^(1/2)/Pi^(1/2)*Int(exp(-1/2* xi^2)*(1+erf(1/2*(-2*rho+2)^(1/2)*xi/(rho+1)^(1/2))),xi = 1e-2..2),rh o=rTst)';\n%;:\n#remaining2:=evalf(%);\n#erf(31622.776593778099168*0.0 1);\nDigits:=remDigits:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 89 "but th e argument in erf approaches 1 with increasing xi and even for 10000 d igits it 'is'" }}{PARA 0 "" 0 "" {TEXT -1 91 "already 1 for xi=0.01. S o i set it to 1. And only need an integral, which is the difference" } }{PARA 0 "" 0 "" {TEXT -1 19 "of error functions:" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 86 "1/4*2^(1/2)/Pi^(1/2)*Int(exp(-1/2*xi^2)*(1+1 ),xi = 1e-2..2);\nremaining2:=evalf(%,96);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 59 "Now add the things and compare with the recursive proce dure" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 157 "remDigits:=Digits: Digits:=96:\n'2*(eval(lowerPart,rho=rTst)+remaining1+remaining2)'; ev alf(%);\n``;\n'N2_as_sum(2,2,rTst,100)'; evalf(%,96);\nDigits:=remDigi ts:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 24 "which makes me happy :-)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 86 "At least almost ... taking -1 instead of -0.999999999 and using A&S 26.3 .16 the result" }}{PARA 0 "" 0 "" {TEXT -1 91 "should be N(-2) - N(2) \+ = - erf( sqrt(2) ), which is the above - but we have opposite sign: " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "eval(N(-x)-N(y), [x=2,y=2]); evalf(%,96); #evalf(-%); identify (%);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 24 " appendix: abs(rho) = 1 " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "I thin k the sign in A&S 26.3.16 actually is false, it has to be (-1) times t heir result," }}{PARA 0 "" 0 "" {TEXT -1 75 "the 3 other cases are cor rect (but my proofs are too ugly to include them):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 332 "N2_absRhoEqualOne:=proc(x,y,rho)\n# limiti ng case abs(rho) = 1 for N2(x,y,rho)\n\nif rho = 1 then\n if x <= y t hen\n return N(x)\n elif y <= x then\n return N(y)\n else\n e nd if;\nelif rho = - 1 then\n if -y <= x then\n return N(x)+N(y)-1 \n elif x <= -y then\n return 0\n else\n end if;\nelse\nend if; \n\nend proc; maplemint(%); ``;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "Abramowitz & Stegun is as follows:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 486 "N2_absRhoEqualOne_AS:=proc(x,y,rho)\n# A&S 26.3.14 f f, N2(x,y,rho) = L(-x,-y,rho), P(t) = N(t), Q(t) = N(-t)\n\nif rho = - 1 then\n if -x-y >= 0 then # 26.3.16\n return 0\n else \+ # 26.3.16\n return N(-x)-N(-(-y)) # = P(-x) - Q( -y)\n end if;\nelif rho = 1 then\n if -y <= -x then # 26.3. 17\n return N(-(-x)) # = Q(-x)\n else # \+ 26.3.18\n return N(-(y)) # = Q(-y)\n end if;\nelse\nend if ;\n\nend proc; maplemint(%); ``;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "compare that procedures using some test data" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 268 "xTst:=1; yTst:='xTst'; rTst:=1-(1e-9); nStep sTst:=100:\n``;\nremDigits:=Digits: Digits:=20:\n'N2_as_sum(xTst,yTst, rTst, nStepsTst)'; %;\n'N2_absRhoEqualOne(xTst,yTst,signum(rTst))'; ev alf(%); ``;\n'N2_absRhoEqualOne_AS(xTst,yTst,signum(rTst))'; evalf(%); \nDigits:=remDigits:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 269 "xT st:=-1; yTst:='xTst'; rTst:=1-(1e-9); nStepsTst:=100:\n``;\nremDigits: =Digits: Digits:=20:\n'N2_as_sum(xTst,yTst,rTst, nStepsTst)'; %;\n'N2_ absRhoEqualOne(xTst,yTst,signum(rTst))'; evalf(%); ``;\n'N2_absRhoEqua lOne_AS(xTst,yTst,signum(rTst))'; evalf(%);\nDigits:=remDigits:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 272 "xTst:=-1; yTst:='xTst'; rTs t:=-(1-(1e-9)); nStepsTst:=100:\n``;\nremDigits:=Digits: Digits:=20:\n 'N2_as_sum(xTst,yTst,rTst, nStepsTst)'; %;\n'N2_absRhoEqualOne(xTst,yT st,signum(rTst))'; evalf(%); ``;\n'N2_absRhoEqualOne_AS(xTst,yTst,sign um(rTst))'; evalf(%);\nDigits:=remDigits:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 318 "xTst:=1; yTst:='xTst'; rTst:=-(1-(1e-2)); nStepsTs t:=100:\n``;\nremDigits:=Digits: Digits:=20:\n'N2_as_sum(xTst,yTst,rTs t, nStepsTst)'; %; ``;\n'reduce_N2(xTst,yTst,rTst)'; evalf(%); ``;\n'N 2_absRhoEqualOne(xTst,yTst,signum(rTst))'; evalf(%); ``;\n'N2_absRhoEq ualOne_AS(xTst,yTst,signum(rTst))'; evalf(%);\nDigits:=remDigits:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}} {MARK "2 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }