{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 } {CSTYLE "" -1 256 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE " " -1 257 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE " " -1 -1 "Helvetica" 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 3" -1 5 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 0 0 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 }{PSTYLE "Maple Outpu t" -1 12 1 {CSTYLE "" -1 -1 "Times" 1 10 0 0 0 1 2 2 2 2 2 2 1 1 1 1 } 1 3 0 0 0 0 3 0 3 0 2 2 0 1 }{PSTYLE "Maple Plot" -1 13 1 {CSTYLE "" -1 -1 "Helvetica" 1 10 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 3 0 3 0 2 2 0 1 }{PSTYLE "Normal" -1 256 1 {CSTYLE "" -1 -1 "Courier" 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 }} {SECT 0 {EXCHG {PARA 5 "" 0 "" {TEXT -1 55 "The Cumulative Normal Dist ribution for Dimensions up to" }}{PARA 5 "" 0 "" {TEXT -1 57 "3 using \+ the qfloat floating-point Library from LCC-WIN32:" }}{PARA 5 "" 0 "" {TEXT -1 28 "Some Tests for Dimension = 3" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "This implementation gives exactnes s over almost the 104 digits which the library provides." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 19 "The system library " }{TEXT 257 10 "qfloat.dll" }{TEXT -1 58 " should be in Window's sys tem directory and the concurrent" }}{PARA 0 "" 0 "" {TEXT 256 11 "cdfn 123.dll" }{TEXT -1 53 " should be placed in the directory of this work sheet." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 13 "AVt, Jan 2006" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "restart;\nk ernelopts(version);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%VMaple~10.02,~ IBM~INTEL~NT,~Nov~8~2005~Build~ID~208934G" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 16 "Digits_lcc:=105;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #>%+Digits_lccG\"$0\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "Di gits:=2*Digits_lcc; # greater precision to check results" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%'DigitsG\"$5#" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "For using the DLL locate \+ its directory to call external functions from there:" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "currentdi r(): myDLL:=cat(%,`\\\\cdfn123.dll`);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&myDLLGQ_pC:\\_Work\\Maple_Work\\_allTheStuff\\Statistik\\Verteil ungen\\TrivariateNormal\\cdfn123.dll6\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 47 "Accessing the DLL functio ns is through strings:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 98 "lccstr:=proc(var)\n convert( evalf (parse(convert(var,string)),Digits_lcc + 10), string);\nend proc:" }}} {SECT 1 {PARA 5 "" 0 "" {TEXT -1 34 "The cumulative normal distributio n" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "Define the cumulative normal distribution within Maple" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 91 "cdf N := x -> 1/2+1/2*erf(1/2*x*2^(1/2));\npdfN := x -> 1/2*1/Pi^(1/2)*exp (-1/2*x^2)*2^(1/2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%cdfNGf*6#%\" xG6\"6$%)operatorG%&arrowGF(,&#\"\"\"\"\"#F.*&F-F.-%$erfG6#,$*&F-F.*&9 $F.F/F-F.F.F.F.F(F(F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%pdfNGf*6#% \"xG6\"6$%)operatorG%&arrowGF(,$*&#\"\"\"\"\"#F/*(%#PiG#!\"\"F0-%$expG 6#,$*&#F/F0F/*$)9$F0F/F/F4F/F0F.F/F/F(F(F(" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 411 "fct_cdfN := define_external(\n 'cdfN_string_Maple ',\n 'C',\n 'x_str'::string[],\n 'ret_str'::string[],\n RETURN::in teger[4],\n LIB=myDLL):\n\ncdfN_lcc:=proc(x)\n local X::string, resu lt::string;\n result:=StringTools:-Fill( `0` , Digits_lcc+10);\n X:= lccstr(x): \n if type(parse(X),numeric) then\n fct_cdfN(X,result); \n return parse(result);\n else\n return 'cdfN_lcc(x)';\n end \+ if;\nend proc: #maplemint(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 102 "This will provide the DLL with memor y space (given as a string y_str) to store the results in the DLL." }} {PARA 0 "" 0 "" {TEXT -1 103 "Since update is 'inplace' this will modi fy the string and its length. As Maple is a symbolic system one" }} {PARA 0 "" 0 "" {TEXT -1 111 "should never call this result directly, \+ since this inconsistency for the same object will crash it (just try i t" }}{PARA 0 "" 0 "" {TEXT -1 61 "and restart ...). So a simple proced ure is used as interface." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 18 "The Bivariate Case" }}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 52 "The bivariate norm al distribution can be written as:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 181 "pdfN2:= (x,y,rho) -> \n 1 /sqrt(1-rho^2)/(2*Pi)*exp(-(x^2-2*rho*x*y+y^2)/(2*(1-rho^2)));\n``;\nc dfN2:= (x,y,rho) -> \n Int(Int( pdfN2(xi,eta,rho), eta=-infinity..y ),xi=-infinity..x);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&pdfN2Gf*6%% \"xG%\"yG%$rhoG6\"6$%)operatorG%&arrowGF*,$*&#\"\"\"\"\"#F1*(-%%sqrtG6 #,&F1F1*$)9&F2F1!\"\"F;%#PiGF;-%$expG6#,$*&,(*$)9$F2F1F1**F2F1F:F1FEF1 9%F1F;*$)FGF2F1F1F1,&F2F1*&F2F1F9F1F;F;F;F1F1F1F*F*F*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&cdfN2Gf* 6%%\"xG%\"yG%$rhoG6\"6$%)operatorG%&arrowGF*-%$IntG6$-F/6$-%&pdfN2G6%% #xiG%$etaG9&/F7;,$%)infinityG!\"\"9%/F6;F;9$F*F*F*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 97 "It is coded withi n the DLL giving 104 decimal points of precision and can be accessed a s follows:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 559 "fct_cdfN2 := define_external(\n 'cdfN2_string_Map le',\n 'C',\n 'x_str'::string[], \n 'y_str'::string[],\n 'r_str':: string[],\n 'ret_str'::string[],\n RETURN::integer[4],\n LIB=myDLL) :\n\ncdfN2_lcc:=proc(x,y,r)\n local X::string, Y::string, R::string, \+ result;\n\n result:=StringTools:-Fill( `0` , Digits_lcc+10);\n X:=lc cstr(x); \n Y:=lccstr(y);\n R:=lccstr(r);\n\n if type(parse(X),nume ric) and type(parse(Y),numeric) and type(parse(R),numeric) then\n f ct_cdfN2(X,Y,R,result); \n return parse(result);\n else\n retur n 'cdfN2_lcc(x,y,r)';\n end if;\nend proc:" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 39 "This is fast, about 10 milli seconds. " }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 53 "But difficult to t est. For that use a re-formulation:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "'cdfN2(x,y,rho)'= 'Int(pdf N(tau)*cdfN((y-rho*tau)/sqrt(1-rho^2)),tau = -infinity .. x)';" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/-%&cdfN2G6%%\"xG%\"yG%$rhoG-%$IntG6$* &-%%pdfNG6#%$tauG\"\"\"-%%cdfNG6#*&,&F(F2*&F)F2F1F2!\"\"F2-%%sqrtG6#,& F2F2*$)F)\"\"#F2F9F9F2/F1;,$%)infinityGF9F'" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 "For dimension 3 a simil ar reduction is available and attributed to D B Owen." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}} {SECT 1 {PARA 5 "" 0 "" {TEXT -1 13 "Dimension = 3" }}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 89 "The trivariate cas e can be handled through semi-definite integrals using bivariate norma l" }}{PARA 0 "" 0 "" {TEXT -1 29 "distributions, cf Alan Genz, " } {TEXT 258 65 "Numerical Computation of Rectangular Bivariate and Triva riate ..." }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 49 "http://www.s ci.wsu.edu/math/faculty/genz/homepage" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 88 "To see that the used integration routine is highly exact as an exa mple just compare what" }}{PARA 0 "" 0 "" {TEXT -1 88 "happens for int egrating over the usual pdfN (that is coded in the DLL only for a test ): " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 107 "For directly inputting constan ts it not clear that they belong to a correlation matrix, so provide a check:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 566 "check_coefficients:=proc(_r12,_r13,_r23)\nlocal R, d et, ev, remDigits;\n\nremDigits:=Digits;\nDigits:=18;\nR := Matrix(\n \+ [[1,_r12,_r13], \n [1,_r23],\n [1]], \n shape=symmetr ic,\n scan=triangular[upper]);\n\nif (LinearAlgebra:-IsDefinite(R, qu ery=positive_definite)) then\n det:=evalf(LinearAlgebra:-Determinant( R), 2*Digits);\n det:=evalf(det,6);\n print(`valid correlation matri x`, `determinant `=det);\n Digits:=remDigits;\n return 0;\nelse\n p rint(`not a correlation matrix!`);\n Digits:=remDigits;\n return -1; \nend if;\nDigits:=remDigits;\nreturn;\nend proc:" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 35 "Provide an interfa ce to the DLL ..." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 906 "fct3 := define_external(\n 'cdfN3_string_M aple',\n 'C',\n 'x1_str'::string[], 'x2_str'::string[], 'x3_str'::st ring[], \n 'r12_str'::string[], 'r13_str'::string[], 'r23_str'::strin g[],\n 'ret_str'::string[],\n RETURN::integer[4],\n LIB=myDLL):\n\n cdfN3_lcc:=proc(x1,x2,x3,r12,r13,r23)\n#local X::string, Y::string, R: :string, result;\nlocal X1, X2, X3, R12, R13, R23, result, lccstr;\n\n lccstr:=proc(var)\n convert( evalf(parse(convert(var,string)),105), s tring);\nend proc:\n\nresult:=StringTools:-Fill( `0` , Digits_lcc+10); \n#X1:=convert(x1,string); X2:=convert(x2,string); X3:=convert(x3,stri ng);\n#R12:=convert(r12,string); R13:=convert(r13,string); R23:=conver t(r23,string);\nX1:=lccstr(x1): X2:=lccstr(x2): X3:=lccstr(x3):\nR12:= lccstr(r12): R13:=lccstr(r13): R23:=lccstr(r23):\n\nif type(parse(X1), numeric) then\n fct3(X1, X2, X3, R12, R13, R23,result);\n return par se(result);\nelse\n return 'x';\nend if;\nend proc:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 74 "To cross check \+ results one can use a formula, which is attributed to Owen:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 275 " 'cdfN3'='Int(pdfN(xi)*F_Owen(xi), xi= -infinity .. x1)';\nF_Owen:= (xi ,x2,x3, rho12,rho13,rho23) -> \n cdfN2(\n\011\011(x2 - rho12 * xi) / \+ sqrt(1 - rho12*rho12),\n\011\011(x3 - rho13 * xi) / sqrt(1 - rho13*rho 13),\n\011\011(rho23 - rho13 * rho12) / sqrt( (1 - rho12*rho12)*(1 - r ho13*rho13) )\n\011\011);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%&cdfN3G -%$IntG6$*&-%%pdfNG6#%#xiG\"\"\"-%'F_OwenGF+F-/F,;,$%)infinityG!\"\"%# x1G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'F_OwenGf*6(%#xiG%#x2G%#x3G%& rho12G%&rho13G%&rho23G6\"6$%)operatorG%&arrowGF--%&cdfN2G6%*&,&9%\"\" \"*&9'F79$F7!\"\"F7-%%sqrtG6#,&F7F7*&F9F7F9F7F;F;*&,&9&F7*&9(F7F:F7F;F 7-F=6#,&F7F7*&FEF7FEF7F;F;*&,&9)F7*&FEF7F9F7F;F7-F=6#*&,&F7F7*&F9F7F9F 7F;F7,&F7F7*&FEF7FEF7F;F7F;F-F-F-" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 88 "That reduces the problem to dimens ion 2. But it will be very slow for 100 Digits and for" }}{PARA 0 "" 0 "" {TEXT -1 83 "improving speed I would use an external implementati on for cdfN2 - which would make" }}{PARA 0 "" 0 "" {TEXT -1 61 "the te st depend on a repeated use my implementation in dim 2." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 64 "But one can reduce t o dimension one through partial integration:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 267 "'Int(pdfN(x)*F(x ),x=-infinity..x1)' =\n'F(x1)*cdfN(x1)-int(D(F)(x)*cdfN(x),x = -infini ty .. x1)';\n``;\n'diff(cdfN2(a1+b1*x, a2+b2*x,r),x)' =\n 'b2*pdfN(b2 *x+a2)*cdfN((a1-r*a2+(b1-r*b2)*x)/(1-r^2)^(1/2)) +\n b1*pdfN(b1*x+a1 )*cdfN((a2-r*a1+(b2-r*b1)*x)/(1-r^2)^(1/2))';\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%$IntG6$*&-%%pdfNG6#%\"xG\"\"\"-%\"FGF*F,/F+;,$%)infi nityG!\"\"%#x1G,&*&-F.6#F4F,-%%cdfNGF8F,F,-%$intG6$*&--%\"DG6#F.F*F,-F :F*F,F/F3" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%%diffG6$-%&cdfN2G6%,&%#a1G\"\"\"*&%#b1GF,%\"xGF,F,,& %#a2GF,*&%#b2GF,F/F,F,%\"rGF/,&*(F3F,-%%pdfNG6#F0F,-%%cdfNG6#*&,(F+F,* &F4F,F1F,!\"\"*&,&F.F,*&F4F,F3F,F@F,F/F,F,F,,&F,F,*$)F4\"\"#F,F@#F@FGF ,F,*(F.F,-F86#F*F,-F;6#*&,(F1F,*&F4F,F+F,F@*&,&F3F,*&F4F,F.F,F@F,F/F,F ,F,FDFHF,F," }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 45 "where the constant have to be to evaluated as" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 170 " L:=\{\nb1 = -r12/(1-r12^2)^(1/2), a1 = x2/(1-r12^2)^(1/2),\nb2 = -r13/ (1-r13^2)^(1/2), a2 = x3/(1-r13^2)^(1/2),\nr = (r23 - r13 * r12) / sqr t( (1 - r12*r12)*(1 - r13*r13) )\};" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #>%\"LG<'/%#b2G,$*&%$r13G\"\"\",&F+F+*$)F*\"\"#F+!\"\"#F0F/F0/%#a2G*&% #x3GF+F,F1/%#b1G,$*&%$r12GF+,&F+F+*$)F:F/F+F0F1F0/%#a1G*&%#x2GF+F;F1/% \"rG*&,&%$r23GF+*&F*F+F:F+F0F+*&F;F+F,F+F1" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 28 "For short let us write it as" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 270 "theIntegrand:= 'cdfN(xi)' * \n '(b2*pdfN(b2*xi+a2)*c dfN((a1-r*a2+(b1-r*b2)*xi)/(1-r^2)^(1/2)) +\n b1*pdfN(b1*xi+a1)*cdfN ((a2-r*a1+(b2-r*b1)*xi)/(1-r^2)^(1/2)))';\n``;\n'cdfN3' =\n'F(x1)*cdfN (x1)-Int(theIntegrand,x = -infinity .. x1)';\n`F(x)` = 'cdfN2(a1+b1*x, a2+b2*x,r)'; " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%-theIntegrandG*&-% %cdfNG6#%#xiG\"\"\",&*(%#b2GF*-%%pdfNG6#,&*&F-F*F)F*F*%#a2GF*F*-F'6#*& ,(%#a1GF**&%\"rGF*F3F*!\"\"*&,&%#b1GF**&F:F*F-F*F;F*F)F*F*F*,&F*F**$)F :\"\"#F*F;#F;FCF*F**(F>F*-F/6#,&*&F>F*F)F*F*F8F*F*-F'6#*&,(F3F**&F:F*F 8F*F;*&,&F-F**&F:F*F>F*F;F*F)F*F*F*F@FDF*F*F*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%&cdfN3G,&*&-% \"FG6#%#x1G\"\"\"-%%cdfNGF)F+F+-%$IntG6$%-theIntegrandG/%\"xG;,$%)infi nityG!\"\"F*F7" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%%F(x)G-%&cdfN2G6%, &%#a1G\"\"\"*&%#b1GF*%\"xGF*F*,&%#a2GF**&%#b2GF*F-F*F*%\"rG" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 6 "Test 1" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 232 "x1 :=\0111.000000000;\nx2 : =\0110.330000000;\nx3 :=\011-0.500000000;\n\nr12 :=\0110.9;\nr13 := \0110.7;\nr23 :=\0110.8;\n\nst:=time():\nif 0 <= check_coefficients(r1 2,r13,r23) then\n cdfN3_lcc(x1,x2,x3,r12,r13,r23);\n #evalf(%,16);\n end if;\n`seconds`=time()-st;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x1 G$\"+++++5!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x2G$\"*+++I$!\"* " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x3G$!*++++&!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r12G$\"\"*!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r13G$\"\"(!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r23G$\" \")!\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrix G/%-determinant~G$\"#o!\"$" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq%=t *oAhT/6VD&4uoN,'HHGumQVP!G)>!o!\\%[qpXN#*o^$*o/-'p#e8P$[6'H!$0\"" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/%(secondsG$\"%>a!\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 47 "Test it w ith the modification of Owen's formula" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "TheIntegrand:='eval(the Integrand, L)';\nTheIntegral:= 'Int(theIntegrand,xi = -infinity .. x1) ';" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%-TheIntegrandG-%%evalG6$%-theI ntegrandG%\"LG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,TheIntegralG-%$In tG6$%-theIntegrandG/%#xiG;,$%)infinityG!\"\"%#x1G" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 29 "Since the integran d contains " }{XPPEDIT 18 0 "cdfN(xi);" "6#-%%cdfNG6#%#xiG" }{TEXT -1 34 " as factor I will cut off against " }{XPPEDIT 18 0 "-infinity;" "6 #,$%)infinityG!\"\"" }{TEXT -1 32 " (as it will be below 1E-105): " } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "TheIntegral:= 'Int(TheIntegrand,xi = -23.0 .. x1, method = _Gqua d)';" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,TheIntegralG-%$IntG6%%-TheI ntegrandG/%#xiG;$!$I#!\"\"%#x1G/%'methodG%'_GquadG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 59 "Then Maple gives \+ the same result as the DLL implementation:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 182 "'TheIntegral'=ev alf[105](TheIntegral):\n`eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdfN(x1) , L) - TheIntegral =`;\nevalf(eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdf N(x1), L) -TheIntegral,105);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%\\oev al(cdfN2_lcc(a1+b1*x1,~a2+b2*x1,r)*cdfN(x1),~L)~-~TheIntegral~=G" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq)=t*oAhT/6VD&4uoN,'HHGumQVP!G)>!o !\\%[qpXN#*o^$*o/-'p#e8P$[6'H!$0\"" }}}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 6 "Test 2" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 81 "If the determinant is close to zero the correlation coe fficients are choosen from" }}{PARA 0 "" 0 "" {TEXT -1 37 "a valid mat rix R (cf the Genz paper):" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 439 "theta1 :=\0110.01;\ntheta2 :=\0110 .01; #1 - 1e-6;\ntheta3 :=\011-0.02;\n\nR:='R': A:='A': `R`= A*A^t;\n` A`=Matrix('[[1,0,0],\n [cos(theta1*Pi),sin(theta1*Pi),0],\n [cos(the ta2*Pi)*cos(theta3*Pi),cos(theta2*Pi)*sin(theta3*Pi),sin(theta2*Pi)]]' );\nrhs(%):\nA:=evalf(%): #evalf(evalhf(%),15);\nLinearAlgebra:-Multip ly(A, LinearAlgebra:-Transpose(A)):\nR:=evalf(evalhf(%),15); A:='A':\n ``;\nr12:=R[1,2];\nr13:=R[1,3];\nr23:=R[2,3];\n\ncheck_coefficients(r1 2,r13,r23):" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta1G$\"\"\"!\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta2G$\"\"\"!\"#" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%'theta3G$!\"#F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"RG*&%\"AG\"\"\")F&%\"tGF'" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #/%\"AG-%'RTABLEG6%\"*KF\"f9-%'MATRIXG6#7%7%\"\"\"\"\"!F/7%-%$cosG6#*& %'theta1GF.%#PiGF.-%$sinGF3F/7%*&-F26#*&%'theta2GF.F6F.F.-F26#*&%'thet a3GF.F6F.F.*&F;F.-F8F@F.-F8F<%'MatrixG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"RG-%'RTABLEG6%\"*C_*f9-%'MATRIXG6#7%7%$\"\"\"\"\"!$\"0KdOgl] ***!#:$\"01W[iU`(**F37%F1F.$\"0v8([rq]**F37%F4F7F.%'MatrixG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r1 2G$\"0KdOgl]***!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r13G$\"01W[iU `(**!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r23G$\"0v8([rq]**!#:" }} {PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-determin ant~G$\"']M(*!#7" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "x1 := \0111.000000000;\nx2 :=\0110.330000000;\nx3 :=\011-0.400000000;" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x1G$\"+++++5!\"*" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%#x2G$\"*+++I$!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 " 6#>%#x3G$!*++++%!\"*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 78 "Here the computation needs not much more time t o achieve the desired accuracy:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 156 "st:=time():\nif 0 <= check_ coefficients(r12,r13,r23) then\n print(`cdfN3 =`);\n cdfN3_lcc(x1,x2 ,x3,r12,r13,r23);\n #evalf(%,16);\nend if;\n`seconds`=time()-st;" }} {PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-determin ant~G$\"']M(*!#7" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%(cdfN3~=G" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dqB0bsq?Aq&4Q^kU(=J*y5y/5p2hT,+c&)f ^lS?q7l\\Vb`S]D`n*Qe#yXM!$0\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%(se condsG$\"%&z'!\"$" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "TheIn tegrand:='eval(theIntegrand, L)';\nTheIntegral:= 'Int(TheIntegrand,xi \+ = -23.0 .. x1, method = _Gquad)';" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#> %-TheIntegrandG-%%evalG6$%-theIntegrandG%\"LG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,TheIntegralG-%$IntG6%%-TheIntegrandG/%#xiG;$!$I#!\" \"%#x1G/%'methodG%'_GquadG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 182 "'TheIntegral'=evalf[105](TheIntegral):\n`eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdfN(x1), L) - TheIntegral =`;\nevalf(eval(cdfN2_lcc(a1+b 1*x1, a2+b2*x1,r)*cdfN(x1), L) -TheIntegral,105);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%\\oeval(cdfN2_lcc(a1+b1*x1,~a2+b2*x1,r)*cdfN(x1),~L)~- ~TheIntegral~=G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dqC0bsq?Aq&4Q^kU (=J*y5y/5p2hT,+c&)f^lS?q7l\\Vb`S]D`n*Qe#yXM!$0\"" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "Again up to the la st decimal place the results coincide." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 6 "Test 3" }}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 443 "theta1 :=\0110.01;\ntheta2 :=\0111 - 1e-5; #1 - 1e-6;\ntheta3 :=\011-0.02;\n\nR:='R': A:='A': `R`= A*A^t;\n`A`=Matr ix('[[1,0,0],\n [cos(theta1*Pi),sin(theta1*Pi),0],\n [cos(theta2*Pi) *cos(theta3*Pi),cos(theta2*Pi)*sin(theta3*Pi),sin(theta2*Pi)]]'):\nrhs (%):\nA:=evalf(%): #evalf(evalhf(%),15);\nLinearAlgebra:-Multiply(A, L inearAlgebra:-Transpose(A)):\nR:=evalf(evalhf(%),15); A:='A':\n``;\nr1 2:=R[1,2];\nr13:=R[1,3];\nr23:=R[2,3];\n\ncheck_coefficients(r12,r13,r 23):" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta1G$\"\"\"!\"#" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta2G$\"&*****!\"&" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%'theta3G$!\"#F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"RG*&%\"AG\"\"\")F&%\"tGF'" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #>%\"RG-%'RTABLEG6%\"*%Q=f9-%'MATRIXG6#7%7%$\"\"\"\"\"!$\"0KdOgl]***!# :$!0ld$zsE!)**F37%F1F.$!0!z6T'>c&**F37%F4F7F.%'MatrixG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r12G$\"0 KdOgl]***!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r13G$!0ld$zsE!)**!# :" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r23G$!0!z6T'>c&**!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-determinant~G$ \"'mP(*!#=" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "x1 :=\0111.00 0000000;\nx2 :=\0110.330000000;\nx3 :=\011-0.400000000;" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%#x1G$\"+++++5!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x2G$\"*+++I$!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%#x3G$!*++++%!\"*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 156 "st:= time():\nif 0 <= check_coefficients(r12,r13,r23) then\n print(`cdfN3 \+ =`);\n cdfN3_lcc(x1,x2,x3,r12,r13,r23);\n #evalf(%,16);\nend if;\n`s econds`=time()-st;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlat ion~matrixG/%-determinant~G$\"'mP(*!#=" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%(cdfN3~=G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq#*zMw`\"z#fNy 'p]3#Qn.)pTtO&3\")\\*\\$yH)Q78x!3:wQ\"p8mxygHvT-U#4n%!$2\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%(secondsG$\"%'y%!\"$" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 58 "For that case Owen 's function is already difficult to use:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 447 "theF := proc(xi,x2, x3, rho12,rho13,rho23)\nlocal a1,b1,a2,b2,corr, result,r;\n\nb1 := -rh o12/(1-rho12^2)^(1/2); a1 := x2/(1-rho12^2)^(1/2);\nb2 := -rho13/(1-rh o13^2)^(1/2); a2 := x3/(1-rho13^2)^(1/2);\ncorr := (rho23 - rho13 * rh o12) / sqrt( (1 - rho12*rho12)*(1 - rho13*rho13) );\n\nr:=corr;\n#[(a1 -r*a2+(b1-r*b2)*x)/(1-r^2)^(1/2),(a2-r*a1+(b2-r*b1)*x)/(1-r^2)^(1/2)] \+ ;\n\nreturn cdfN2_lcc(a1+b1*xi, a2+b2*xi,corr); #cdfN2(a1+b1*x, a2+b2* x,corr);\nend proc;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%%theFGf*6(%#x iG%#x2G%#x3G%&rho12G%&rho13G%&rho23G6)%#a1G%#b1G%#a2G%#b2G%%corrG%'res ultG%\"rG6\"F5C)>8%,$*&9'\"\"\",&F8$*&9%F< F=FB>8',$*&9(F<,&F8&*&9&F8(*&,&9)F<*&FKF8*FTO-%*cdf N2_lccG6%,&FDF<*&F8F<9$F " 0 "" {MPLTEXT 1 0 85 "theF(xi,x2,x3, r12,r13,r23):\ntst:= unapply(%,xi): # Owen's function for the test data\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "remDigits:=Digits: Digits:=14;\nplot(tst , 0..0.6); #plot(tst, 0.353685..0.353686);\nmid:=0.3536856;\ntst(mid): evalf(%);\nDigits:=remDigits:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'D igitsG\"#9" }}{PARA 13 "" 1 "" {GLPLOT2D 484 277 277 {PLOTDATA 2 "6%-% 'CURVESG6$7io7$\"\"!$\"/gw[5bF%*!#C7$$\"/++]#HyI\"!#:$\"/'*>*4Y*zN!#B7 $$\"/+]([kdW#F/$\"/VD3%zW5\"!#A7$$\"/++v;\\DPF/$\"/<<]QNuPF87$$\"/++D< q8]F/$\"/-PQC[[7!#@7$$\"/+]P\"*y&H'F/$\"/7$$\"/++lN?c7Ffn$\"/KsM]k2hFi n7$$\"/+]U$e6P\"Ffn$\"/\\pKTn)Q\"!#=7$$\"/++&>q0]\"Ffn$\"/$\\?BZ(pLFdo 7$$\"/++DM^I;Ffn$\"/IO()p_!)yFdo7$$\"/++0ytbX_z>Q$Fdp7$$\"/++XDn/?Ffn$\"/t-i*Q9E(Fdp7$$\"/++!y?#>@Ff n$\"/\\=M`8T8!#;7$$\"/+v3wY_AFfn$\"/hp>'GQj#Fdq7$$\"/++IOTqBFfn$\"/9^# )[KCYFdq7$$\"/+v3\">)*\\#Ffn$\"/*4D&>2k#)Fdq7$$\"/+DEP/BEFfn$\"/dTUd;' Q\"F/7$$\"/+](o:;v#Ffn$\"/B[e`B#H#F/7$$\"/+v$)[opGFfn$\"/k.[&[8_$F/7$$ \"/]7t;OLHFfn$\"/GwfDM#Q%F/7$$\"/+]i%Qq*HFfn$\"/0@'HecS&F/7$$\"/]i]2=j IFfn$\"/#=_*y**fmF/7$$\"/+vQIKHJFfn$\"/YB:ViG\")F/7$$\"/++&4+p=$Ffn$\" /#evs(R&f*F/7$$\"/+D^rZWKFfn$\"/cE)34[7\"Ffn7$$\"/DJ]&pbF$Ffn$\"/;4j.3 A7Ffn7$$\"/]P\\>m1LFfn$\"/aG:,6D8Ffn7$$\"/vV[VvPLFfn$\"/%*f(4tRV\"Ffn7 $$\"/+]Zn%)oLFfn$\"/ZzO!>([:Ffn7$$\"/]7Q#o4S$Ffn$\"/&yG=lMn\"Ffn7$$\"/ +vG(*3LMFfn$\"/V9\\z]/=Ffn7$$\"/]P>7@lMFfn$\"/QNxo!=%>Ffn7$$\"/++5FL( \\$Ffn$\"/XHMiG&3#Ffn7$$\"/v$f,XI^$Ffn$\"/flg!)od@Ffn7$$\"/](=Kd(GNFfn $\"/#Ffn7$$\"/+vL>=gNFfn$\"/< -REJ`?Ffn7$$\"/voRU*ed$Ffn$\"/7JjJ19>Ffn7$$\"/]iXlg\"f$Ffn$\"/zjKyv!y \"Ffn7$$\"/Dc^)=tg$Ffn$\"/GC[9Y`;Ffn7$$\"/+]d6.BOFfn$\"/%G'ol?K:Ffn7$$ \"/iS@OBQOFfn$\"/cgq&Q1U\"Ffn7$$\"/DJ&3OMl$Ffn$\"/#**>4.ZJ\"Ffn7$$\"/) =#\\&Q'oOFfn$\"/bKU3N97Ffn7$$\"/]785%Qo$Ffn$\"/T5@(4&>6Ffn7$$\"/v$4%fC 9PFfn$\"/#R))>)of%*F/7$$\"/+vo3lWPFfn$\"/+N(3\\6$zF/7$$\"/v=d.TyPFfn$ \"/!\\l$p`hkF/7$$\"/]iX)p@\"QFfn$\"/s[*eWF@&F/7$$\"/D1M$Hf%QFfn$\"/;Y \")fvjTF/7$$\"/+]A))ozQFfn$\"/+#*=Mq#H$F/7$$\"/+DEwNSRFfn$\"/Q'oI?`5#F /7$$\"/++Ik-,SFfn$\"/0I$[BBI\"F/7$$\"/+]FL!e1%Ffn$\"/=#\\+&G:vFdq7$$\" /++D-eITFfn$\"/y%on=F<%Fdq7$$\"/](=sx#*=%Ffn$\"/j'fu+uO#Fdq7$$\"/+v=_( zC%Ffn$\"/Z6sUc+8Fdq7$$\"/](o3Z@J%Ffn$\"/m-olT5lFdp7$$\"/++b*=jP%Ffn$ \"/lp*e*HMJFdp7$$\"/+v3/3(\\%Ffn$\"/%oXQm#>rFdo7$$\"/+vB4JBYFfn$\"/^^@ \"4;I\"Fdo7$$\"/++DVsYZFfn$\"/ZD`I2H@Fin7$$\"/+v=n#f([Ffn$\"/)R\"zpjGF FN7$$\"/++!)RO+]Ffn$\"/hM`-3KKFC7$$\"/+]_!>w7&Ffn$\"/wU1!)H;JF87$$\"/+ v)Q?QD&Ffn$\"/Y0$[dFfn$\"/8*fIK ?Z$!#G7$$\"/+D6EjpeFfn$\"/X*[SYE_\"!#H7$$\"\"'!\"\"$\"/)Ryd*e)[%!#J-%' COLOURG6&%$RGBG$\"#5F^dl$F(F(Fhdl-%+AXESLABELSG6$Q!6\"F\\el-%%VIEWG6$; FhdlF\\dl%(DEFAULTG" 1 2 0 1 10 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 0 "Curve 1" }}}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$midG$\"( co`$!\"(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"/4N(>0)pA!#9" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "Ev en if the function is actually analytic it has a numerical peak." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "To avoid \+ problem in applying the modification one re-arranges the integrand to \+ be positive" }}{PARA 0 "" 0 "" {TEXT -1 91 "(zero correspond to a peak through D(F) ), i.e. rho12, rho13 have the same sign, which I do" }} {PARA 0 "" 0 "" {TEXT -1 34 "by hand here (and sketch the way):" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "`sigma^(-1): 1 -> 3, 2 -> 1, 3 -> 2`;\n`sigma: 1 -> 2, 2 -> 3, 3 -> 1`;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%Csigma^(-1):~1~->~3,~ 2~->~1,~3~->~2G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%Csigma:~~~~~~1~->~ 2,~2~->~3,~3~->~1G" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "`sigm a gives x_new` = 'x3, x1, x3';" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%2s igma~gives~x_newG6%%#x3G%#x1GF&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "`tau` = `sigma^(-1)`;\n'r_new[i,j]' = r_old[tau(i),tau(j)];\n' r_new[i,j]' = r_old[tau(j),tau(i)];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #/%$tauG%+sigma^(-1)G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%&r_newG6$% \"iG%\"jG&%&r_oldG6$-%$tauG6#F'-F-6#F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%&r_newG6$%\"iG%\"jG&%&r_oldG6$-%$tauG6#F(-F-6#F'" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "'r_new[1,2]' = r[1,3]; #rho[3,1];\n 'r_new[1,3]' = r[2,3]; #rho[3,2];\n'r_new[2,3]' = r[1,2];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%&r_newG6$\"\"\"\"\"#&%\"rG6$F'\"\"$" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/&%&r_newG6$\"\"\"\"\"$&%\"rG6$\"\"#F( " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%&r_newG6$\"\"#\"\"$&%\"rG6$\"\" \"F'" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "# new values\nr12 : = -0.998026727935765;\nr13 := -0.995561964111790;\nr23 := 0.9995065603 65732;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r12G$!0ld$zsE!)**!#:" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r13G$!0!z6T'>c&**!#:" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%$r23G$\"0KdOgl]***!#:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "# new values\nx1 :=\011-0.400000000;\nx2 :=\0111 .000000000;\nx3 :=\0110.330000000;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%#x1G$!*++++%!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x2G$\"+++++5! \"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x3G$\"*+++I$!\"*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 40 "For the D LL that change does not matter:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 180 "st:=time():\nif 0 <= check_ coefficients(r12,r13,r23) then\n print(`cdfN3 =`);\n cdfN3_lcc(x1,x2 ,x3,r12,r13,r23);\n #evalf(%,16);\nend if;\nnew_result_from_DLL:=%:\n `seconds`=time()-st;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correl ation~matrixG/%-determinant~G$\"'mP(*!#=" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%(cdfN3~=G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq#*zMw`\"z#fNy 'p]3#Qn.)pTtO&3\")\\*\\$yH)Q78x!3:wQ\"p8mxygHvT-U#4n%!$2\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%(secondsG$\"%'y%!\"$" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 47 "but Owen's functio n now looks better to handle:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "theF(xi,x2,x3, r12,r13,r23): \ntst:=unapply(%,xi):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "re mDigits:=Digits: Digits:=14;\nplot(tst, -2.. x1);\nDigits:=remDigits: " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'DigitsG\"#9" }}{PARA 13 "" 1 " " {GLPLOT2D 485 255 255 {PLOTDATA 2 "6%-%'CURVESG6$7jn7$!\"#$\"/x/=`h \")\\!#$)7$$!/LL`X7l>!#8$\"/X0llA+K!#!)7$$!/nm8'zZ$>F/$\"/x9>/%)[z!#y7 $$!/LLbNl+>F/$\"/Q;Zq/sM!#v7$$!/LL(G,j'=F/$\"/hJJ#e[Q\"!#s7$$!/nm*G7@$ =F/$\"/!e!R(e6r%!#q7$$!/LL7ZT+=F/$\"/97w/IG$*!#o7$$!/++2Pfn!#l7$$!/LLw,lLl:F/$\"/G0hAYM*R\"F/$\"/8%pe;rz' !#V7$$!/nm)p*)yO\"F/$\"/9\"3%3>TG!#T7$$!/++r:QL8F/$\"/g*\\!GA/:!#R7$$! /++t;_+8F/$\"//peP-Fe!#Q7$$!/nm;eBm7F/$\"/$oX=+%GB!#O7$$!/nm(p]ZB\"F/$ \"/52v6ePh!#N7$$!/LLV(*y+7F/$\"/0DJ;>]=!#L7$$!/LLcQ^l6F/$\"/26+[@fb!#K 7$$!/++$41[8\"F/$\"/GSJgR6'*!#J7$$!/LLn3k,6F/$\"/]ite$p&=!#H7$$!/++WzP n5F/$\"/T#R#zU$[$!#G7$$!/++e$eQ.\"F/$\"/t7;0=5a!#F7$$!/++NkU,5F/$\"/[= rIRLo!#E7$$!.++9jTl*F/$\"/SGb7^r**!#D7$$!.nm='fI$*F/$\"/FrsSC>)*!#C7$$ !.++S>^)*)F/$\"/uoJAHZ**!#B7$$!.LL3m?n)F/$\"/fpbLCXs!#A7$$!.++7;)H$)F/ $\"/r,w(y/i&!#@7$$!.nmd&y2!)F/$\"/![#efSRM!#?7$$!.++(37$$!.nmYo?M(F/$\"/AODV&)>5!#=7$$!.++vGv*pF/$\"/6.6EV*)[Fjy7$$!.LL0 'plmF/$\"/s>aYFj>!#<7$$!.LL>\\jK'F/$\"/r/1t@8sFez7$$!/++X23eh!#9$\"/RZ >5*\\J\"!#;7$$!.nmH7)*)fF/$\"/$*p=%zwK#Fa[l7$$!/NLo5>NeF^[l$\"/z]Q3aLQ Fa[l7$$!.++%)p0o&F/$\"/m0#ec/;'Fa[l7$$!/NLjkN.bF^[l$\"/!e\"*e-(H5!#:7$ $!.nm3VhK&F/$\"/iNg3Qn;Ff\\l7$$!/++I8kn^F^[l$\"/U&3UY')\\#Ff\\l7$$!.LL dR\"4]F/$\"/uK:!*\\_OFf\\l7$$!/lmTO9S[F^[l$\"/]?*R*>H`Ff\\l7$$!.++rZ6n %F/$\"/H:f3#\\c(Ff\\l7$$!/++vBF!f%F^[l$\"/&zC/w.'))Ff\\l7$$!/++SqR4XF^ [l$\"/Wm?mUJ5F^[l7$$!/++0<_GWF^[l$\"/u8yuS$>\"F^[l7$$!.++PYwM%F/$\"/x& zCeDP\"F^[l7$$!/+]xZtgUF^[l$\"/rw)pWYe\"F^[l7$$!/++&=BQ<%F^[l$\"/hgKt: <=F^[l7$$!/+v)Qn.8%F^[l$\"/m!*)RK5%>F^[l7$$!/+]#f6p3%F^[l$\"/gH\\\">*p ?F^[l7$$!/+D'zbM/%F^[l$\"/Ul(GPP?#F^[l7$$!*++++%!\"*$\"/:JHAQUBF^[l-%' COLOURG6&%$RGBG$\"#5!\"\"$\"\"!F[blFjal-%+AXESLABELSG6$Q!6\"F_bl-%%VIE WG6$;$F(F[blF^al%(DEFAULTG" 1 2 0 1 10 0 2 9 1 4 2 1.000000 46.000000 45.000000 0 0 "Curve 1" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "TheIntegrand:='eval(th eIntegrand, L)';\nTheIntegral:= 'Int(TheIntegrand,xi = -23.0 .. x1, me thod = _Gquad)';" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%-TheIntegrandG-% %evalG6$%-theIntegrandG%\"LG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,The IntegralG-%$IntG6%%-TheIntegrandG/%#xiG;$!$I#!\"\"%#x1G/%'methodG%'_Gq uadG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 182 "'TheIntegral'=eval f[105](TheIntegral):\n`eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdfN(x1), \+ L) - TheIntegral =`;\nevalf(eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdfN( x1), L) -TheIntegral,105);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%\\oeval (cdfN2_lcc(a1+b1*x1,~a2+b2*x1,r)*cdfN(x1),~L)~-~TheIntegral~=G" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#$\"cqqZjP:z#fNy'p]3#Qn.)pTtO&3\")\\*\\ $yH)Q78x!3:wQ\"p8mxygHvT-U#4n%!$1\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "and up to the last places that \+ approach gives the same" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 6 "Test 4" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 443 "theta1 :=\0111e-3;\ntheta2 :=\0111 - 1e-6; #1 - 1e-6 ;\ntheta3 :=\011-1e-3;\n\nR:='R': A:='A': `R`= A*A^t;\n`A`=Matrix('[[1 ,0,0],\n [cos(theta1*Pi),sin(theta1*Pi),0],\n [cos(theta2*Pi)*cos(th eta3*Pi),cos(theta2*Pi)*sin(theta3*Pi),sin(theta2*Pi)]]'):\nrhs(%):\nA :=evalf(%): #evalf(evalhf(%),15);\nLinearAlgebra:-Multiply(A, LinearAl gebra:-Transpose(A)):\nR:=evalf(evalhf(%),15); A:='A':\n``;\nr12:=R[1, 2];\nr13:=R[1,3];\nr23:=R[2,3];\n\ncheck_coefficients(r12,r13,r23):" } }{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta1G$\"\"\"!\"$" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta2G$\"'******!\"'" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'theta3G$!\"\"!\"$" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #/%\"RG*&%\"AG\"\"\")F&%\"tGF'" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\" RG-%'RTABLEG6%\"*KF\"f9-%'MATRIXG6#7%7%$\"\"\"\"\"!$\"0e=?l]*****!#:$! 0Bp>l]*****F37%F1F.$!0-7&3E!)****F37%F4F7F.%'MatrixG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r12G$\"0e= ?l]*****!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r13G$!0Bp>l]*****!#: " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r23G$!0-7&3E!)****!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-determinant~G$ \"'BU(*!#A" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "x1 :=\0111.00 0000000;\nx2 :=\0110.330000000;\nx3 :=\011-0.400000000;" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%#x1G$\"+++++5!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x2G$\"*+++I$!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%#x3G$!*++++%!\"*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 73 "For that case Owen's function can not be used, it would evaluate to zero:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "theF(xi,x2,x3, r12,r13,r23):\ntst:= unapply(%,xi):\ntst(xi): op(%): evalf(%,8): cdfN2_lcc(%);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%*cdfN2_lccG6%,&$\")WU]5!\"&\"\"\"*&$\")%)3$=$ F)F*%#xiGF*!\"\",&$\")5Ct7F)F/*&$\")o3$=$F)F*F.F*F*$\")]******!\")" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 86 "G enerally the DLL uses sorting on the correlation coefficients and take s care that the" }}{PARA 0 "" 0 "" {TEXT -1 60 "integral is provided w ith a 'positive integrand' as follows:" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 808 "arrangeXR:=proc(x,r)\n local R,X,swapR,swapX;\n\nX := [x[1],x[2],x[3]];\nR := [r[1],r[2],r[3] ];\n\nif R[2] < R[1] then # swapR them\n swapR:= R[2];\n R[2]:=R[1]; \n R[1]:=swapR;\n\n swapX:= X[2];\n X[2]:=X[3];\n X[3]:=swapX;\nen d if;\n#print(R);\n\nif R[3] < R[2] then # swapR them\n swapR:= R[3]; \n R[3]:=R[2];\n R[2]:=swapR;\n\n swapX:=X[1];\n X[1]:=X[2];\n X[ 2]:=swapX;\nend if;\n#print(R);\n\nif R[2] < R[1] then # swapR them\n \+ swapR:=R[2];\n R[2]:=R[1];\n R[1]:=swapR;\n\n swapX:= X[2];\n X[2 ]:=X[3];\n X[3]:=swapX;\nend if;\n#print([X,R]);\n\nif ( R[3] < 0 or \+ 0 < R[1] or R[1] = 0 or R[2] = 0) then\nelse\n if 0 < R[2] then\nprin t(\"we do it \", [X,R]);\n swapR:=R[1];\n R[1]:=R[2];\n R[2]: =R[3];\n R[3]:=swapR;\n\n swapX:=X[1];\n X[1]:=X[3];\n X[3 ]:=X[2];\n X[2]:=swapX;\n end if;\nend if;\n\nreturn [X,R];\nend p roc;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%*arrangeXRGf*6$%\"xG%\"rG6&% \"RG%\"XG%&swapRG%&swapXG6\"F.C)>8%7%&9$6#\"\"\"&F46#\"\"#&F46#\"\"$>8 $7%&9%F5&FAF8&FAF;@$2&F>F8&F>F5C(>8&FF>FFFG>FGFJ>8'&F1F8>FO&F1F;>FQFN@ $2&F>F;FFC(>FJFU>FUFF>FFFJ>FN&F1F5>FenFO>FOFN@$FEC(>FJFF>FFFG>FGFJ>FNF O>FOFQ>FQFN@%5552FU\"\"!2FeoFG/FGFeo/FFFeoF.@$2FeoFFC+-%&printG6$Q*we~ do~it~F.7$F1F>>FJFG>FGFF>FFFU>FUFJ>FNFen>FenFQ>FQFO>FOFNOF`pF.F.F." }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 30 "I simply use that for the test" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 166 "X:=[x1,x2,x3]: R:=[r12,r13, r23]:\nLL:=arrangeXR(X,R):\n\n# new values\nr12 := LL[2][1];\nr13 := L L[2][2];\nr23 := LL[2][3];\n\nx1 :=\011LL[1][1];\nx2 :=\011LL[1][2];\n x3 :=\011LL[1][3];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r12G$!0Bp>l]* ****!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r13G$!0-7&3E!)****!#:" } }{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r23G$\"0e=?l]*****!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x1G$!*++++%!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x2G$\"+++++5!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%#x3G$\"*+++I$!\"*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 180 "st: =time():\nif 0 <= check_coefficients(r12,r13,r23) then\n print(`cdfN3 =`);\n cdfN3_lcc(x1,x2,x3,r12,r13,r23);\n #evalf(%,16);\nend if;\nn ew_result_from_DLL:=%:\n`seconds`=time()-st;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-determinant~G$\"'BU(*!#A " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%(cdfN3~=G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq$3'Qz3>nr1WS(ewn]EUR\\x7)p8e;0E1'*R1mU%\\s%48iz%)) y'>)fFD>i@#)!$P\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%(secondsG$\"%wV !\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "where Owen's function becomes somewhat nicer at least:" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 164 "theF(xi,x2,x3, r12,r13,r23):\ntst:=unapply(%,xi):\ntst(xi): op( %): evalf(%,8): cdfN2_lcc(%);\n\nremDigits:=Digits: Digits:=14:\nplot( tst, -0.41.. x1);\nDigits:=remDigits:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%*cdfN2_lccG6%,&$\")D5$=$!\"&\"\"\"*&$\")o3$=$F)F*%#xiGF*F*,&$\")q 9__!\"'F**&$\")$G:f\"F)F*F.F*F*$\")()******!\")" }}{PARA 13 "" 1 "" {GLPLOT2D 485 253 253 {PLOTDATA 2 "6%-%'CURVESG6$7en7$$!#T!\"#$\"/0g`' RP)>!#]7$$!/Le%G?y4%!#9$\"/Dtvjt\"4$F-7$$!/U&esBf4%F1$\"/w8\")GDWXF-7$ $!/$3s%3z$4%F1$\"/lGK`$***pF-7$$!/$e/$Qk\"4%F1$\"/0'RS2,3\"!#\\7$$!/U5 =q]*3%F1$\"/C`h\\Gh;FC7$$!/3_>f_(3%F1$\"/r;cmrtCFC7$$!/vo1YZ&3%F1$\"/@ #G'\\\">t$FC7$$!/3-OJN$3%F1$\"/wh;lW.dFC7$$!/v$*o%Q73%F1$\"/(*H#3g[p)F C7$$!/n\"RFj!zSF1$\"/ny4\\/S8!#[7$$!/$e4OZr2%F1$\"/m#f5X&f>F\\o7$$!/+v '\\!*\\2%F1$\"/4sNdQ-IF\\o7$$!/+DwZ#G2%F1$\"/\\n]:&Gg%F\\o7$$!/+D.xtqS F1$\"/fYG*o+%pF\\o7$$!/3-TC%)oSF1$\"/s-1dp15!#Z7$$!/n\"4z)emSF1$\"/(Ga \\I[c\"Ffp7$$!/nm`'zY1%F1$\"/o%4.B:F#Ffp7$$!/v=t)eC1%F1$\"/rQ4J9+NFfp7 $$!/n;1J\\gSF1$\"/'HVAin7&Ffp7$$!/v=[jLeSF1$\"/%3-)zH%y(Ffp7$$!/Dc/EGc SF1$\"/(y+jPt:\"!#Y7$$!/<_SF1$\"/R)[M t;b#Fer7$$!/$e*e$\\+0%F1$\"/6Cl\"p<$QFer7$$!/3-;Y%y/%F1$\"/^ltoLQeFer7 $$!/D\"3QDf/%F1$\"/L#Qtp^T)Fer7$$!/$3Ub_Q/%F1$\"/%>jRWwC\"!#X7$$!/+]@6 rTSF1$\"/&4WP;>(=Fdt7$$!/]PZhhRSF1$\"/#[/gJ3y#Fdt7$$!/v=_\"*ePSF1$\"/n !*)4vS2%Fdt7$$!/]i>&Q`.%F1$\"/\\MJy2=iFdt7$$!/n;EiJLSF1$\"/**4&f'*>3*F dt7$$!/+D'*p:JSF1$\"/EMw4Uf8!#W7$$!/3-8/?HSF1$\"/*3!*Q^s&>Fcv7$$!/+v]8 1FSF1$\"/BZ8x@7HFcv7$$!/U&)f'[]-%F1$\"/cyA)Q\"GUFcv7$$!/v$z\"[%H-%F1$ \"/&3q0:kB'Fcv7$$!/n\"z#z)3-%F1$\"/\\DvUL4\"*Fcv7$$!/voaXt=SF1$\"/c\\z \"))GN\"!#V7$$!/LL+1m;SF1$\"/RURf)z(>Fbx7$$!/$eCoRX,%F1$\"/\\-E'GO\"HF bx7$$!/U&oKOC,%F1$\"/\"yu#yLtUFbx7$$!/+]hN]5SF1$\"/2Tui$*pgFbx7$$!//SF1$\"/Ut\"*45'*=F az7$$!/QM-#*o.SF1$\"/&p/w\"Qw?Faz7$$!/+]JP=.SF1$\"/+sYvltAFaz7$$!/ilg# yE+%F1$\"/Fz&\\;&*[#Faz7$$!/D\")*ys@+%F1$\"/L.OQpDFFaz7$$!/%fBfH;+%F1$ \"/]*f&zJ/IFaz7$$!/i!\\R'3,SF1$\"/iFDy<6LFaz7$$!/JX(>V0+%F1$\"/.Q&f4\" \\OFaz7$$!*++++%!\"*$\"/u`? " 0 "" {MPLTEXT 1 0 342 "TheIntegrand:='eval(theIntegrand, L)';\nTh eIntegral:= 'Int(TheIntegrand,xi = -23.0 .. x1, method = _Gquad)';\n`` ;\n'TheIntegral'=evalf[105](TheIntegral):\n`eval(cdfN2_lcc(a1+b1*x1, a 2+b2*x1,r)*cdfN(x1), L) - TheIntegral =`;\nevalf(eval(cdfN2_lcc(a1+b1* x1, a2+b2*x1,r)*cdfN(x1), L) -TheIntegral,105);\n`computational error` = % -new_result_from_DLL;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%-TheIn tegrandG-%%evalG6$%-theIntegrandG%\"LG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,TheIntegralG-%$IntG6%%-TheIntegrandG/%#xiG;$!$I#!\"\"%#x1G/%' methodG%'_GquadG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%!G" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#%\\oeval(cdfN2_lcc(a1+b1*x1,~a2+b2*x1,r)*cdfN(x1 ),~L)~-~TheIntegral~=G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"`q?m)e!R( [A%yROBI$*>%R\\x7)p8e;0E1'*R1mU%\\s%48iz%))y'>)fFD>i@#)!$L\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%4computational~errorG$!B$3'=8?8G%==?w%Hawr #!$P\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 87 "The computational error is less than 1E-104 which is the \+ exactness of qfloat (it is not" }}{PARA 0 "" 0 "" {TEXT -1 31 "in digi ts, but decimal places)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 6 "Test 5" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 86 "The correlation coefficients need \+ not be close to 1 to give an ugly matrix neither the" }}{PARA 0 "" 0 " " {TEXT -1 60 "determinant need to be extremly small to cause difficul ties:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 93 "r12:=1/2;\nr13:=1/2;\nr23:=-1/2 + 10^(-3); r23 := eva lf(r23):\n\ncheck_coefficients(r12,r13,r23):" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r12G#\"\"\"\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%$r13G#\"\"\"\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$r23G#!$*\\\" %+5" }}{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-d eterminant~G$\"'+*\\\"!\")" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "x1 :=\0111.000000000;\nx2 :=\0110.330000000;\nx3 :=\011-0.40000000 0;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x1G$\"+++++5!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x2G$\"*+++I$!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x3G$!*++++%!\"*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 43 "The DLL asserts an result in moder ate time:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 180 "st:=time():\nif 0 <= check_coefficients(r12,r13,r2 3) then\n print(`cdfN3 =`);\n cdfN3_lcc(x1,x2,x3,r12,r13,r23);\n #e valf(%,16);\nend if;\nnew_result_from_DLL:=%:\n`seconds`=time()-st;" } }{PARA 11 "" 1 "" {XPPMATH 20 "6$%9valid~correlation~matrixG/%-determi nant~G$\"'+*\\\"!\")" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%(cdfN3~=G" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq?^F.+(Hq)GI(*eCj'=2%Q?3#=kNTO-VOV @5.GZ?=/Pjg8Lqds8.SR'HU\"!$0\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%(s econdsG$\"%;%)!\"$" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "TheIn tegrand:='eval(theIntegrand, L)';\nTheIntegral:= 'Int(evalf(TheIntegra nd),xi = -23.0 .. x1)';" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%-TheInteg randG-%%evalG6$%-theIntegrandG%\"LG" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 #>%,TheIntegralG-%$IntG6$-%&evalfG6#%-TheIntegrandG/%#xiG;$!$I#!\"\"%# x1G" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 98 "And looking at the integrand Maple should be able to compute th at (otherwise split it into pieces," }}{PARA 0 "" 0 "" {TEXT -1 100 "f or r23 towards -1/2 it becomes quite ugly ... I lost patients with Ma ple to integrate then and did" }}{PARA 0 "" 0 "" {TEXT -1 34 "not work out an example for that):" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 194 "remDigits:=Digits: Digits:=14:\npl ot(TheIntegrand,xi = -3.0 .. x1);\nDI:=evalf(diff(TheIntegrand,xi)): e xpand(%): Tryhard(%): combine(%,exp): evalf(%,14):\nfsolve(%=0, xi=-1. .0);\nDigits:=remDigits:" }}{PARA 13 "" 1 "" {GLPLOT2D 485 243 243 {PLOTDATA 2 "6%-%'CURVESG6$7do7$$!\"$\"\"!$!3%y`3S3P7s\"!#@7$$!3cLLL$Q 6G\"H!#<$!3eScs\"ROv3:%p'**>F1$!3`OJLlTiRdFM7$$!3.+++]5*H\">F1$!3/$G8CSDEN (FM7$$!3,+++I\"3&H=F1$!3]L_WX/pd#*FM7$$!3OLL$3k(p` 7$$!3smmmO;bj;F1$!3OVx8F#[*G9Fep7$$!3!ommm9'=(e\"F1$!3G4Iy/-VE9F1$!3p!>KqwvPb#Fep7$$! 3\"*****\\FRXL8F1$!33Xb#pp1i3$Fep7$$!3)*****\\#=/8D\"F1$!3jN_>5TlnOFep 7$$!3immmT&*el6F1$!3rhz'4+*fcVFep7$$!3pmm;Wn(o3\"F1$!3:%)RH$e$em]Fep7$ $!3PLLLeV(>+\"F1$!32R**y@$\\x\"fFep7$$!3iOLL3k%y8*!#=$!3yTF$=\\7m*oFep 7$$!31-++DB:q$)Ffs$!3n;<)exal#yFep7$$!3XNLL$o@5a(Ffs$!3pDRz;3s4*)Fep7$ $!3T,+++'[Wo'Ffs$!3nVS@qC)4,\"Ffs7$$!3;/++]*ek%eFfs$!32&[\"=z!ec8\"Ffs 7$$!3:.++v3mN]Ffs$!3%)>Mx9$oAE\"Ffs7$$!3b.++]ySNTFfs$!3`2=rF<\\39Ffs7$ $!3[pmmm/\\ELFfs$!3u#fc8N]Oa\"Ffs7$$!3z++++&)ziCFfs$!3e]O\"Q7r-p\"Ffs7 $$!3)zmmT&=[r?Ffs$!3;3m:Ffs$!3Xe.@r\"e$*y\"Ffs7$$!3a,+D\")*eiY\"Ffs$!378IeAep zdIVg*)*4O\"Ffs7$$!3q1+++ISX#)Fep$!3^=Zt /KS+7Ffs7$$!3XWL3_v0RsFep$!3[*oPHWy_.\"Ffs7$$!3?#omT57FB'Fep$!3Kg$[3Gk ,l)Fep7$$!3&*>+DcmOE_Fep$!3KUOqjmG&)pFep7$$!3qdLL37-?UFep$!3&f2_fwL-W& Fep7$$!3X&p;/wvO@$Fep$!36AZBpzrySFep7$$!3?L+]7.L2AFep$!3w9\\@U?2RHFep7 $$!3&4P$ek[)4?\"Fep$!3Wq8hYwtK?Fep7$$!3)p3nm;%RY>FM$!3B]a&y')*zZ8Fep7$ $\"3Ri*****\\-#4>Fep$!33/1L`F/g\\FM7$$\"3[Lmm;W/8SFep$!3S5s$z#4S+:FM7$ $\"3.p***\\Pl\\1&Fep$!3'\\LaChCRl(F-7$$\"3c/LLLj)o6'Fep$!3,Q\\*=Mh1r$F -7$$\"36Smm\"H2)orFep$!3A\\*[x;r)3lmmm5:xuFfsFi_l7$$\"3&4++]sI@K)FfsFi _l7$$\"3)3++]2%)38*FfsFi_l7$$\"\"\"F*Fi_l-%'COLOURG6&%$RGBG$\"#5!\"\"$ F*F*Fjal-%+AXESLABELSG6$Q#xi6\"Q!F_bl-%%VIEWG6$;$!#IFial$\"+++++5!\"*% (DEFAULTG" 1 2 0 1 10 0 2 9 1 4 2 1.000000 45.000000 46.000000 0 0 "Cu rve 1" }}}{PARA 11 "" 1 "" {XPPMATH 20 "6#$!/4:>$pTl\"!#9" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "Again hav e patients to let Maple compute it up to 105 Digits and then compare" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "evalf(TheIntegrand): normal(%): evalf(%):\nTheIntegral:= 'Int (%,xi = -23.0 .. x1)':\nTheI_105:=evalf(%,105);\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%)TheI_105G$!dq@^F.+(Hq)GI(*eCj'=2%Q?3#=kNTO-VOV@5.GZ? =/Pjg8Lqds8.SR'HU\"!$0\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 188 "`eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdfN(x1), L) - TheIntegral \+ =`;\nevalf(eval(cdfN2_lcc(a1+b1*x1, a2+b2*x1,r)*cdfN(x1), L) -TheI_105 ,105);\n`computational error` = % -new_result_from_DLL;" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#%\\oeval(cdfN2_lcc(a1+b1*x1,~a2+b2*x1,r)*cdfN(x1 ),~L)~-~TheIntegral~=G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"dq@^F.+(H q)GI(*eCj'=2%Q?3#=kNTO-VOV@5.GZ?=/Pjg8Lqds8.SR'HU\"!$0\"" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#/%4computational~errorG$\"\"\"!$0\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 "which sho ws that within the exactness of qfloat the result is correct." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "1 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }{RTABLE_HANDLES 145912732 }{RTABLE M7R0 I6RTABLE_SAVE/145912732X,%'sfloatG6"6"[gl!"%!!!#*"$"$$"""""!$"0e=?l]*****!#:$!0 Bp>l]*****F,F*F'$!0-7&3E!)****F,F-F/F'6" }