*! Date : 17 Aug 2021 *! Version : 1.6 *! Authors : Michael J Grayling & Adrian P Mander /* 16/04/15 v1.0 Basic version complete. 01/06/15 v1.1 Extended to use variable re-ordering. 20/09/15 v1.2 Removed dependency on integrate() being pre-installed and tailored integration to our particular requirements. 21/09/15 v1.3 Made small changes for improved efficiency. 08/09/17 v1.4 Changed so df = 0 or df = . corresponds to mvnormal. Added support for non-integer df. 03/01/19 v1.5 Minor changes for speed. Converted delta and sigma to be optional with internal defaults. 17/08/21 v1.6 Fixed a small printing error. */ program define mvt, rclass version 15.0 syntax , LOWer(numlist miss) UPPer(numlist miss) [DELta(numlist) /// Sigma(string) df(real 1) INTegrator(string) SHIfts(integer 12) /// SAMples(integer 1000) ALPha(real 3)] ///// Check input variables //////////////////////////////////////////////////// if ("`integrator'" == "") { local integrator "pmvnormal" } local lenlower:list sizeof lower local lenupper:list sizeof upper local lendelta:list sizeof delta if (`lenlower' > 100) { di "{error}Only multivariate t distributions of dimension up to 100 are supported." exit(198) } if (`lenlower' != `lenupper') { di "{error}Vector of lower limits (lower) and vector of upper limits (upper) must be of equal length." exit(198) } forvalues i = 1/`lenlower' { local loweri:word `i' of `lower' local upperi:word `i' of `upper' if ((`loweri' != .) & (`upperi' != .)) { if (`loweri' >= `upperi') { di "{error}Each lower integration limit (in lower) must be strictly less than the corresponding upper limit (in upper)." exit(198) } } } if ((`lendelta' != 0) & (`lenlower' != `lendelta')) { di "{error}Vector of lower limits (lower) and vector of non-centrality parameters (delta) must be of equal length." exit(198) } if ("`sigma'" != "") { if (colsof(`sigma') != rowsof(`sigma')) { di "{error}Scale matrix Sigma (sigma) must be square." exit(198) } if (`lenlower' != colsof(`sigma')) { di "{error}Vector of lower limits (lower) must be the same length as the dimension of scale matrix Sigma (sigma)." exit(198) } cap mat chol = cholesky(`sigma') if (_rc > 0) { di "{error}Scale matrix Sigma (sigma) must be symmetric positive-definite." exit(198) } mat sigma = (`sigma') } else { mat sigma = I(`lenlower') } if (`df' < 0) { di "{error}Degrees of freedom (df) must be greater than or equal to zero." exit(198) } if (("`integrator'" != "mvnormal") & ("`integrator'" != "pmvnormal")) { di "{error}Choice of integrator (integrator) must be mvnormal or pmvnormal." exit(198) } if (`shifts' < 1) { di "{error}Number of shifts of the Quasi-Monte-Carlo integration algorithm to use (shifts) must be a strictly positive integer." exit(198) } if (`samples' < 1) { di "{error}Number of samples to use in each shift of the Quasi-Monte-Carlo integration algorithm (samples) must be a strictly positive integer." exit(198) } if (`alpha' <= 0) { di "{error}Chosen Monte-Carlo confidence factor (alpha) must be strictly positive." exit(198) } ///// Perform main computations //////////////////////////////////////////////// local matalower "" foreach l of local lower { if "`matalower'" == "" local matalower "`l'" else local matalower "`matalower',`l'" } mat lower = (`matalower') local mataupper "" foreach l of local upper { if "`mataupper'" == "" local mataupper "`l'" else local mataupper "`mataupper',`l'" } mat upper = (`mataupper') if (`lendelta' != 0) { local matadelta "" foreach l of local delta { if "`matadelta'" == "" local matadelta "`l'" else local matadelta "`matadelta',`l'" } mat delta = (`matadelta') } else { mat delta = J(1, `lenlower', 0) } mata: mvt_void(`df', "`integrator'", `shifts', `samples', `alpha') ///// Output /////////////////////////////////////////////////////////////////// mat returns = returns return scalar error = returns[2, 1] return scalar integral = returns[1, 1] di "{txt}integral = {res}" returns[1, 1] "{txt}" di "{txt}error = {res}" returns[2, 1] "{txt}" end ///// Mata ///////////////////////////////////////////////////////////////////// mata: void mvt_void(df, integrator, shifts, samples, alpha) { lower = st_matrix("lower") upper = st_matrix("upper") delta = st_matrix("delta") Sigma = st_matrix("sigma") st_matrix("returns", mvt_mata(lower, upper, delta, Sigma, df, integrator, shifts, samples, alpha)) } real scalar integrate(real scalar lower, real scalar upper, real vector arg) { if (((lower == .) & (upper == .)) | ((lower == 0) & (upper == .)) | ((lower != .) & (upper != .))) { return(integrate_main(lower, upper, arg)) } else if ((lower == .) & (upper != .)) { return(integrate_main(0, upper, arg) + integrate_main(0, ., arg)) } else if ((lower != 0) & (upper == .)) { return(integrate_main(lower, 0, arg) + integrate_main(0, ., arg)) } else { return(integrate_main(lower, upper, arg)) } } real matrix integrate_main(real lower, real upper, real vector arg) { if ((lower != .) & (upper != .)) { rw = ( .9997137268, .9984919506, .9962951347, .9931249370, .9889843952, .9838775407, .9778093585, .9707857758, .9628136543, .9539007829, .9440558701, .9332885350, .9216092981, .9090295710, .8955616450, .8812186794, .8660146885, .8499645279, .8330838799, .8153892383, .7968978924, .7776279096, .7575981185, .7368280898, .7153381176, .6931491994, .6702830156, .6467619085, .6226088602, .5978474702, .5725019326, .5465970121, .5201580199, .4932107892, .4657816498, .4378974022, .4095852917, .3808729816, .3517885264, .3223603439, .2926171880, .2625881204, .2323024818, .2017898641, .1710800805, .1402031372, .1091892036, .0780685828, .0468716824, .0156289844, -.0156289844, -.0468716824, -.0780685828, -.1091892036, -.1402031372, -.1710800805, -.2017898641, -.2323024818, -.2625881204, -.2926171880, -.3223603439, -.3517885264, -.3808729816, -.4095852917, -.4378974022, -.4657816498, -.4932107892, -.5201580199, -.5465970121, -.5725019326, -.5978474702, -.6226088602, -.6467619085, -.6702830156, -.6931491994, -.7153381176, -.7368280898, -.7575981185, -.7776279096, -.7968978924, -.8153892383, -.8330838799, -.8499645279, -.8660146885, -.8812186794, -.8955616450, -.9090295710, -.9216092981, -.9332885350, -.9440558701, -.9539007829, -.9628136543, -.9707857758, -.9778093585, -.9838775407, -.9889843952, -.9931249370, -.9962951347, -.9984919506, -.9997137268 \ .0007346345, .0017093927, .0026839254, .0036559612, .0046244501, .0055884280, .0065469485, .0074990733, .0084438715, .0093804197, .0103078026, .0112251140, .0121314577, .0130259479, .0139077107, .0147758845, .0156296211, .0164680862, .0172904606, .0180959407, .0188837396, .0196530875, .0204032326, .0211334421, .0218430024, .0225312203, .0231974232, .0238409603, .0244612027, .0250575445, .0256294029, .0261762192, .0266974592, .0271926134, .0276611982, .0281027557, .0285168543, .0289030896, .0292610841, .0295904881, .0298909796, .0301622651, .0304040795, .0306161866, .0307983790, .0309504789, .0310723374, .0311638357, .0312248843, .0312554235, .0312554235, .0312248843, .0311638357, .0310723374, .0309504789, .0307983790, .0306161866, .0304040795, .0301622651, .0298909796, .0295904881, .0292610841, .0289030896, .0285168543, .0281027557, .0276611982, .0271926134, .0266974592, .0261762192, .0256294029, .0250575445, .0244612027, .0238409603, .0231974232, .0225312203, .0218430024, .0211334421, .0204032326, .0196530875, .0188837396, .0180959407, .0172904606, .0164680862, .0156296211, .0147758845, .0139077107, .0130259479, .0121314577, .0112251140, .0103078026, .0093804197, .0084438715, .0074990733, .0065469485, .0055884280, .0046244501, .0036559612, .0026839254, .0017093927, .0007346345) sum = rw[2, ]:*sf((upper - lower)/2*rw[1, ] :+ (upper + lower)/2, arg) return((upper - lower)/2*quadrowsum(sum)) } else if (lower == 0 & upper == .) { rw = (374.9841128000, 355.2613119000, 339.4351019000, 325.6912634000, 313.3295340000, 301.9858553000, 291.4401336000, 281.5463283000, 272.2011700000, 263.3281685000, 254.8686293000, 246.7762410000, 239.0136298000, 231.5500680000, 224.3598948000, 217.4213933000, 210.7159729000, 204.2275596000, 197.9421331000, 191.8473694000, 185.9323602000, 180.1873909000, 174.6037612000, 169.1736398000, 163.8899456000, 158.7462485000, 153.7366875000, 148.8559014000, 144.0989700000, 139.4613646000, 134.9389050000, 130.5277232000, 126.2242308000, 122.0250921000, 117.9271991000, 113.9276512000, 110.0237356000, 106.2129115000, 102.4927949000, 98.8611460500, 95.3158573500, 91.8549432600, 88.4765308200, 85.1788512100, 81.9602322200, 78.8190913200, 75.7539294700, 72.7633254300, 69.8459306400, 67.0004645200, 64.2257101200, 61.5205102900, 58.8837639800, 56.3144229900, 53.8114889200, 51.3740103300, 49.0010802100, 46.6918335400, 44.4454451100, 42.2611274800, 40.1381290600, 38.0757323700, 36.0732524100, 34.1300351200, 32.2454560000, 30.4189187700, 28.6498541500, 26.9377187300, 25.2819938700, 23.6821847600, 22.1378194500, 20.6484479700, 19.2136415800, 17.8329919500, 16.5061104700, 15.2326276000, 14.0121922500, 12.8444711800, 11.7291484900, 10.6659250900, 9.6545182440, 8.6946611140, 7.7861023780, 6.9286058290, 6.1219500310, 5.3659279860, 4.6603468360, 4.0050275820, 3.3998048270, 2.8445265430, 2.3390538500, 1.8832608260, 1.4770343300, 1.1202738350, .8128912841, .5548109376, .3459691810, .1863141021, .0758036120, .0143861470 \ 7.596410e-96, 3.146290e-94, 9.411750e-96, 7.09970e-106, 3.82990e-102, 2.953620e-97, 2.133670e-97, 2.228220e-97, 1.326430e-95, 5.932280e-99, 1.291530e-99, 1.15800e-100, 1.657940e-95, 4.53820e-100, 1.815580e-90, 1.437780e-93, 2.020360e-91, 1.279760e-88, 6.704800e-86, 2.884880e-83, 1.037900e-80, 3.152470e-78, 8.153750e-76, 1.809860e-73, 3.471850e-71, 5.792630e-69, 8.454960e-67, 1.085370e-64, 1.231380e-62, 1.240240e-60, 1.113570e-58, 8.947340e-57, 6.456270e-55, 4.197750e-53, 2.466810e-51, 1.313980e-49, 6.361270e-48, 2.806040e-46, 1.130480e-44, 4.168860e-43, 1.410140e-41, 4.383790e-40, 1.254830e-38, 3.313060e-37, 8.081630e-36, 1.824200e-34, 3.815860e-33, 7.407430e-32, 1.336210e-30, 2.242650e-29, 3.506270e-28, 5.112370e-27, 6.959200e-26, 8.853230e-25, 1.053590e-23, 1.174020e-22, 1.226010e-21, 1.200850e-20, 1.104090e-19, 9.536250e-19, 7.743090e-18, 5.914430e-17, 4.252610e-16, 2.880120e-15, 1.838350e-14, 1.106500e-13, 6.283520e-13, 3.368210e-12, 1.705060e-11, 8.154800e-11, 3.686330e-10, 1.575600e-09, 6.369710e-09, 2.436430e-08, 8.820060e-08, 3.022640e-07, 9.808340e-07, 3.014270e-06, 8.774310e-06, .0000241958, .0000632109, .0001564521, .0003668548, .0008148716, .0017143197, .0034149800, .0064389510, .0114854424, .0193678281, .0308463086, .0463401336, .0655510093, .0870966385, .1083141121, .1254070908, .1340433397, .1303566130, .1121151033, .0796767462, .0363926059) sum = rw[2, ]:*exp(rw[1, ]):*sf(rw[1, ], arg) return(quadrowsum(sum)) } else if (lower == . & upper == .) { rw = ( 13.4064873400, 12.82379975000, 12.3429642200, 11.9150619400, 11.5214154000, 11.15240439000, 10.8022607500, 10.4671854200, 10.1445099400, 9.83226980800, 9.5289658230, 9.2334208900, 8.9446892170, 8.66199616800, 8.3846969400, 8.1122473110, 7.8441823840, 7.58010080800, 7.3196528220, 7.0625310600, 6.8084633530, 6.55720703200, 6.3085443610, 6.0622788330, 5.8182321350, 5.57624164900, 5.3361583600, 5.0978451050, 4.8611750920, 4.62603063600, 4.3923020790, 4.1598868550, 3.9286886830, 3.69861685900, 3.4695856360, 3.2415136800, 3.0143235800, 2.78794142400, 2.5622964020, 2.3373204640, 2.1129479960, 1.88911553700, 1.6657615090, 1.4428259700, 1.2202503910, .99797743610, .7759507615, .5541148236, .3324146923, .11079587240, -.1107958724, -.3324146923, -.5541148236, -.77595076150, -.9979774361, -1.2202503910, -1.4428259700, -1.6657615090, -1.8891155370, -2.1129479960, -2.3373204640, -2.5622964020, -2.7879414240, -3.0143235800, -3.2415136800, -3.4695856360, -3.6986168590, -3.9286886830, -4.1598868550, -4.3923020790, -4.6260306360, -4.8611750920, -5.0978451050, -5.3361583600, -5.5762416490, -5.8182321350, -6.0622788330, -6.3085443610, -6.5572070320, -6.8084633530, -7.0625310600, -7.3196528220, -7.5801008080, -7.8441823840, -8.1122473110, -8.3846969400, -8.6619961680, -8.9446892170, -9.2334208900, -9.5289658230, -9.8322698080, -10.1445099400, -10.4671854200, -10.8022607500, -11.1524043900, -11.5214154000, -11.9150619400, -12.3429642200, -12.8237997500, -13.4064873400 \ 5.908070e-79, 1.972860e-72, 3.083030e-67, 9.019220e-63, 8.518880e-59, 3.459480e-55, 7.191530e-52, 8.597560e-49, 6.420730e-46, 3.185220e-43, 1.100470e-40, 2.748780e-38, 5.116230e-36, 7.274570e-34, 8.067430e-32, 7.101810e-30, 5.037790e-28, 2.917350e-26, 1.394840e-24, 5.561030e-23, 1.865000e-21, 5.302320e-20, 1.286830e-18, 2.682490e-17, 4.829840e-16, 7.548900e-15, 1.028870e-13, 1.227880e-12, 1.287900e-11, 1.191300e-10, 9.747920e-10, 7.075860e-09, 4.568130e-08, 2.629100e-07, 1.351800e-06, 6.221520e-06, .0000256762, .0000951716, .0003172920, .0009526922, .0025792733, .0063030003, .0139156652, .0277791274, .0501758127, .0820518274, .1215379868, .1631300305, .1984628503, .2188926296, .2188926296, .1984628503, .1631300305, .1215379868, .0820518274, .0501758127, .0277791274, .0139156652, .0063030003, .0025792733, .0009526922, .0003172920, .0000951716, .0000256762, 6.221520e-06, 1.351800e-06, 2.629100e-07, 4.568130e-08, 7.075860e-09, 9.747920e-10, 1.191300e-10, 1.287900e-11, 1.227880e-12, 1.028870e-13, 7.548900e-15, 4.829840e-16, 2.682490e-17, 1.286830e-18, 5.302320e-20, 1.865000e-21, 5.561030e-23, 1.394840e-24, 2.917350e-26, 5.037790e-28, 7.101810e-30, 8.067430e-32, 7.274570e-34, 5.116230e-36, 2.748780e-38, 1.100470e-40, 3.185220e-43, 6.420730e-46, 8.597560e-49, 7.191530e-52, 3.459480e-55, 8.518880e-59, 9.019220e-63, 3.083030e-67, 1.972860e-72, 5.908070e-79) sum = rw[2, ]:*exp(rw[1, ]:^2):*sf(rw[1, ], arg) return(quadrowsum(sum)) } } real colvector mvt_mata(real vector lower, real vector upper, real vector delta, real matrix Sigma, real scalar df, string integrator, real scalar shifts, real scalar samples, real scalar alpha) { if ((df == 0) | (df == .)) { if (integrator == "mvnormal") { for (i = 1; i <= k; i++) { if (lower[i] == .) { lower[i] = -8e307 } if (upper[i] == .) { upper[i] = 8e307 } } return((mvnormalcv(lower, upper, delta, vech(Sigma)') \ .)) } else { return(pmvnormal_mata(lower, upper, delta, Sigma, shifts, samples, alpha)) } } else { k = rows(Sigma) if (k == 1) { if ((lower == .) & (upper == .)) { I = 1 } else if ((lower != .) & (upper == .)) { I = 1 - t(df, lower - delta) } else if ((lower == .) & (upper != .)) { I = t(df, upper - delta) } else { I = t(df, upper - delta) - t(df, lower - delta) } E = 0 } else { a = lower - delta b = upper - delta C = J(k, k, 0) u = y = J(1, k - 1, 0) atilde = (a[1], J(1, k - 2, 0)) btilde = (b[1], J(1, k - 2, 0)) sqrt_Sigma11 = sqrt(Sigma[1, 1]) for (i = 1; i <= k - 1; i++) { args = J(1, k - i + 1, 0) for (j = 1; j <= k - i + 1; j++){ s = j + i - 1 if (i > 1) { y2 = sum(y[1::(i - 1)]:^2) Cy = sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]) C2 = sum(C[s, 1::(i - 1)]:^2) if ((a[s] != .) & (b[s] != .)) { args[j] = t(df + i - 1, sqrt((df + i - 1)/(df + y2))*((b[s] - Cy)/ sqrt(Sigma[s, s] - C2))) - t(df + i - 1, sqrt((df + i - 1)/(df + y2))*((a[s] - Cy)/ sqrt(Sigma[s, s] - C2))) } else if ((a[s] == .) & (b[s] != .)) { args[j] = t(df + i - 1, sqrt((df + i - 1)/(df + y2))*((b[s] - Cy)/ sqrt(Sigma[s, s] - C2))) } else if ((b[s] == .) & (a[s] != .)) { args[j] = 1 - t(df + i - 1, sqrt((df + i - 1)/(df + y2))*((a[s] - Cy)/ sqrt(Sigma[s, s] - C2))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } else { if ((a[s] != .) & (b[s] != .)) { args[j] = t(df, b[s]/sqrt_Sigma11) - t(df, a[s]/sqrt_Sigma11) } else if ((a[s] == .) & (b[s] != .)) { args[j] = t(df, b[s]/sqrt_Sigma11) } else if ((b[s] == .) & (a[s] != .)) { args[j] = 1 - t(df, a[s]/sqrt_Sigma11) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } } ii = ww =. minindex(args, 1, ii, ww) m = i - 1 + ii[1] if (i != m) { tempa = a tempb = b tempa[i] = a[m] tempa[m] = a[i] a = tempa tempb[i] = b[m] tempb[m] = b[i] b = tempb tempSigma = Sigma tempSigma[i, 1::k] = Sigma[m, 1::k] tempSigma[m, 1::k] = Sigma[i, 1::k] Sigma = tempSigma Sigma[1::k, i] = tempSigma[1::k, m] Sigma[1::k, m] = tempSigma[1::k, i] } if (i > 1) { if (i != m) { tempC = C tempC[i, 1::k] = C[m, 1::k] tempC[m, 1::k] = C[i, 1::k] C = tempC C[1::k, i] = tempC[1::k, m] C[1::k, m] = tempC[1::k, i] } C[i, i] = sqrt(Sigma[i, i] - sum(C[i, 1::(i - 1)]:^2)) for (s = i + 1; s <= k; s++){ C[s, i] = (Sigma[s, i] - sum(C[i, 1::(i - 1)]:*C[s, 1::(i - 1)]))/C[i, i] } Cy = sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]) atilde[i] = (a[i] - Cy)/C[i, i] btilde[i] = (b[i] - Cy)/C[i, i] } else { C[i, i] = sqrt(Sigma[i, i]) C[2::k, i] = Sigma[2::k, i]/C[i, i] } arg = (df, i) u[i] = (gamma((df + 1)/2)/(gamma((df + i - 1)/2)*((df + i - 1)*pi())^0.5))* integrate(atilde[i], btilde[i], arg)/ (t(df + i - 1, btilde[i]) - t(df + i - 1, atilde[i])) if (i == 1) { y[i] = u[i] } else { y[i] = u[i]*sqrt((df + sum(y[1::(i - 1)]:^2))/(df + i - 1)) } } C[k, k] = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2)) sqrt_primes = sqrt((2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541)[1::(k - 1)]) I = V = 0 d = J(1, k, 0) e = J(1, k, 1) f = J(1, k, 0) if (a[1] != .) { d[1] = t(df, a[1]/C[1, 1]) } if (b[1] != .) { e[1] = t(df, b[1]/C[1, 1]) } f = (e[1] - d[1], J(1, k - 1, 0)) for (i = 1; i <= shifts; i++) { Ii = 0 Delta = runiform(1, k - 1) for (j = 1; j <= samples; j++) { u[1] = invt(df, d[1] + abs(2*mod(j*sqrt_primes[1] + Delta[1], 1) - 1)* (e[1] - d[1])) y[1] = u[1] for (l = 2; l <= k; l++) { if ((a[l] != .) & (b[l] != .)) { y2 = sum(y[1::(l - 1)]:^2) Cy = sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]) d[l] = t(df + l - 1, sqrt((df + l - 1)/(df + y2))* ((a[l] - Cy)/C[l, l])) e[l] = t(df + l - 1, sqrt((df + l - 1)/(df + y2))* ((b[l] - Cy)/C[l, l])) } else if ((a[l] != .) & (b[l] == .)) { d[l] = t(df + l - 1, sqrt((df + l - 1)/(df + sum(y[1::(l - 1)]:^2)))* ((a[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l])) } else if ((a[l] == .) & (b[l] != .)) { e[l] = t(df + l - 1, sqrt((df + l - 1)/(df + sum(y[1::(l - 1)]:^2)))* ((b[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l])) } f[l] = (e[l] - d[l])*f[l - 1] if (l < k) { u[l] = invt(df + l - 1, d[l] + abs(2*mod(j*sqrt_primes[l] + Delta[l], 1) - 1)* (e[l] - d[l])) y[l] = u[l]*sqrt((df + sum(y[1::(l - 1)]:^2))/(df + l - 1)) } } Ii = Ii + (f[k] - Ii)/j } del = (Ii - I)/i I = I + del V = (i - 2)*V/i + del^2 E = alpha*sqrt(V) } } } return((I \ E)) } real colvector pmvnormal_mata(real vector lower, real vector upper, real vector mean, real matrix Sigma, real scalar shifts, real scalar samples, real scalar alpha) { k = rows(Sigma) if (k == 1) { if ((lower == .) & (upper == .)) { I = 1 } else if ((lower != .) & (upper == .)) { I = 1 - normal((lower - mean)/sqrt(Sigma)) } else if ((lower == .) & (upper != .)) { I = normal((upper - mean)/sqrt(Sigma)) } else { sqrt_Sigma = sqrt(Sigma) I = normal((upper - mean)/sqrt_Sigma) - normal((lower - mean)/sqrt_Sigma) } E = 0 } else if (k == 2) { if (lower[1] == .) lower[1] = -8e307 if (lower[2] == .) lower[2] = -8e307 if (upper[1] == .) upper[1] = 8e307 if (upper[2] == .) upper[2] = 8e307 sqrt_Sigma = sqrt((Sigma[1, 1], Sigma[2, 2])) a = (lower - mean):/(sqrt_Sigma[1], sqrt_Sigma[2]) b = (upper - mean):/(sqrt_Sigma[1], sqrt_Sigma[2]) r = Sigma[1, 2]/(sqrt_Sigma[1]*sqrt_Sigma[2]) I = binormal(b[1], b[2], r) + binormal(a[1], a[2], r) - binormal(b[1], a[2], r) - binormal(a[1], b[2], r) E = 0 } else { a = lower - mean b = upper - mean C = J(k, k, 0) zero_k_min_2 = J(1, k - 2, 0) zero_k_min_1 = (zero_k_min_2, 0) y = zero_k_min_1 atilde = (a[1], zero_k_min_2) btilde = (b[1], zero_k_min_2) sqrt_Sigma11 = sqrt(Sigma[1, 1]) args = J(1, k, 1) for (j = 1; j <= k; j++) { if ((a[j] != .) & (b[j] != .)) { args[j] = normal(b[j]/sqrt_Sigma11) - normal(a[j]/sqrt_Sigma11) } else if ((a[j] == .) & (b[j] != .)) { args[j] = normal(b[j]/sqrt_Sigma11) } else if ((b[j] == .) & (a[j] != .)) { args[j] = 1 - normal(a[j]/sqrt_Sigma11) } } ii = ww = . minindex(args, 1, ii, ww) if (ii[1] != 1) { tempa = a tempb = b tempa[1] = a[ii[1]] tempa[ii[1]] = a[1] a = tempa tempb[1] = b[ii[1]] tempb[ii[1]] = b[1] b = tempb tempSigma = Sigma tempSigma[1, ] = Sigma[ii[1], ] tempSigma[ii[1], ] = Sigma[1, ] Sigma = tempSigma Sigma[, 1] = tempSigma[, ii[1]] Sigma[, ii[1]] = tempSigma[, 1] C[1, 1] = sqrt(Sigma[1, 1]) C[2::k, 1] = Sigma[2::k, 1]/C[1, 1] } else { C[, 1] = (sqrt_Sigma11 \ Sigma[2::k, 1]/sqrt_Sigma11) } if (atilde[1] != btilde[1]) { y[1] = (normalden(atilde[1]) - normalden(btilde[1]))/ (normal(btilde[1]) - normal(atilde[1])) } for (i = 2; i <= k - 1; i++) { args = J(1, k - i + 1, 1) i_vec = 1::(i - 1) for (j = 1; j <= k - i + 1; j++) { s = j + i - 1 if ((a[s] != .) & (b[s] != .)) { Cy = sum(C[s, i_vec]:*y[i_vec]) denom = sqrt(Sigma[s, s] - sum(C[s, i_vec]:^2)) args[j] = normal((b[s] - Cy)/denom) - normal((a[s] - Cy)/denom) } else if ((a[s] == .) & (b[s] != .)) { args[j] = normal((b[s] - sum(C[s, i_vec]:*y[i_vec]))/ sqrt(Sigma[s, s] - sum(C[s, i_vec]:^2))) } else if ((b[s] == .) & (a[s] != .)) { args[j] = 1 - normal((a[s] - sum(C[s, i_vec]:*y[i_vec]))/ sqrt(Sigma[s, s] - sum(C[s, i_vec]:^2))) } } ii = ww = . minindex(args, 1, ii, ww) m = i - 1 + ii[1] if (i != m) { tempa = a tempb = b tempa[i] = a[m] tempa[m] = a[i] a = tempa tempb[i] = b[m] tempb[m] = b[i] b = tempb tempSigma = Sigma tempSigma[i, ] = Sigma[m, ] tempSigma[m, ] = Sigma[i, ] Sigma = tempSigma Sigma[, i] = tempSigma[, m] Sigma[, m] = tempSigma[, i] tempC = C tempC[i, ] = C[m, ] tempC[m, ] = C[i, ] C = tempC C[, i] = tempC[, m] C[, m] = tempC[, i] } C[i, i] = sqrt(Sigma[i, i] - sum(C[i, i_vec]:^2)) i_vec2 = (i + 1)::k C[i_vec2, i] = (Sigma[i_vec2, i] - rowsum(J(k - i, 1, C[i, i_vec]):*C[i_vec2, i_vec]))/ C[i, i] Cy = sum(C[i, i_vec]:*y[i_vec]) atilde[i] = (a[i] - Cy)/C[i, i] btilde[i] = (b[i] - Cy)/C[i, i] if (atilde[i] != btilde[i]) { y[i] = (normalden(atilde[i]) - normalden(btilde[i]))/ (normal(btilde[i]) - normal(atilde[i])) } } C[k, k] = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2)) I = V = 0 if (a[1] != .) { d = J(samples, 1, (normal(a[1]/C[1, 1]), J(1, k - 1, 0))) } else { d = J(samples, k, 0) } if (b[1] != .) { e = J(samples, 1, (normal(b[1]/C[1, 1]), J(1, k - 1, 1))) } else { e = J(samples, 1, J(1, k, 1)) } f = (e[, 1] - d[, 1], J(samples, k - 1, 0)) y = J(samples, k - 1, 0) Delta = runiform(shifts, k - 1) samples_sqrt_primes = (1::samples)*sqrt((2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541)[1::(k - 1)]) Ii = J(1, shifts, 0) for (i = 1; i <= shifts; i++) { for (l = 2; l <= k; l++) { l_vec = 1::(l - 1) y[, l - 1] = invnormal(d[, l - 1] + abs(2*mod(samples_sqrt_primes[, l - 1] :+ Delta[i, l - 1], 1) :- 1):* (e[, l - 1] - d[, l - 1])) if ((a[l] != .) & (b[l] != .)) { Cy = rowsum(J(samples, 1, C[l, l_vec]):*y[, l_vec]) d[, l] = normal((a[l] :- Cy)/C[l, l]) e[, l] = normal((b[l] :- Cy)/C[l, l]) f[, l] = (e[, l] :- d[, l]):*f[, l - 1] } else if ((a[l] != .) & (b[l] == .)) { d[, l] = normal((a[l] :- rowsum(J(samples, 1, C[l, l_vec]):*y[, l_vec]))/C[l, l]) f[, l] = (1 :- d[, l]):*f[, l - 1] } else if ((a[l] == .) & (b[l] != .)) { e[, l] = normal((b[l] :- rowsum(J(samples, 1, C[l, l_vec]):*y[, l_vec]))/C[l, l]) f[, l] = e[, l]:*f[, l - 1] } else { f[, l] = f[, l - 1] } } for (j = 1; j <= samples; j++) { Ii[i] = Ii[i] + (f[j, k] - Ii[i])/j } del = (Ii[i] - I)/i I = I + del V = (i - 2)*V/i + del^2 E = alpha*sqrt(V) } } return((I \ E)) } real vector sf(real rowvector s, real vector arg) { return(s:*(1 :+ (s:^2)/(arg[1] + arg[2] - 1)):^(-(arg[1] + arg[2])/2)) } end