Rev 4161 | Go to most recent revision | Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
20 | reyssat | 1 | !if $wims_read_parm!=slib_header |
2 | !goto proc |
||
3 | !endif |
||
4 | |||
5 | slib_title=Generation of Gaussian random data |
||
6 | slib_parms=3\ |
||
7 | ,m\ |
||
8 | 0,mu\ |
||
9 | 1,sigma |
||
10 | |||
11 | slib_author=Adeline Grelot and Bernadette PERRIN-RIOU |
||
12 | slib_out=random data of m reals following a normal law of parameters mu and sigma. |
||
13 | slib_comment=Box Muller method |
||
14 | slib_example=10,0,2 |
||
15 | !exit |
||
16 | |||
17 | :proc |
||
18 | !distribute item $wims_read_parm into slib_m,slib_mu,slib_sigma |
||
19 | !default slib_mu=0 |
||
20 | !default slib_sigma=1 |
||
21 | |||
22 | slib_m2= $[floor($slib_m/2)] |
||
23 | slib_U1=!random 10^-20,1 repeat $[$slib_m2+1] |
||
24 | slib_U2=!random 0,1 repeat $[$slib_m2+1] |
||
25 | |||
26 | slib_w=!values sqrt(-2*log(x)) for x in $slib_U1 |
||
27 | slib_w1=!values cos(2*Pi*x) for x in $slib_U2 |
||
28 | slib_w2=!values sin(2*Pi*x) for x in $slib_U2 |
||
29 | |||
30 | slib_V=$empty |
||
31 | !for slib_i=1 to $[$slib_m2 +1] |
||
32 | slib_v=!item $slib_i of $slib_w |
||
33 | slib_v1=!item $slib_i of $slib_w1 |
||
34 | slib_v2=!item $slib_i of $slib_w2 |
||
35 | slib_v1=$[$slib_sigma*$slib_v1*$slib_v+$slib_mu] |
||
36 | slib_v2=$[$slib_sigma*$slib_v2*$slib_v+$slib_mu] |
||
37 | |||
38 | !if $slib_i<=$slib_m2 |
||
39 | slib_V=!append item $slib_v1 to $slib_V |
||
40 | slib_V=!append item $slib_v2 to $slib_V |
||
41 | !endif |
||
42 | !if $slib_i>$slib_m2 and $[2*$slib_m2]!=$slib_m |
||
43 | slib_V=!append item $slib_v1 to $slib_V |
||
44 | !endif |
||
45 | !next slib_i |
||
46 | |||
47 | !if $slib_m >1 |
||
48 | slib_V=!shuffle $slib_V |
||
49 | !endif |
||
50 | |||
51 | slib_out=!trim $slib_V |