{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 1 10 128 0 0 1 0 2 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 1 11 0 0 0 1 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 0 0 128 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 0 128 0 0 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "2D Output" 2 20 "" 0 0 0 0 128 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "Non-proportional" -1 256 "Courier" 1 11 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 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 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 6 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1 " -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 24 6 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times" 1 15 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 20 4 3 0 3 0 2 2 0 1 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 3 0 0 6 0 3 36 1 0 2 2 0 1 }{PSTYLE "Maple \+ Output" -1 12 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 3 0 0 6 0 3 36 1 0 2 2 0 1 }{PSTYLE "Maple Plot" -1 13 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 6 0 1 0 1 0 2 2 0 1 }{PSTYLE "Title" -1 18 1 {CSTYLE "" -1 -1 "Times" 1 22 0 128 128 1 2 1 2 2 2 2 1 1 1 1 }3 1 0 0 14 12 1 0 1 0 2 2 19 1 } {PSTYLE "Author" -1 19 1 {CSTYLE "" -1 -1 "Times" 1 16 0 0 128 1 2 1 2 2 2 2 1 1 1 1 }3 1 0 0 8 0 1 0 1 0 2 2 261 1 }{PSTYLE "Address" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 128 1 1 2 2 2 2 2 1 1 1 1 }3 1 0 0 6 0 1 0 1 0 2 2 263 1 }{PSTYLE "Date" -1 257 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 128 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 8 6 1 0 1 0 2 2 0 1 }{PSTYLE "Address - email" -1 258 1 {CSTYLE "" -1 -1 "Verdana" 1 10 0 0 128 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 4 0 1 0 1 0 2 2 262 1 }{PSTYLE "Normal Indented" -1 259 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 18 6 0 1 0 1 0 2 2 0 1 }{PSTYLE "Address" -1 260 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 128 1 1 2 2 2 2 2 1 1 1 1 }3 1 0 0 6 0 1 0 1 0 2 2 263 1 }{PSTYLE "Date" -1 261 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 128 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 8 6 1 0 1 0 2 2 0 1 }{PSTYLE "Address - email" -1 262 1 {CSTYLE "" -1 -1 "Verdana" 1 10 0 0 128 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 4 0 1 0 1 0 2 2 262 1 }{PSTYLE "Normal Indented" -1 263 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 18 6 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 125 "A Method for Directly Ge nerating a Gaussian Distribution with Nonunit Variance and Nonzero Mea n from Uniform Random Deviates " }}{PARA 19 "" 0 "" {TEXT -1 15 "Marc \+ A. Murison" }}{PARA 260 "" 0 "" {TEXT -1 74 "Astronomical Applications Department\nU.S. Naval Observatory\nWashington, DC" }}{PARA 262 "" 0 "" {TEXT -1 57 "murison@aa.usno.navy.mil\nhttp://aa.usno.navy.mil/muri son/" }}{PARA 261 "" 0 "" {TEXT -1 11 "8 May, 2000" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 15 "1. Introduction" }}{EXCHG {PARA 263 "" 0 "" {TEXT -1 196 "It is universally the case that, when an algorithm is pr esented that takes uniform random deviates and produces a Gaussian dis tribution, the Gaussian distribution has zero mean and unit variance: \+ " }{XPPEDIT 18 0 "G(x) = exp(-x^2/2)/sqrt(2*Pi);" "6#/-%\"GG6#%\"xG*&- %$expG6#,$*&F'\"\"#F.!\"\"F/\"\"\"-%%sqrtG6#*&F.F0%#PiGF0F/" }{TEXT -1 20 ". For nonzero mean " }{XPPEDIT 18 0 "mu;" "6#%#muG" }{TEXT -1 22 " and nonunit variance " }{XPPEDIT 18 0 "sigma^2;" "6#*$%&sigmaG\" \"#" }{TEXT -1 31 ", the Gaussian distribution is " }{XPPEDIT 18 0 "G( x,mu,sigma) = exp(-(x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*Pi));" "6#/-%\" GG6%%\"xG%#muG%&sigmaG*&-%$expG6#,$*&,&F'\"\"\"F(!\"\"\"\"#*&F3F1*$F)F 3F1F2F2F1*&F)F1-%%sqrtG6#*&F3F1%#PiGF1F1F2" }{TEXT -1 21 ", with norma lization " }{XPPEDIT 18 0 "int(G(x,mu,sigma),x = -infinity .. infinity ) = 1;" "6#/-%$intG6$-%\"GG6%%\"xG%#muG%&sigmaG/F*;,$%)infinityG!\"\"F 0\"\"\"" }{TEXT -1 156 ". How does one obtain a distribution with a g iven variance from a distribution with unit variance? The Gaussian fu nction has a useful invariant relation: " }{XPPEDIT 18 0 "G(x,mu,sigma ) = G(y,0,1)/sigma;" "6#/-%\"GG6%%\"xG%#muG%&sigmaG*&-F%6%%\"yG\"\"!\" \"\"F/F)!\"\"" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "y = (x-mu)/sigma ;" "6#/%\"yG*&,&%\"xG\"\"\"%#muG!\"\"F(%&sigmaGF*" }{TEXT -1 46 ". He nce, if we have a set of random deviates " }{TEXT 257 1 "y" }{TEXT -1 205 " which are distributed as a Gaussian with unit variance and zero \+ mean, we may convert those deviates to a Gaussian distribution with no nzero mean and nonunit variance by the simple trick of multiplying by \+ " }{XPPEDIT 18 0 "sigma;" "6#%&sigmaG" }{TEXT -1 12 " and adding " } {XPPEDIT 18 0 "mu;" "6#%#muG" }{TEXT -1 82 ". However, one might stil l wish for an algorithm that takes uniform deviates and " }{TEXT 258 8 "directly" }{TEXT -1 63 " produces the desired nonunit variance Gaus sian distribution. " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 34 "2. A Modified Box-Muller Algorithm" }}{EXCHG {PARA 263 " " 0 "" {TEXT -1 271 "The usual algorithm for producing a Gaussian dist ribution from uniform random deviates makes use of the Box-Muller tran sformation. We will modify the usual algorithm so that the Box-Muelle r transformation produces the desired nonunit-variance Gaussian distri bution. If " }{XPPEDIT 18 0 "z[1];" "6#&%\"zG6#\"\"\"" }{TEXT -1 5 " \+ and " }{XPPEDIT 18 0 "z[2];" "6#&%\"zG6#\"\"#" }{TEXT -1 78 " are inde pendent and uniformly distributed in the range 0 to 1, then consider \+ " }{XPPEDIT 18 0 "y[1];" "6#&%\"yG6#\"\"\"" }{TEXT -1 5 " and " } {XPPEDIT 18 0 "y[2];" "6#&%\"yG6#\"\"#" }{TEXT -1 10 ", given by" }}} {EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "restart;" "6#%(restartG" }}} {EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "y[1] = sqrt(-2*ln(z[1]))*cos(2*P i*z[2]);" "6#/&%\"yG6#\"\"\"*&-%%sqrtG6#,$*&\"\"#F'-%#lnG6#&%\"zG6#F'F '!\"\"F'-%$cosG6#*(F.F'%#PiGF'&F36#F.F'F'" }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "y[2] = sqrt(-2*ln(z[1]))*sin(2*Pi*z[2]);" "6#/&%\"yG6#\"\"#*&-%% sqrtG6#,$*&F'\"\"\"-%#lnG6#&%\"zG6#F.F.!\"\"F.-%$sinG6#*(F'F.%#PiGF.&F 36#F'F.F." }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%\"yG6#\"\"\"*&-%%sqrtG 6#,$-%#lnG6#&%\"zGF&!\"#F'-%$cosG6#,$*&%#PiGF'&F16#\"\"#F'F;F'" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/&%\"yG6#\"\"#*&-%%sqrtG6#,$-%#lnG6#&% \"zG6#\"\"\"!\"#F3-%$sinG6#,$*&%#PiGF3&F1F&F3F'F3" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 81 "This is the Box-Muller transformation, which conve rts a 2D uniform distribution (" }{XPPEDIT 18 0 "z[1],z[2];" "6$&%\"zG 6#\"\"\"&F$6#\"\"#" }{TEXT -1 40 ") to a bivariate Gaussian distributi on (" }{XPPEDIT 18 0 "y[1],y[2];" "6$&%\"yG6#\"\"\"&F$6#\"\"#" }{TEXT -1 199 ") with zero mean and unit variance. (The proof that this tran sformation yields the desired Gaussian distribution is a special case \+ of the proof of the modified transformation, which is given below.)" } }}{EXCHG {PARA 263 "" 0 "" {TEXT -1 83 "The trick to making the Box-Mu ller transformation produce a Gaussian distribution (" }{XPPEDIT 18 0 "x[1],x[2];" "6$&%\"xG6#\"\"\"&F$6#\"\"#" }{TEXT -1 82 ") with nonunit variance is to use the invariant relation mentioned above. Setting" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "isolate(y = (x-mu)/sigma,x); " "6#-%(isolateG6$/%\"yG*&,&%\"xG\"\"\"%#muG!\"\"F+%&sigmaGF-F*" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"xG,&*&%\"yG\"\"\"%&sigmaGF(F(%#muG F(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 7 "we have" }}}{EXCHG {PARA 0 " " 0 "" {XPPEDIT 19 1 "x[1] = sqrt(-2*ln(z[1]))*cos(2*Pi*z[2])*sigma+mu ;" "6#/&%\"xG6#\"\"\",&*(-%%sqrtG6#,$*&\"\"#F'-%#lnG6#&%\"zG6#F'F'!\" \"F'-%$cosG6#*(F/F'%#PiGF'&F46#F/F'F'%&sigmaGF'F'%#muGF'" }}{PARA 0 " " 0 "" {XPPEDIT 19 1 "x[2] = sqrt(-2*ln(z[1]))*sin(2*Pi*z[2])*sigma+mu ;" "6#/&%\"xG6#\"\"#,&*(-%%sqrtG6#,$*&F'\"\"\"-%#lnG6#&%\"zG6#F/F/!\" \"F/-%$sinG6#*(F'F/%#PiGF/&F46#F'F/F/%&sigmaGF/F/%#muGF/" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#/&%\"xG6#\"\"\",&*(-%%sqrtG6#,$-%#lnG6#&%\"zGF&! \"#F'-%$cosG6#,$*&%#PiGF'&F26#\"\"#F'F&%\"XG6#\"\"\",&*(-%% sqrtG6#,$*&\"\"#F'-%#lnG6#&%\"zG6#F'F'!\"\"F'-%$cosG6#*(F/F'%#PiGF'&F4 6#F/F'F'%&sigmaGF'F'%#muGF'" }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "X[2] := \+ sqrt(-2*ln(z[1]))*sin(2*Pi*z[2])*sigma+mu:" "6#>&%\"XG6#\"\"#,&*(-%%sq rtG6#,$*&F'\"\"\"-%#lnG6#&%\"zG6#F/F/!\"\"F/-%$sinG6#*(F'F/%#PiGF/&F46 #F'F/F/%&sigmaGF/F/%#muGF/" }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "sqrt(((x[ 1]-mu)/sigma)^2+((x[2]-mu)/sigma)^2) = simplify(sqrt(((X[1]-mu)/sigma) ^2+((X[2]-mu)/sigma)^2));" "6#/-%%sqrtG6#,&*$*&,&&%\"xG6#\"\"\"F.%#muG !\"\"F.%&sigmaGF0\"\"#F.*$*&,&&F,6#F2F.F/F0F.F1F0F2F.-%)simplifyG6#-F% 6#,&*$*&,&&%\"XG6#F.F.F/F0F.F1F0F2F.*$*&,&&FB6#F2F.F/F0F.F1F0F2F." }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/*$-%%sqrtG6#,&*&*$),&&%\"xG6#\"\"\"F0 %#muG!\"\"\"\"#F0F0*$)%&sigmaGF3F0F2F0*&*$),&&F.6#F3F0F1F2F3F0F0*$F5F0 F2F0F0*&-F&6#F3F0-F&6#,$-%#lnG6#&%\"zGF/F2F0" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 19 "which we solve for " }{XPPEDIT 18 0 "z[1]" "6#&%\"zG6# \"\"\"" }{TEXT -1 1 "," }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "isola te(%,z[1]):" "6#-%(isolateG6$%\"%G&%\"zG6#\"\"\"" }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "location(%,select(has,op(rhs(%)),mu)):" "6#-%)locationG 6$%\"%G-%'selectG6%%$hasG-%#opG6#-%$rhsG6#F&%#muG" }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "subsop(% = csquare(op(%,%%),[x[1], x[2]]),%%):" "6#-%'s ubsopG6$/%\"%G-%(csquareG6$-%#opG6$F'%#%%G7$&%\"xG6#\"\"\"&F16#\"\"#F. " }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "subs(A = x[1]-mu,B = x[2]-mu,expand (subs(x[1] = mu+A,x[2] = mu+B,%)));" "6#-%%subsG6%/%\"AG,&&%\"xG6#\"\" \"F,%#muG!\"\"/%\"BG,&&F*6#\"\"#F,F-F.-%'expandG6#-F$6%/&F*6#F,,&F-F,F 'F,/&F*6#F4,&F-F,F0F,%\"%G" }}{PARA 0 "" 0 "" {XPPEDIT 19 1 "eqn_z1 := %:" "6#>%'eqn_z1G%\"%G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%\"zG6#\" \"\"*&-%$expG6#,$*&*$),&&%\"xG6#\"\"#F'%#muG!\"\"F4F'F'*$)%&sigmaGF4F' F6#F6F4F'-F*6#,$*&*$),&&F2F&F'F5F6F4F'F'*$F8F'F6F:F'" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 11 "and we have" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "x[2]/x[1] = X[2]/X[1];" "6#/*&&%\"xG6#\"\"#\"\"\"&F&6#F )!\"\"*&&%\"XG6#F(F)&F/6#F)F," }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/*&&% \"xG6#\"\"#\"\"\"&F&6#F)!\"\"*&,&*(-%%sqrtG6#,$-%#lnG6#&%\"zGF+!\"#F)- %$sinG6#,$*&%#PiGF)&F8F'F)F(F)%&sigmaGF)F)%#muGF)F),&*(F0F)-%$cosGFF'F:)F@F'F:F:*,F'F:)F5F'F:FFF:,&*$FFF:F:*$)FF:F@F:F:*&F%'eqn_z2G%\"%G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/&%\"zG6#\"\"#,$*&-%'arctanG6#*&,&&%\"xGF&!\"\"%#m uG\"\"\"F3,&&F06#F3F1F2F3F1F3%#PiGF1#F3F'" }}}{EXCHG {PARA 263 "" 0 " " {TEXT -1 69 "Now, the fundamental transformation law of probabilitie s states that " }{XPPEDIT 18 0 "abs(p(y)*dy) = abs(p(x)*dx);" "6#/-%$a bsG6#*&-%\"pG6#%\"yG\"\"\"%#dyGF,-F%6#*&-F)6#%\"xGF,%#dxGF," }{TEXT -1 5 ", or " }{XPPEDIT 18 0 "p(y) = p(x)*abs(dx/dy);" "6#/-%\"pG6#%\"y G*&-F%6#%\"xG\"\"\"-%$absG6#*&%#dxGF,%#dyG!\"\"F," }{TEXT -1 8 ", wher e " }{XPPEDIT 18 0 "p(x)*dx;" "6#*&-%\"pG6#%\"xG\"\"\"%#dxGF(" }{TEXT -1 25 " is the probability that " }{TEXT 259 1 "x" }{TEXT -1 14 " lies between " }{TEXT 260 1 "x" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "x+dx; " "6#,&%\"xG\"\"\"%#dxGF%" }{TEXT -1 6 ", and " }{XPPEDIT 18 0 "p(y)*d y" "6#*&-%\"pG6#%\"yG\"\"\"%#dyGF(" }{TEXT -1 25 " is the probability \+ that " }{TEXT 261 1 "y" }{TEXT -1 14 " lies between " }{TEXT 262 1 "y " }{TEXT -1 5 " and " }{XPPEDIT 18 0 "y+dy;" "6#,&%\"yG\"\"\"%#dyGF%" }{TEXT -1 28 ". The extension of this to " }{TEXT 263 1 "n" }{TEXT -1 15 " dimensions is " }{XPPEDIT 18 0 "p(Y)*dY = p(X)*J(X,Y)*dY;" "6# /*&-%\"pG6#%\"YG\"\"\"%#dYGF)*(-F&6#%\"XGF)-%\"JG6$F.F(F)F*F)" }{TEXT -1 8 ", where " }{XPPEDIT 18 0 "p(X)" "6#-%\"pG6#%\"XG" }{TEXT -1 42 " is the joint probability distribution of " }{TEXT 264 1 "X" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "p(Y)" "6#-%\"pG6#%\"YG" }{TEXT -1 42 " is th e joint probability distribution of " }{TEXT 265 1 "Y" }{TEXT -1 2 ", \+ " }{XPPEDIT 18 0 "X = (x[1], x[2] .. x[n]);" "6#/%\"XG6$&%\"xG6#\"\"\" ;&F'6#\"\"#&F'6#%\"nG" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "Y = (y[1], y[2 ] .. y[n]);" "6#/%\"YG6$&%\"yG6#\"\"\";&F'6#\"\"#&F'6#%\"nG" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "dY = dy[1]*dy[2] .. dy[n];" "6#/%#dYG;*&&%#d yG6#\"\"\"F*&F(6#\"\"#F*&F(6#%\"nG" }{TEXT -1 6 ", and " }{XPPEDIT 18 0 "J(X,Y)" "6#-%\"JG6$%\"XG%\"YG" }{TEXT -1 49 " is the Jacobi determi nant. For two dimensions, " }{XPPEDIT 18 0 "J(X,Y) = det(matrix([[dif f(x[1],y[1]), diff(x[1],y[2])], [diff(x[2],y[1]), diff(x[2],y[2])]])); " "6#/-%\"JG6$%\"XG%\"YG-%$detG6#-%'matrixG6#7$7$-%%diffG6$&%\"xG6#\" \"\"&%\"yG6#F7-F26$&F56#F7&F96#\"\"#7$-F26$&F56#FA&F96#F7-F26$&F56#FA& F96#FA" }{TEXT -1 10 ". Hence, " }}}{EXCHG {PARA 0 "" 0 "" {MPLTEXT 1 0 645 "JacobiDeriv := proc( eqs::\{name=algebraic,list(name=algebrai c)\}, vars::\{name,list(name)\} )\n local M, Nrows, Ncols, k, j, expr , var;\n if type(eqs,name=algebraic) then\n Nrows := 1;\n else\n \+ Nrows := nops(eqs);\n end if;\n if type(vars,name) then\n Ncol s := 1;\n else\n Ncols := nops(vars);\n end if;\n M := matrix(Nr ows,Ncols);\n for k from 1 to Nrows do\n if Nrows=1 then\n ex pr := rhs(eqs);\n else\n expr := rhs(eqs[k]);\n end if;\n \+ for j from 1 to Ncols do\n if Ncols=1 then\n var := vars ;\n else\n var := vars[j];\n end if;\n M[k,j] := diff( expr, var );\n od;\n od;\n det(M);\nend proc:" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "J(Z,X) = JacobiDeriv([eqn_z1, eqn_z2],[ x[1], x[2]]);" "6#/-%\"JG6$%\"ZG%\"XG-%,JacobiDerivG6$7$%'eqn_z1G%'eqn _z2G7$&%\"xG6#\"\"\"&F16#\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-% \"JG6$%\"ZG%\"XG,$*&*&-%$expG6#,$*&*$),&&%\"xG6#\"\"\"!\"\"%#muGF7\"\" #F7F7*$)%&sigmaGF:F7F8#F8F:F7-F-6#,$*&*$),&&F56#F:F8F9F7F:F7F7*$FF7F7*&%#PiGF7F" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "This s ays that for a uniform distribution (" }{XPPEDIT 18 0 "z[1],z[2];" "6$ &%\"zG6#\"\"\"&F$6#\"\"#" }{TEXT -1 23 ") the distribution of (" } {XPPEDIT 18 0 "x[1],x[2];" "6$&%\"xG6#\"\"\"&F$6#\"\"#" }{TEXT -1 59 " ) is a symmetric bivariate Gaussian distribution with mean " } {XPPEDIT 18 0 "mu;" "6#%#muG" }{TEXT -1 14 " and variance " }{XPPEDIT 18 0 "sigma^2;" "6#*$%&sigmaG\"\"#" }{TEXT -1 3 ". " }{TEXT 266 5 "Q. E.D" }{TEXT -1 1 "." }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 18 "3. Nume rical Tests" }}{EXCHG {PARA 0 "" 0 "" {MPLTEXT 1 0 767 "Gaussian := mo dule()\n option package;\n export BM_kernel, BM, G, init;\n\n init \+ := proc()\n global __y2, __used;\n __y2 := 0;\n __used := t rue;\n end proc;\n\n BM_kernel := proc( mu, sigma )\n local x1, x 2, y1, r;\n global __used, __y2;\n if __used=true then\n r \+ := rand(1000000)/1000000;\n x1 := r();\n x2 := \+ 2*Pi*r();\n y1 := evalhf(sqrt(-2*ln(x1))*cos(x2)*sigma+mu);\n __y2 := evalhf(sqrt(-2*ln(x1))*sin(x2)*sigma+mu);\n __use d := false;\n y1;\n else\n __used := true;\n __y2;\n end if;\n end proc;\n\n BM := proc( mu, sigma, n::posint )\n \+ [seq(BM_kernel(mu,sigma),k=1..n)];\n end proc;\n\n G := proc( x, mu, sigma )\n exp(-(x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*Pi));\n end pr oc;\n\nend module:" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "with(Gaus sian):" "6#-%%withG6#%)GaussianG" }}}{EXCHG {PARA 263 "" 0 "" {TEXT -1 46 "In the plot below, the red curve is a plot of " }{XPPEDIT 18 0 "G(x,0,1)" "6#-%\"GG6%%\"xG\"\"!\"\"\"" }{TEXT -1 143 ", while the his togram is the result of running 5000 uniform random deviates through t he modified Box-Muller transformation just described with " }{XPPEDIT 18 0 "mu = 0;" "6#/%#muG\"\"!" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "sig ma = 1;" "6#/%&sigmaG\"\"\"" }{TEXT -1 3 ". " }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "with(stats[statplots]):" "6#-%%withG6#&%&statsG6#%*s tatplotsG" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "h1 := histogram(BM (0, 1, 5000),color = blue):" "6#>%#h1G-%*histogramG6$-%#BMG6%\"\"!\"\" \"\"%+]/%&colorG%%blueG" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "g1 : = plot(G(x,0,1),x = -4 .. 4,color = red):" "6#>%#g1G-%%plotG6%-%\"GG6% %\"xG\"\"!\"\"\"/F+;,$\"\"%!\"\"F1/%&colorG%$redG" }}}{EXCHG {PARA 0 " " 0 "" {XPPEDIT 19 1 "display(h1,g1,view = [-6 .. 10, 0 .. .4],labels \+ = [\"x\", \"pdf\"]);" "6#-%(displayG6&%#h1G%#g1G/%%viewG7$;,$\"\"'!\" \"\"#5;\"\"!$\"\"%F./%'labelsG7$Q\"x6\"Q$pdfF8" }}{PARA 13 "" 1 "" {GLPLOT2D 597 371 371 {PLOTDATA 2 "6.-%)POLYGONSG647&7$$!3#)4WShmdBO!# <\"\"!7$F($\"+\"el]J$!#67$$!+]7\"*o9!\"*F-7$F1F+7&F47$F1$\"+8?d#)=!#57 $$!+A6\\*3\"F3F77$F;F+7&F=7$F;$\"+<#*HMEF97$$!+s$GM=)F9F@7$FCF+7&FE7$F C$\"+n(44N$F97$$!+)46=0'F9FH7$FKF+7&FM7$FK$\"+XIRZKF97$$!+SXC_QF9FP7$F SF+7&FU7$FS$\"+)zB7Y$F97$$!+#)[c)y\"F9FX7$FenF+7&Fgn7$Fen$\"+o)f$\\QF9 7$$\"+'R@Jq'!#7Fjn7$F]oF+7&F`o7$F]o$\"+h718(4qp^: !#?7$$!3_LLL[9cgJF*$\"3m:$=!\\gn-FF\\v7$$!3smmmhN2-IF*$\"3(oz$H#QbVS%F \\v7$$!3!******\\`oz$GF*$\"3>%*y>uQ+7rF\\v7$$!3!omm;)3DoEF*$\"3[u2T:Iz M6!#>7$$!3?+++:v2*\\#F*$\"3#QF&Hlj(ov\"Faw7$$!3BLLL8>1DBF*$\"37K)y;=(> tEFaw7$$!3kmmmw))yr@F*$\"3?J*=**p/Jx$Faw7$$!3;+++S(R#**>F*$\"3h'*z_n2J 2aFaw7$$!30++++@)f#=F*$\"3(H>2>'=uJvFaw7$$!3-+++gi,f;F*$\"3DVZ_Ar]25!# =7$$!3qmmm\"G&R2:F*$\"3AluPQI&3G\"F`y7$$!3XLLLtK5F8F*$\"3W*o#*fH\\Pl\" F`y7$$!3_LLL$yP2D\"F*$\"3Jf8f!p1[#=F`y7$$!3eLLL$HsV<\"F*$\"3/Z\")eH$\\ =+#F`y7$$!3!pmmT2Tb3\"F*$\"3v&4AHT;K@#F`y7$$!3+-++]&)4n**F`y$\"3:>`8/% owU#F`y7$$!3cpmmTnt/GF`y7$$!3G)*****\\&y!pmF`y$\"3gV-2cj'R>$F`y7$$!3Y******\\O3E]F`y $\"3!)Q[@DF0;NF`y7$$!3NKLLL3z6LF`y$\"3mvD88W`wPF`y7$$!3/LLLeGmCDF`y$\" 3W7_?lgGkQF`y7$$!3sLLL$)[`P\"=bRF`y7$$!3m0++]-6&)))Faw$\"3-$*[dclqtRF`y7$$!3')RLLe4**RYFaw$ \"3w+R`(fI^)RF`y7$$!3gSnmmmr[RF\\v$\"3Q!oq=q\"R*)RF`y7$$\"3n=L$3Fr)4=F aw$\"33%=mIYp())RF`y7$$\"3S6LL3Uh9SFaw$\"3;wKH(>4i)RF`y7$$\"38/L$e9d$> iFaw$\"3#)z=9OYr\")RF`y7$$\"3'oHLL3+TU)Faw$\"371PdEBHvRF`y7$$\"3AGL$ef eLG\"F`y$\"3w8!oZi/n&RF`y7$$\"3yELL$=2Vs\"F`y$\"3<>dR-SbIRF`y7$$\"3Khm mm7+#\\#F`y$\"3!e7Q>A`u'QF`y7$$\"3)e*****\\`pfKF`y$\"3A$zPDL/Iy$F`y7$$ \"36HLLLm&z\"\\F`y$\"3n_?Goi+NNF`y7$$\"3>(******z-6j'F`y$\"3PF?6iE/-KF `y7$$\"3q\"******4#32$)F`y$\"33fPs\")GGDGF`y7$$\"3q#****\\AF`y7$$\"3G******H%=H<\"F*$\"3p%oH4/o_+#F`y7$$\"3 qKLLo,\"QD\"F*$\"3R<_l%H*z<=F`y7$$\"35mmm1>qM8F*$\"3p(R**\\94rj\"F`y7$ $\"3%)*******HSu]\"F*$\"3CU%*z?hw!G\"F`y7$$\"3'HLL$ep'Rm\"F*$\"3;wN3\" >QD***Faw7$$\"3')******R>4N=F*$\"3SH<6IZ=2uFaw7$$\"3#emm;@2h*>F*$\"30. !*QjXDTaFaw7$$\"3]*****\\c9W;#F*$\"3B#yB)**[\"R$QFaw7$$\"3Lmmmmd'*GBF* $\"3OH:[qH-\\EFaw7$$\"3j*****\\iN7]#F*$\"3;SLCoCUZ:nEF *$\"3S,!f%4Y7Q6Faw7$$\"35LLL.a#o$GF*$\"3v/N=S#4^8(F\\v7$$\"3ammm^Q40IF *$\"3!y=^:5\"ekVF\\v7$$\"3y******z]rfJF*$\"376hw7W\"*4FF\\v7$$\"3gmmmc %GpL$F*$\"3([:!=tr*Q_\"F\\v7$$\"3/LLL8-V&\\$F*$\"3=#zn%[cVn))Fgt7$$\"3 =+++XhUkOF*$\"3V'))p3H7B%[Fgt7$$\"3=+++:o%#h2G-%*histogramG6$-%#BMG6%\"\"#-%%sqrt G6#\"\"&\"%+]/%&colorG%%blueG" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "g2 := plot(G(x,2,sqrt(5)),x = -6 .. 10,color = red):" "6#>%#g2G-%%p lotG6%-%\"GG6%%\"xG\"\"#-%%sqrtG6#\"\"&/F+;,$\"\"'!\"\"\"#5/%&colorG%$ redG" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "display(h2,g2,view = [- 6 .. 10, 0 .. .4],labels = [\"x\", \"pdf\"]);" "6#-%(displayG6&%#h2G%# g2G/%%viewG7$;,$\"\"'!\"\"\"#5;\"\"!$\"\"%F./%'labelsG7$Q\"x6\"Q$pdfF8 " }}{PARA 13 "" 1 "" {GLPLOT2D 597 371 371 {PLOTDATA 2 "6.-%)POLYGONSG 647&7$$!3'eu\")fyG_V&!#<\"\"!7$F($\"+C-0]B(F\"[<\"F3FP7$FSF+7&FU7$FS$\"+JDE*o\"F<7$$\"+i2l(f\"F3FX7$FenF +7&Fgn7$Fen$\"+SeVzF3Fjn7$F]oF+7&F_o7$F]o$\"+X#z9l\"F< 7$$\"+wZdJCF3Fbo7$FeoF+7&Fgo7$Feo$\"+*=nxo\"F<7$$\"+!>)yaGF3Fjo7$F]pF+ 7&F_p7$F]p$\"+r+f`:F<7$$\"+)y_XJ$F3Fbp7$FepF+7&Fgp7$Fep$\"+fx^m8F<7$$ \"+CyDPQF3Fjp7$F]qF+7&F_q7$F]q$\"+nDbk6F<7$$\"+wTh]WF3Fbq7$FeqF+7&Fgq7 $Feq$\"+'>u66)F/7$$\"+]NBJ`F3Fjq7$F]rF+7&F_r7$F]r$\"+N['yx\"F/7$$\"3)H wA%\\R*)[$*F*Fbr7$FerF+-%'COLOURG6&%$RGBG$F+F+F\\s$\"*++++\"!\")-%*THI CKNESSG6#\"\"#-%'SYMBOLG6$%'CIRCLEG\"\")-%&STYLEG6#%%LINEG-%'CURVESG6' 7hn7$$!\"'F+$\"3s`'Q9SCW'H!#@7$$!3OLLLLbC^cF*$\"3J1\\$R$))y;^Fft7$$!3? mmmOhzZ`F*$\"3'GV'*e)f6m!)Fft7$$!3LLLL`b`1]F*$\"3v\\hg8mX;8!#?7$$!3#HL LL(G,jYF*$\"3IfUa[%Qa5#Ffu7$$!30nmm'*G7@VF*$\"371q\\^=&>G$Ffu7$$!3XLLL Br9/SF*$\"3w7aJ4Io][Ffu7$$!3C+++qq$fn$F*$\"3@](3$ec*o6(Ffu7$$!3fLLLj<] OLF*$\"3[M*y@WhU.\"!#>7$$!3Q+++I]:)*HF*$\"3w!4kh?-sY\"F`w7$$!3YmmmEQ7] EF*$\"39(yIqi#o_?F`w7$$!3HLLL`xdVBF*$\"3oYpTZwI/FF`w7$$!3I+++![z%)*>F* $\"3OC7hKzY1OF`w7$$!35++++U'>l\"F*$\"3B%o4c1]7q%F`w7$$!3/+++?D.=8F*$\" 31A'zc6*GLfF`w7$$!3SLLLj0z95F*$\"3svsR/IY*=(F`w7$$!3*)ommma1Ul!#=$\"3: Vq[gX3?))F`w7$$!3)=nmm'eW([$Fgy$\"38KSP]EkF5Fgy7$$\"33+'*******G!e'Ffu $\"3qxaU^43*>\"Fgy7$$\"3wDLLL,.6KFgy$\"3o1g.'[$*eM\"Fgy7$$\"3Y.+++H%=m 'Fgy$\"3]c:\\!))[L\\\"Fgy7$$\"33,+++F$y%**Fgy$\"3W$*p!3VaEh\"Fgy7$$\"3 `LLLL=kP8F*$\"3:lbSoMa23*\\yF*$\"3Q-**QC$)e# y\"Fgy7$$\"3_mmmmD5#*>F*$\"391Hq))H6%y\"Fgy7$$\"3Amm;%G#H!3#F*$\"3#f7: ')GuHy\"Fgy7$$\"3%fmm;+#[o@F*$\"3K\"pr'\\o1zl#F*$\"3))R$)GWL))40'>;Fgy7$$\"3W******f0AELF*$\"3!G4f*)**oj\\\"Fgy7$$\"3 M)*****>kThOF*$\"35#[TNAvPN\"Fgy7$$\"3a)****\\.wN#QF*$\"35GUO%f\"Rz7Fg y7$$\"3u)*****\\ct&)RF*$\"3$QLd.)Gv-7Fgy7$$\"3m)****\\D'ylTF*$\"3%>2/# Hp7;6Fgy7$$\"3e)*****fo$eM%F*$\"3'e-%)zRY!H5Fgy7$$\"3RlmmO.i2XF*$\"3;K mPXhN8&*F`w7$$\"3?KLL8QSpYF*$\"3IXdXVT,\\()F`w7$$\"3p*******f!)[,&F*$ \"3;.<;0F2*=(F`w7$$\"3%fmmm\"R$zK&F*$\"3)=NgIntV*eF`w7$$\"3s******zQ=q cF*$\"3*4zX]][*QYF`w7$$\"3mJLLBW@#*fF*$\"3o[?,bzcCOF`w7$$\"3.******H\" H)GjF*$\"3zQ@sR&>\"RFF`w7$$\"3mKLLL:$zl'F*$\"35SlB@*=y.#F`w7$$\"3E**** **\\7Z-qF*$\"3C^X#4p#)3Y\"F`w7$$\"32nmmYRIMtF*$\"37\\Fu())*oO5F`w7$$\" 3?mmm13ltwF*$\"3_V\\\\jmQNrFfu7$$\"33LLL.x=5!)F*$\"3v-Mu7tg:[Ffu7$$\"3 d******f,V>$)F*$\"3?\"*R7/=)*)G$Ffu7$$\"3?LLL8p&Qn)F*$\"3a\"G]Tc2_2#Ff u7$$\"33mmmE/'3**)F*$\"3=`%**HKgcM\"Ffu7$$\"3Q+++!H_)G$*F*$\"3OqeS!pEN H)Fft7$$\"3O+++ION_'*F*$\"3kO&pV?>\"3^Fft7$$\"#5F+Fdt-Fir6&F[sF]sF\\sF \\sF`sFdsFisF`s-%(SCALINGG6#%.UNCONSTRAINEDGFds-%+AXESLABELSG6%Q\"x6\" Q$pdfFigl-%%FONTG6%%*HELVETICAG%(OBLIQUEG\"#=-%&TITLEG6$Q!6\"F[hl-%*AX ESTICKSG6%%(DEFAULTGFihl-F\\hl6%F^hlF_hl\"#7Fis-%*AXESSTYLEG6#%&FRAMEG Fhr-%%VIEWG6$;FbtF]gl;F\\s$\"\"%!\"\"" 1 6 4 1 8 2 2 6 1 3 2 1.000000 45.000000 45.000000 0 0 "Curve 1" "Curve 2" }}}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 31 "The peak of the curve should be" }}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "('1/(sqrt(5)*sqrt(2*Pi))') = evalf(1/(sqrt(5)*sqrt (2*Pi)));" "6#/.*&\"\"\"F&*&-%%sqrtG6#\"\"&F&-F)6#*&\"\"#F&%#PiGF&F&! \"\"-%&evalfG6#*&F&F&*&-F)6#F+F&-F)6#*&F/F&F0F&F&F1" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#/*&\"\"\"F%*&-%%sqrtG6#\"\"&F%-F(6#,$%#PiG\"\"#F%!\" \"$\"+;T7%y\"!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 133 "This shows t hat the modified algorithm is indeed generating the correct Gaussian d istribution with nonunit variance and nonzero mean." }}}}{EXCHG {PARA 0 "" 0 "" {XPPEDIT 19 1 "" "6#%#%?G" }}}}{MARK "1 1 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 1 1 2 33 1 1 }