1

Obliczenia transportu 

promieniowania metodą 

Monte Carlo 

Podstawy

 

2

Podstawy

• Metoda Monte Carlo – dział 

matematyki – wykorzystujące 

stochastyczne rozwiązywanie równań

• Podejście eksperymentalne (“Metoda 

prób statystycznych”)

• Podejście Monte Carlo w celu estymacji 

fizycznie mierzalnych wielkości rozbija 

się na dwa kroki:

1.

Opracowanie eksperymentu numerycznego 

którego wartość oczekiwana  odpowiada 

fizycznej wartości mierzonej wielkości -  

2.

Wykonanie obliczeń – uzyskanie estymacji 

mierzonej wielkości. 

x

xˆ

 

3

Podstawy

• Pierwszy krok – trudność zależy od komplikacji 

systemu fizycznego

• Jeśli system fizyczny jest z natury 

stochastyczny – pierwszy etap jest oczywisty – 

symulacja matematyczna winna być odbiciem 

zwierciadlanym fizycznej sytuacji. Nazywamy 

to symulacją analogową – obliczenia są 

idealnym analogiem sytuacji fizycznej

• Transport promieniowania cząstek neutralnych 

jest takim typem sytuacji stochastycznej

• CO należy zrobić? – odzwierciedlić w modelu 

stochastyczne decyzje podejmowane przez 

przyrodę – tak aby można było je „mierzyć”

 

4

Podstawy

• W procesach które nie są inherentnie 

stochastyczne, model eksperymentu 

numerycznego jest bardziej złożony i 

wymaga:
1. 

Wyprowadzenia równania (np. 

transportu ciepła, Boltzmana) z rozwiązania 

którego można wyciągać wnioski odnośnie 

estymowanych wielkości.

2. Rozwinąć i zastosować metodę Monte 

Carlo do rozwiązania równania. 

 

5

Pierwszy przykład: obliczyć 

• Zadanie: przeprowadzić numeryczną estymację 

 metodą tarczy strzelniczej

• Znamy stosunek powierzchni koła do kwadratu

• Na bazie tego można zaprojektować 

eksperyment dostarczający wartości 

oczekiwanej liczby 

 

kolo}

 

 w

ie

Pr{trafien

4

2

2

2









r

r

 

6

przykład

1. Wylosować punkt wewnątrz kwadratu: 

A. Wybrać losowo liczbę pomiędzy -1 i 1 dla 

położenia we współrzędnej x 

B. Wybrać tak samo dla współrzędnej y

2. Sprawdzić czy trafiono koło i zliczać wyniki:
      trafienie 

s=4 (dlaczego 4?)

pudło:  s=0

3. Przeliczyć eksperyment wiele razy. Końcowa 

estymacja wyliczona ze wzoru:

1

2

2



 y

x

N

s

N

i

i







1

ˆ



 

7

Podstawowe spojrzenie na 

proces MC

Podstawowe spojrzenie na 

proces MC

• Proces Monte Carlo to „czarna skrzynka” – 

wchodzi sekwencja liczb losowych

• Wychodzi sekwencja wartości estymowanych

• Czasami wartości estymowane są przybliżone ale 

przy wydłużeniu ciągu próba staje się 

wystarczająco dobra.

• Potrzebna jest statystyczna ocena próby

 

8

Formuły statystyczne

Formuły statystyczne

• Trzy najważniejsze formuły 

statystyczne:

– Estymowana wartość oczekiwana
– Estymowana wariancja próby
– Estymowana wariancja wartości 

oczekiwanej

– Należy je dobrzy rozróżniać

 

9

Estymator wartości 

oczekiwanej

Estymator wartości 

oczekiwanej

• Najlepszym estymatorem     jest wartość 

średnia

• Zasadnicza cecha – jest nieobciążony – 

spełnia centralne twierdzenie graniczne (z 

N dąży do wartości prawdziwej x)

xˆ

N

x

x

N

i

i







1

ˆ

 

10

Wartości średnie z rozkładów 

ciągłych

Wartości średnie z rozkładów 

ciągłych

• Rozkład ciągły f(x), w przedziale (a,b) (tzn: 

x=a to x=b). 

• Prawdziwa wartość średnia      :

• Założyliśmy że f(x) jest prawdziwym 

rozkładem prawdopodobieństwa tj. 

spełniającym warunki:

x

 

dx

x

f

x

x

b

a





 

1

)

(

b)

(a,

 

przedziale

 

 w

,

0

0









dx

x

f

x

f

 

11

Wartość średnia z rozkładów 

dyskretnych

Wartość średnia z rozkładów 

dyskretnych

• W przypadku rozkładów dyskretnych z M 

wartościami x

i

 o różnych 

prawdopodobieństwach 

• Wartość średnia:

• Warunki na prawdopodobieństwa:

i

p







M

i

i

i

x

p

x

1

1

 

,

0

1









M

i

i

i

p

p

 

12

przykład: problem 

przykład: problem 

• Zadanie  obliczenia     , zawiera 

rozkład dwumianowy (dwie 
odpowiedzi)

• odpowiedź 1 = trafienie w koło:  

• odpowiedź 2 = pudło:  

• Wartość oczekiwana: 



















0

2146

.

0

4

7854

.

0

1

M

i

i

i

x

p

7854

.

0

4

/

4

1

1









p

x

     

2146

.

0

4

/

1

0

2

2











p

x

     

 

13

Estymowana wariancja 

próby

Estymowana wariancja 

próby

• Wariancja – oczekiwany średni błąd 

kwadratowy

• Estymacja wariancji w MC:

 

 





















M

i

i

i

b

a

x

x

p

dx

x

x

x

f

x

1

2

2

2

  

        

  



 

 



























































N

i

N

i

i

N

i

i

i

n

N

x

N

x

N

N

N

x

x

x

S

x

1

2

1

1

2

2

2

2

1

1

ˆ

 

14

Odchylenie standardowe 

próby

Odchylenie standardowe 

próby

• Odchylenie standardowe – 

pierwiastek z wariancji.

• Analogicznie dla wartości 

estymowanych:

• Z własności statystycznych rozkładu 

próby wynika możliwość 
występowania różnic pomiędzy 
indywidualnymi estymacjami.

 

 

n

n

x

S

x

S

2



 

15

przykład: problem 

przykład: problem 

• Wykorzystując poprzednie dane:

• Jest to rozkład bardzo nie-normalny

 









 

642

.

1

697

.

2

697

.

2

2646

.

0

4

7184

.

0

2

2

1

2

2

































x

x

x

p

x

M

i

i

i

        

        

 

16

Estymowana wariancja wartości 

oczekiwanej

Estymowana wariancja wartości 

oczekiwanej

• Znaczenie: miara zaufania do 

otrzymanej wartości średniej z 
pomiaru N próbek 

• Można formułować szereg 

estymatorów a następnie podawać 
je analizie statystycznej

• Przyjmujemy rozwiązania 

standardowe z teorii 
prawdopodobieństwa 

 

17

wariancja wartości 

oczekiwanej

wariancja wartości 

oczekiwanej

• wariancja wartości oczekiwanej 

wykorzystująca wariancję rozkładu 
prawdziwego:

• Musimy jednak wariancję rozkładu 

prawdziwego zastąpić jej estymatorem z próby:

 

 

N

x

x

2

2

ˆ







 

 

 

 

 

 

N

x

S

x

S

x

N

x

S

x

S

x

i

i













ˆ

ˆ

ˆ

ˆ

2

2

2

 

18

przykład: problem 

przykład: problem 

• Nasz przykład dla próby N=10,000 – 

wartość odchylenia standardowego 

średniej wynosi:

• - WIDZIMY podstawową zasadę 

procesu Monte Carlo          :   

zwiększenie dokładności o jeden rząd 

wielkości wymaga zwiększenia 100 

krotnego ilości symulowanych historii 

(co najmniej)

 

 

 

01642

.

0

000

,

10

642

.

1

ˆ

ˆ











N

x

S

x

S

x

i

2

1



N

 

19

Liczby pseudolosowe

• Konieczność wykorzystania liczb losowych 

w MC

• Co to są liczby pseudolosowe
• Co lepsze liczby losowe czy pseudolosowe
• Liniowe generatory przystające 

(kongruentne)

• Wydłużanie okresu generatora – układy 

kombinowane generatorów

• Generatory quasi-losowe

 

20

Ważne reguły

• Nie dysponujemy prawdziwie losowymi 

liczbami (i tak naprawdę nie chcemy ich)

• Liczby losowe z komputera nie są prawdziwie 

losowe (gdyż są powtarzalne) 

• Powtarzalność jest bardzo istotna dla potrzeb 

walidacyjnych. Także jest konieczna w 

przypadku niektórych metod perturbacyjnych 

MC.

• Najbardziej powszechne generatory to liniowe 

generatory przystawalne (kongruentne) (LGP)- 

generują serie powiązanych liczb całkowitych:

m

i

b

ai

i

n

n

n

n

1

1

1

       

m

 

mod

|

)

 

(















 

21

możliwe problemy

•

Periodyczność: ciąg liczb ma skończoną długość, a 
wiec w końcu nastąpią powtórzenia – grozi złamaniem 
podstawowej zasady MC (zwiększania precyzji ze 
zwiększaniem próby) 

•

Błędny wybór liczb startowych a,b, i m może 
prowadzić do powstania korelacji pomiędzy kolejnymi 
liczbami (bez uzasadnienia losowego. Np. po bardzo 
małej liczbie zawsze następuje inna bardzo mała 
liczba) 

•

Można pokazać że pełen okres m będzie pokryty przy 
wyborze:

–

a(mod 4)=1

–

m = potęga 2

–

b - nieparzyste.

•

Inne problemy: wpływ na pokrycie zagadnienia, 
losowość rozkładu cyfr (w zależności od pozycji)

 

 

22

Przykład: m=8, a=5, b=3, 

i0=1

#

i

ai+b mod 

m



0

1

8

0

0

1

0

3

3

.375

2

3

18

2

.25

3

2

13

5

.625

4

5

28

4

.5

5

4

23

7

.875

6

7

38

6

.75

7

6

33

1

.125

 

23

Wydłużanie okresu generatora

• Nawet przy pełnym pokryciu okresu przez generator 

liczb losowych LGP mamy problemy z okresem:

– Okres jest nie większy niż m 
– M jest liczba typu INEGER w systemie 

komputerowym a:

gdzie N jest liczbą bitów w typie INTEGER (zwykle 16 

lub 32) 

• Rozwiązanie problemu: wydłużanie okresu generatora – 

liniowa kombinacja nieskorelowanych generatorów LGP

• Okres jest iloczynem okresów składowych

1

2





N

m





N



















2

1

FRAC

 

24

KENO random number 

generator

DOUBLE PRECISION FUNCTION FLTRN()
INTEGER X, Y, Z, X0, Y0, Z0, XM, YM, ZM
DOUBLE PRECISION V, X1, Y1, Z1
EXTERNAL RANDNUM
COMMON /QRNDOM/ X, Y, Z
SAVE /QRNDOM/
PARAMETER ( MX = 157, MY = 146, MZ = 142 )
PARAMETER ( X0 = 32363, Y0 = 31727, Z0 = 31657 )
PARAMETER ( X1 = X0, Y1 = Y0, Z1 = Z0 )
X      = MOD(MX*X,X0)
Y      = MOD(MY*Y,Y0)
Z      = MOD(MZ*Z,Z0)
V      = X/X1 + Y/Y1 + Z/Z1
FLTRN  = V - INT(V)
RETURN
END

 

25

Generatory quasi-losowe

• Oprócz generatorów pseudolosowych 

wykorzystywane też są generatory quasi-
losowe:

–

Nie są one losowe – trzeba uważnie je stosować

–

Nie są stosowane na ogół w Monte Carlo do 
transportu promieniowania – (dla nas mają 
mniejsze znaczenie) ale mogą być przydatne do 
celów specjalnych. 

•

Bywają bardzo istotne o ogólnych zastosowaniach 
Monte Carlo.

•

-polegają na generowaniu uporządkowanej 
sekwencji liczb pokrywającej przedział i 
teoretycznie mogącej zastąpić liczby pseuro-losowe 
w MC 

 

26

Wybieranie liczb losowych 

z rozkładów

• Rozkład jednorodny
• Rozkład dyskretny
• Rozkład ciągły

– Metoda wyboru bezpośredniego
– Metoda odrzucania
– inne

• Przykład – rozkład Gaussa

 

27

Mapowanie ->x 

wykorzystując p(x)

• Podstawowa technika wykorzystuje 

odwzorowanie jednoznaczne  

y->x

– y= z przedziału (0,1) oraz x z przedziału (a,b)
– Będziemy też wykorzystywać mapowanie odwrotne

 

28

Mapowanie

mamy:
•  (a)=0
•  (b)=1
• Funkcja jest niemalejąca w (a,b)

Zagadnienie sprowadza się do:
• Znalezienia funkcji (x)
• Odwrócenia funkcji (x) aby dostać x(), 

jest to zamiana wylosowanych liczb 

pseudolosowych w liczby podlegające 

zadanemu rozkładowi f(x)

 

29

Mapowanie

• Muszą być spełnione warunki:

f(x)

x

F

dx

x

f

x

dx

x

f

a

x

x

f

dx

x

d

dx

x

f

d

x

x

x

x

a

x

a

 

ta

dystrybuan

 ,

)

(

'

)

'

(

)

(

'

)

'

(

)

(

)

(

)

(

)

(

)

(

}

)

d

,

(

 

 

Pr{

}

)

d

,

(

 

 

Pr{







































 

30

Rozkład ciągły: metoda wyboru 

bezpośredniego

(ogólna procedura)

•

Obliczmy dystrybuantę:

•

 Losujemy liczbę pseudolosową 

•

Przyrównujemy dystrybuantę do wylosowanej 

liczby pseudolosowej

•

Wykorzystujemy formułę odwrotną by 

przekształcić  w x:





x

a

dx

x

f

x

F

'

)

'

(

)

(

)

(x

F





)

(

1







F

x

 

31

Typowy rozkład gęstości 

prawdopodobieństwa

Dystrybuanta

Funkcja odwrotna 

dystrybuanty

 

 

32

Rozkład jednorodny

• wybrać x jednorodnie z (a,b):

• Tworzymy dystrybuantę.

 

a

b

dx

dx

x

f

x

f

x

f

b

a

x

x

f

b

a

b

a

















1

1

1

)

(

~

)

(

~

)

(

,

1

)

(

~

 ,

 

          

a

b

a

x

dx

a

b

dx

x

f

x

F

x

a

x

a

















'

1

'

)

'

(

)

(

 

33

Rozkład jednorodny (cd)

•

Przyrównanie:

•

Odwracanie do x():

•

Przykład: wybrać jednorodnie liczbę z 

przedziału (-1,1):

a

b

a

x









)

(

a

b

a

x









1

2

))

1

(

1

(

1





















 

34

Rozkład dyskretny

• Dokonujemy N wyborów stanu i, 

każdy opisany 
prawdopodobieństwem

• Formujemy dystrybuantę:

i

p

)

(

)

(

)

(

)

(

2

2

1

1

N

N

x

p

x

p

x

p

x

f

































































N

i

N

i

x

x

x

p

x

x

x

p

p

x

x

x

p

x

x

dx

x

f

x

F

1

3

2

2

1

2

1

1

1

,

1

,

,

,

0

'

)

'

(

)

(



 

35

Rozkład dyskretny

• Przyrównanie:

• Odwracanie do x():

)

(x

F





 

































































1

,

,

,

,

1

1

1

1

1

2

1

1

2

1

1

1

ξ

p

x

p

ξ

p

x

p

p

p

x

p

x

F

x

N-

i

i

N

n

i

i

n-

i

i

n

 

if

 

 

if

 

 

if

 

0

 

if

 





 

36

Rozkład ciągły: metoda wyboru 

bezpośredniego

• Formujemy dystrybuantę:

 
• Przyrównanie:

• Odwracanie do x(): 

• Przykład:





x

a

dx

x

f

x

F

'

)

'

(

)

(

)

(x

F





)

(

1







F

x











x

e

x

f

x

0

,

)

(

~

 

37

Rozkład ciągły: metoda 

odrzucania

Ogólna procedura:

1. Znajdź wartość maksymalną rozkładu

:

2. Wybierz losowo

3. Losujemy liczbę pseudolosową 

4. Z

atrzymaj       jeśli 

 

w innym przypadku odrzuć

 

)

,

(

 

gdzie

 ,

 

kazdego

 

dla

 

)

(

~

max

b

a

x

x

x

f

f





 

a

b

x

p

a,b

x

i







1

)

(

 :

 

max

~

)

(

f

x

f

i





i

x

 

38

Przeskalowany rozkład gęstości 

prawdopodobieństwa 

(do metody odrzucania)

Typowy rozkład 

gęstości 

prawdopodobieństwa

sens metody odrzucania - strzelamy 

do tarczy: szare pole odrzucamy, 

białe zostawiamy

 

39

Rozkład ciągły: metoda odrzucania z 

funkcją

 

Podobna zasada ale z wykorzystaniem 
funkcji niejednorodnej dystrybuanty p(x) 
przeskalowanej stałą C 

1. Znajdź taką wartość C że:
2. Wybierz losowo                           z rozkładem 

p(x)

3. Losujemy liczbę pseudolosową 
4. Zatrzymaj       jeśli: 

Uwaga: sprowadza się do zwykłego odrzucania 

gdy p(x) jest stałe oraz: 

 

a,b

x

x

f

x

Cp





 ;

 

)

(

)

(

 

a,b

x

i



)

(

)

(

i

i

x

Cp

x

f





i

x

max

~

)

(

f

a

b

C







 

40

Rozkład ciągły: metoda odrzucania - zadania

A) Użyj metody odrzucania 

dla rozkładu:

wykorzystując
1. Stałe ograniczenie
2. Funkcję ograniczającą:

B) Inne rozkłady

 

1

,

0

,

)

(

~





x

x

x

f

n

2

0

)

(

~







x

e

x

f

x

 ,

2

0

 ,

1

)

(

~









x

kx

x

p

















,

,

)

(

~

2

2

x

e

x

f

x

 

41

Budujemy dwuwymiarowy analogowy 

kod Monte Carlo : 2D MC

• Przegląd wyborów jakie musimy wykonywać 

w czasie symulacji MC transportu cząstki.

1.Początkowa pozycja startowa 
2.Początkowy kierunek lotu 
3.Początkowa energia cząstki  
4.Odległość do następnego zderzenia 
5.Typ zderzenia 
6.Wynik rozpraszania 

• Zagadnienie zliczania strumienia 
• Zaczynamy od kodu homogenicznego 2D

 

42

Decyzja 1: Początkowa pozycja 

startowa

• Zwykle jest to wybór w oparciu o 

wielowymiarowe rozkłady parametrów w 
przestrzeni w określonej objętości. 

• Matematycznie: należy zdefiniować funkcję 

rozkładu cząstek źródłowych w odpowiednim 
systemie współrzędnych (wykorzystując 
istniejące symetrie) a następnie dokonać 
niezależnych wyborów liczb losowych, 
niezależnie dla każdego wymiaru obecnego w 
rozkładzie. 

• Możliwe klasyczne układy współrzędnych: 

kartezjański, cylindryczny, sferyczny,
a także złożone (toroidalne)

 

43

Układ kartezjański

Różniczkowy element objętości: 
 

dz

dy

dx

dV 

 

44

Układ kartezjański (2)

• Chcemy wybrać punkt z płaskim 

rozkładem prawdopodobieństwa 
(każdy element objętości jest 
równie prawdopodobny) 













0

1

0

0

1

0

0

1

0

1

)

(

1

)

(

1

)

(

)

(

)

(

)

(

z

z

z

z

z

f

y

y

y

y

y

f

x

x

x

x

x

f

dz

z

f

dy

y

f

dx

x

f

dV









































      

      

    

 

45

Układ cylindryczny 

(prawoskrętny)

• Różniczkowy element objętości: 
 

dz

dr

rd

dV





 

46

Układ cylindryczny  

(2)

• Płaski rozkład prawdopodobieństwa (każdy 

element objętości jest równie prawdopodobny) 

• Transformacja do układu kartezjańskiego

 





0

1

0

0

1

)

(

)

(

2

1

)

(

)

(

)

(

)

(

z

z

z

z

z

f

r

r

r

r

f

f

dz

z

f

dr

r

f

d

f

dV









































      

      

    

z

z

r

y

r

x











sin

cos

 

47

Układ sferyczny

• Różniczkowy element objętości: 

 



 

dr

rd

d

r

dV







 sin

 

48

Układ sferyczny (2)

• Płaski rozkład prawdopodobieństwa 

(każdy element objętości jest równie 
prawdopodobny) 

• Transformacja do układu kartezjańskiego

 























































2

1

)

(

)

(

2

1

cos

sin

)

(

)

(

)

(

)

(

3

0

2

1

f

r

r

r

r

f

f

d

f

dr

r

f

d

f

dV

      

      

    

















cos

sin

sin

cos

sin

r

z

r

y

r

x

 

49

Wybór ze źródeł złożonych

• W sytuacji złożonego źródła należy 

zastosować strategię mieszania 
prawdopodobieństwa. Możliwa jest wtedy 
kombinacja źródeł z rozkładami o różnej 
gęstości, kształcie, rozmiarze. 

1.

Cząstkowe źródło wybieramy wykorzystując 
całkowitą wydajność źródła 
(czasteczki/sekundę) 

2.

W wylosowanym źródle losujemy punkt 
startowy wykorzystując jego funkcje rozkładu.

 

50

Źródła z niejednorodnym 

rozkładem przestrzennym

• W przypadku niejednorodnym 

cząstkowe rozkłady (danej 
zmiennej przestrzennej) winny być 
przemnożone przez odpowiednie 
funkcje rozkładu. 

• Przykład:

 

51

Decyzja 2: Początkowy 

kierunek lotu

Wybór kierunku bazuje 

na 

prawdopodobieństwie 

wejścia w różniczkowy 

element kąta bryłowego 

      wypełniającego 

pełny kąt bryłowy 

(odpowiednik 

powierzchni sfery o 

jednostkowym 

promieniu)



d











d

d

d

sin

 

52

Decyzja 2: Początkowy kierunek 

lotu 

(2)

• Uwaga: wybór osi Z jest całkowicie dowolny (w 

zależności od potrzeby (istniejących symetrii)

• Na ogół definiujemy                    , wtedy 

element kąta bryłowego wynosi: 

• Znak minus bo       opada gdy         rośnie





 cos











d

d

d





 

53

Decyzja 2: Początkowy kierunek 

lotu 

(3)

• Daje nam to rozkład gęstości 

prawdopodobieństwa: 

• W ogólności w metodzie Monte Carlo 

posługujemy się cosinusami kierunkowymi:









































2

1

)

(

2

1

1

)

(

)

(

)

(

f

f

d

f

d

f

d

      

    



























z

y

x

sin

1

cos

1

2

2

 

54

Decyzja 3: Początkowa 

energia cząstki

    W ogólności początkowa energia cząstki może 

bazować na rozkładzie ciągłym, dyskretnym 

bądź wykorzystywać wielogrupowe spektrum 

źródła. 

• Ciągły: wybór bezpośredni bądź metodą 

odrzucania 

• Dyskretny : (zwykle kwanty ) energia zwykle 

jest sprzężona z wydatkiem 

• Wielogrupowe: Multigroup: źródło w każdej 

grupie jest wycałkowanym rozkładem p-stwa 

na przestrzeni grupy. Zatem wartości 

indywidualnych grup mają charakter dyskretny 

i tak są traktowane.

i

p

 

55

Decyzja 4: Oczekiwana odległość do 

następnego zderzenia

• Dla ośrodka nieskończonego: p-stwa 

zderzenia na dx  w odległości x – p-stwo  

dotarcia do x bez kolizji razy p-stwo 

kolizji na dx

• Zatem funkcja rozkładu gęstości p-stwa: 

(jest już znormalizowana w               )   

 

 





dx

e

dx

x

f

t

x

t











 

 

x

t

t

e

x

f









 







,

0

 

56

Oczekiwana odległość do następnego 

zderzenia (2)

• dystrybuanta : 

• Funkcja odwrotna: 

• W ujęciu długości optycznej:
  mamy:  

 

x

t

e

x

F







1

 





t

x











1

ln

x

t



















ln

   

lub

   

1

ln











 

57

Decyzja 5: Typ zderzenia

• Gdy zderzenie następuje w pewnym punkcie musimy 

wybrać typ reakcji na podstawie wartości 

makroskopowych p. czynnych. 

• Stąd mamy prawdopodobieństwa:

• Wybór typu reakcji wykorzystując w/w p-stwa do 

budowy rozkładu dyskretnego i odpowiedniej 

procedury w tym przypadku   

 

 

 

 

E

E

E

E

c

f

s

t















 

 

 

 

 

 

E

E

p

E

E

p

E

E

p

t

c

c

t

f

f

t

s

s



















     

     

 

58

Decyzja 6: Wynik 

rozpraszania

• Wynik rozpraszania cząstki z energią 

początkową E dany jest rozkładem:

   

gdzie M = materiał, „primy” są związane ze stanem po zderzeniu

• losujemy         wykorzystując:

• losujemy         wykorzystując:























ˆ

ˆ

,

0

0

       

E

E

M
s

0



 





E

d

E

E

f

M
s



















0

0

0

,

E



 





i

M
s

E

E

E

f

0

,











 

59

Wynik rozpraszania (2)

• W pewnych przypadkach rozpraszania występuje 

jednoznaczna zależność pomiędzy kątem 
rozproszenia oraz utratą energii. Wtedy jedna ze 
zmiennych jest zależna od drugiej – wystarczy 
losować tylko jedną.

• W przypadku wielogrupowym zależności kątowe 

w rozpraszaniu pomiędzy grupami są 
reprezentowane przez rozwinięcie Legendre’a w 
funkcji kąta rozpraszania, lub za pomocą 
przedziałów „równo-prawdopodobnych”

 

60

Zliczanie strumienia:

estymator zderzeń (collision estimator)

•

Występują dwa podstawowe sposoby zliczania estymatora 

strumienia (w skończonym obszarze): estymator zderzeń 

(collision estimator) oraz estymator długości trasy 

(tracklength estimator). 

•

Estymator zderzeń bazuje na obliczaniu przyrostu 

strumienia na podstawie występujących (symulowanych) 

zderzeń.

•

Z gęstości reakcji:

•

mamy po każdym zderzeniu 
(infinitezymalny) przyrost strumienia:

•

Można też budować estymatory bazujące na 

poszczególnych reakcjach (konsekwencją – zwiększenie 

błędu – ale wprowadzamy korelacje)

 

   

r

r

r

RR

g

tg

g













 

 

r

r

tg

g









1



 

.

,

,

,

 ,

 

etc

el

s

a

x

r

xg







 

61

Zliczanie strumienia:

e

stymator długości trasy (tracklength 

estimator).

 

•

definicja

•

Zalety: 

–

Możliwa estymacja w próżni

–

Umożliwia określenie co powinno się zdarzyć (np. 

ze zbyt małym p-stwem by pojawiło się przy małym 

rozmiarze próby)

Oba są powszechnie używane o obliczeniach 

reaktorowych, k-eff, - zagadnienia 

krytyczności

c

cx

cx

RR







c

V

c

 

komórce

 

 w

przebyta

 

droga

 

całałkowi





c

 

62

Techniki Redukcji Wariancji

•

Cel: zmniejszyć wariancję estymowanych 
parametrów bez zwiększania rozmiaru próby 

•

Reguły postępowania (oszukiwanie rozkładów)

•

Najczęściej używane techniki:

–

Unikanie absorpcji (ze zmiana wag w miejscu 
absorpcji)

–

Rozszczepianie cząstek (nie mylić z 
rozszczepianiem atomów) 

–

Wymuszanie zderzeń 

–

Rosyjska ruletka

–

Transformacja eksponencjalna 

–

Modyfikacja (biasing) źródła (pozycji, energii lub 
kierunku) 

 

63

Cel stosowania

•

Aby uzyskać lepszy program Monte Carlo – 

lepszy to znaczy generujący wyniki z mniejszą 

wariancją

•

Mniejszą wariancja w ujęciu czarnej skrzynki:

 

oznacza że x-sy są  do siebie najbardziej 

zbliżone (jak to możliwe)

•

Skuteczny dla lokalnych wielkości – gorzej w 

przypadku globalnych parametrów

 

64

Reguły postępowania (oszukiwanie 

rozkładów)

•

Mapowanie prowadzimy wykorzystując (zamiast 

naturalnego) zmodyfikowany rozkład 

prawdopodobieństwa – do postaci przez nas pożądanej 

tak by zagęścić wybór w obszarach naszego 

zainteresowania. 

•

Zgodność z naturą zachowujemy zmieniając wagę 

cząstek  

•

Waga cząstki odpowiada części cząstki jaka w 

naturalnym rozkładzie pojawiłaby się w danej sytuacji.

•

Jeśli naturalny rozkład daje gęstość p-stwa:           a 

my chcemy użyć rozkładu o gęstości:               to 

musimy zastosować korektę wagi:

 

•

Konieczny wymóg: wszystkie możliwe punkty x 

dostępne w naturalnym rozkładzie muszą być 

dostępne w nowym – zmodyfikowanym rozkładzie.  

 

 

x

p

x

p

w

corr

~



 

x

p

 

x

p

~

 

65

Przykład: wykorzystanie rozkładu 

płaskiego

• Z lenistwa wykorzystujemy do mapowania 

najprostszy rozkład płaski

• Nie prowadzi do zniekształcenia wyników 

jeśli skorygujemy wagi następująco:

• Chociaż wynik nie będzie zniekształcony, 

postępowanie takie na ogół nie jest 

rozsądne. (potrzebna jest odpowiednia 

motywacja do właściwej zmiany rozkładu)

 

 

 





 



a

b

x

p

a

b

x

p

x

p

x

p

w

corr











1

~

 

66

Unikanie absorpcji 

• Najczęściej używany technika redukcji 

wariancji. W wielu kodach transportowych jest 
ona wymuszona (nie może być wyłączona 
przez użytkownika) 

• Pociąga zmianę wag w miejscu absorpcji

 

• Opis ogólny: Podstawowa zasada polega na 

uniemożliwieniu absorpcji cząsteczki; wtedy 
cząsteczka będzie żyć dłużej i mieć więcej 
szans na bycie zliczaną w innych procesach. 

• Losowy wybór wymagający dopasowania:

- typ zderzenia

 

67

Unikanie absorpcji (2)

• Zarys matematyczny: 

–

Załóżmy że są możliwe dwa stany wyjściowe 

zderzenia: rozpraszanie i absorpcja. W stanie 

naturalnym mamy dystrybucje (tego co powinno 

być):  

–

rozpraszania z odpowiednim rozkładem p.stwa: 

–

absorpcji z odpowiednim rozkładem p.stwa:

•

Nasze decyzje (mapowanie) podejmujemy jednak 

tak by zawsze wypadło rozpraszanie (ze 100% 

p.stwem) 

•

Wymagana korekcja wagi cząstki rozproszonej:

t

s

p







1

 

 

t

s

t

s

corr

x

p

x

p

w















00

.

1

~

t

a

p







2

 

68

Unikanie absorpcji (3)

• Uwaga: zmieniliśmy naturalnie możliwy wybór 

(absorpcja) NIEMOŻLIWYM. Jest to możliwe 

tylko dlatego że cząstka w wypadku absorpcji i 

tak przestaje żyć (nie wnosi żadnych nowych 

wkładów do odpowiedzi systemu) 

• Procedura ta oprócz zmniejszenia wariancji (na 

historię) gwarantuje zwiększenie czasu obliczeń 

na skutek wydłużenia każdej historii.

• Czasami to się nie musi opłacać 
• Inne spojrzenie na tę technikę: Dzielimy 

cząstkę na dwie części: ułamek który się 

rozprasza i ułamek który jest absorbowany. 

Kontynuujemy śledzenie tylko pierwszej 

podczas gdy druga umiera.

 

69

Rozszczepianie cząstek (nie 

atomów)

• Stanowi podstawę innych technik redukcji 

wariancji

• Nie wypełnia wzorca modyfikacji 

dystrybucji p.stwa. Polega na unikaniu 
dyskretnych decyzji w zamian za śledzenie 
wszystkich możliwych opcji decyzyjnych

• Ponieważ probabilistyczne decyzje są 

unikane, wpływ tych decyzji na wariancje 
też znika.

 

70

Rozszczepianie cząstek

 (2)

• Dwie typowe sytuacje w transporcie cząstek z 

wykorzystaniem rozszczepienia (historii)

1. W przypadku zderzenia, historia jednej cząstki 

przebiega dalej niezmienionym śladem (tylko z 

mniejszą wagą) aż do zakończenia życia a wtedy 

wracamy do punktu rozszczepienia by prowadzić 

historię cząstki rozproszonej aż do kończ jej życia. 

2. Cząstka jest rozszczepiana na dwie lub więcej w 

sytuacji wejścia w obszar o większej ważności (tam 

gdzie chcemy zmniejszyć wariancje - na ogół by ją 

wyrównać do poziomu średniego) 

• Matematycznie idea ta polega na śledzeniu 

wszystkich możliwych opcji dyskretnych rozkładów 

p.stwa z odpowiednim ważeniem ich historii 

proporcjonalnie do wartości dyskretnych p.stw.

 

71

Rozszczepianie cząstek

 (3)

• Przykład: zakładamy historię cząstki której 

obecna waga wynosi 0.5 napotyka na dyskretną 

decyzję z trzema możliwymi stanami:

– stan 1 z 60% prawdopodobieństwem;
– stan 2 z 30% prawdopodobieństwem;
– stan 3 z 10% prawdopodobieństwem.

• Unikamy wyboru, a w zamian śledzimy każdą 

opcję.

– Kontynuujemy historię w stanie 1 z wagą 0.3 

zapisując w banku (by kontynuować później) historię 

w stanie 2 z wagą 0.15 oraz historię 3 z wagą 0.05.

• Po skończeniu każdej historii wracamy do 

banku by podjąc pozostające tam cząstki aż 

opróżnimy go całkowicie.

 

72

Wymuszanie zderzeń

• Ogólny opis: Wymuszamy zderzenie w 

szczególnym obszarze drogi cząstki. 
Wykonamy to w szczególny sposób aby 
zapobiec ucieczce cząstki (bez zderzenia) 
z tego obszaru. 

• Utrzymuje cząstki dłużej przy życiu, zatem 

zwiększa ich kontrybucje statystyczne do 
wszystkich zliczeń. 

•

Losowy wybór wymagający dopasowania:
- dystans do najbliższego zderzenia

 

73

Wymuszanie zderzeń (2)

• Zarys matematyczny: zasadniczo jest to 

technika rozszczepienia cząstki na część 

ulegającą

 zderzeniu w określonym obszarze 

oraz  na część oraz 

nie ulegającą

 zderzeniu 

w nim. 

• Załóżmy że odległość do granicy wynosi 0 

(średnich dróg swobodnych) 

• Naturalny rozkład p.stwa ma postać:

–

P.stwo ucieczki:

–

P.stwo uniknięcia ucieczki:

• Nasza decyzja to uniknięcie ucieczki ze 100% 

p.stwem

0





e

p

escape

0

1









 e

p

escape

non

 

74

Wymuszanie zderzeń 

(2)

• Korekta wagi wynosi: 

• Rozszczepieniowa natura tej techniki powoduje 

że musimy wziąć pod uwagę że część cząstki 

ucieka poza dany obszar. Ta część powinna być 

uwzględniona przy zliczaniu ucieczki poza 

najbliższa granicę, w ilości:           

     , gdzie 

w jest wagą cząstki przed przejściem. 

• Także tutaj efektywność może się nie poprawić 

(tzn. straty z wydłużenia czasu obliczeń 

przewyższą zysk z mniejszej wariancji)

–

Faktycznie ta technika nie jest zbyt często używana. 

Aczkolwiek czasami występuje konieczność jej użycia 

(np. symulacja odpowiedzi małego detektora)

 

 

0

0

1

1

1

~



















e

e

x

p

x

p

w

i

i

corr

0





we

 

75

Rosyjska ruletka 

• Podobnie do rozszczepienia 

matematyczne narzędzie wykorzystywane 
w technikach zmniejszania wariancji

• Idea: kombinacja kilku cząstek w jedną 

przy pomocy selektywnej eliminacji (jak 
w rosyjskiej ruletce) 

• Matematycznie oznacza to zabijanie 

cząstek z dużym p.stwem 1-p (90-99%) 
gdzie p oznacza szansę przeżycia.  

• Jeśli cząstka przeżyje, jej waga jest 

zwiększona o czynnik 1/p

 

76

Rosyjska ruletka (2)

• Narzędzie to jest wykorzystywane w przypadku 

stosowania innych technik redukcji wariancji do 

zakańczania historii cząstek których waga  jest na 

tyle mała że tracą znaczenie.

• Narzędzie to jest konieczne w przypadku 

stosowania jednocześnie techniki unikania absorpcji 

oraz wymuszania zderzeń (na całym obszarze) gdyż 

jest to wtedy jedyna sposób na zakończenie historii 

cząstek bo użyte techniki wyeliminowały naturalne 

drogi zakończenia historii. 

• W praktyce jest wykorzystywane zawsze gdy waga 

cząstki spadnie poniżej dolnego progu odcięcia 

wagi..  

–

Technika ta nie redukuje wariancji wprost, tylko pozwala na 

oszczędzenie czasu obliczeniowego poprzez przerywanie 

historii o nikłym znaczeniu – słowem pozwala zwiększyć 

efektywność.

 

77

Transformacja 

eksponencjalna

• Ogólny opis: idea ta polega na ułatwieniu 

cząstce podróży w pożądanym przez nas 

kierunku, poprzez rozrzedzenie materiału 

w tym kierunku oraz zagęszczenie w 

kierunkach niepożądanych.  

– Nie powoduje że pożądany kierunek będzie 

częściej wybierany a jedynie ułatwi penetrację 

materiału w tym kierunku. 

• Losowy wybór wymagający dopasowania:

- dystans do najbliższego zderzenia

 

78

Transformacja eksponencjalna 

(2)

• Zarys matematyczny: zmieniamy 

całkowity przekrój czynny oraz 

uzależniamy go od kierunku lotu:

gdzie p jest ogólnym parametrem 

definiowanym przez użytkownika (0<p<1) 

oraz 

d

 jest wyróżnionym kierunkiem

 

• Ograniczenia zmienności p do przedziału 

(0<p<1) powoduje że zmieniony przekrój 

czynny będzie większy od zera i mniejszy 

od dwukrotności prawdziwego przekroju 

czynnego.

  



d

t

t

p 















ˆ

ˆ

1

ˆ

*

 

79

Transformacja eksponencjalna 

(3)

• Rozważmy przypadek cząstki 

lecącej w kierunku zewnętrznej 

granicy. Mamy dwa przypadki: 

– Cząstka osiągnie granicę
– Cząstka się rozproszy 

• Naturalny rozkład p-stw ma postać: 

– P.stwo ucieczki:  

– P.stwo zderzenia w odległości s:

0





e

p

escape

s

t

t

e

s

p









)

(

 

80

Transformacja eksponencjalna

 

(4)

• Wykorzystywany rozkład ma postać:

– P.stwo ucieczki:  
– P.stwo zderzenia w odległości s:

• Korekta wagi: 

– W przypadku ucieczki:  

– W przypadku zderzenia w odległości s: 





d

p

escape

e

p













ˆ

ˆ

1

0



 









s

p

d

t

d

t

e

p

s

p

























ˆ

ˆ

1

ˆ

ˆ

1





d

d

p

p

corr

e

e

e

w

























ˆ

ˆ

ˆ

ˆ

1

0

0

0



















d

p

p

d

t

t

corr

p

e

e

p

e

w

d

d













































ˆ

ˆ

1

ˆ

ˆ

1

ˆ

ˆ

ˆ

ˆ

1

0

0

0







 

81

Modyfikacja źródła

• Jest to najbardziej ogólna technika 

redukcji wariancji. Zasada jest 
prosta aczkolwiek nie dająca 
jasnych kryteriów postępowania: 
zamiast wykorzystywać naturalny 
rozkład cząstek źródłowych 
wykorzystujemy jakiś inny rozkład o 
którym myślimy że będzie lepszy.

• Losowy wybór wymagający dopasowania:

- pozycja cząstki, energia i kierunek lotu

 

82

Modyfikacja źródła (2)

• Zarys matematyczny:

–

Zasadniczo nie ma ścisłych wskazówek jak 
postępować – musimy wykorzystywać ogólne 
reguły oszukiwania rozkładów. 

–

Jeśli fizyka dyktuje rozkład p.stwa:

a my chcemy użyć:

możemy to zrobić stosując korektę wagi:

 

                                                                

                                                                              

            

 

 

x

p

x

p

w

corr

~



 

x

p

 

x

p

~

 

83

Wybór zmodyfikowanego rozkładu

• W ogólności chcemy zmodyfikować 

naturalny rozkład tak by faworyzować 
wybory bardziej ważne. 

• Pytanie „Co jest bardziej ważne?”
• Odpowiedź:

WAŻNOŚĆ= SPODZIEWANA KONTRYBUCJA
jako kontrybucję rozumiemy zliczenia estymatora 
wybranej wielkości fizycznej – np. sygnału 
detektora

•  Dlatego nasze zadanie polega na takiej 

modyfikacji dystrybucji aby kontrybucje do 
naszych zliczeń były największe.

 

84

Wybór zmodyfikowanego rozkładu 

(2)

• Jeśli znamy funkcję ważności 

(importance) możliwych wyborów I(x), 
wtedy optymalną teoretycznie 
dystrybucją będzie:

• Możliwe podejścia można podzielić na 

kategorie: 

–

(1) heurystyczną (bazującą na intuicji) 

–

(2) eksperymentalną 

–

(3) wykorzystującą sprzężony strumień

 

   

x

I

x

p

x

p

~

~

 

85

Podejście heurystyczne

•

Modyfikujemy naturalny rozkład tak by 

faworyzować wybory cząstek o których myślimy 

że są bardziej ważne dla naszego zagadnienia.

–

Przykład 1:  położenie cząstki – jeśli interesuje nas 

ucieczka  przez prawą stronę – z preferencją 

wybieramy położenia bliższe prawej strony 

–

Przykład 2:  kierunek lotu -  jeśli interesuje nas 

ucieczka  przez lewą stronę – z preferencją wybieramy 

kierunki w lewo 

–

Przykład 3:  energia cząstki- jeśli interesuje nas 

głęboka penetracja – preferujemy wybory energii dla 

których całkowity przekrój czynny jest mały

•

Bazujemy tutaj na naszej intuicji – ocenie jak 

bardzo należy faworyzować ważniejsze cząstki. 

Jest to metoda prób i błędów – przydatne bywa 

doświadczenie.

 

86

Podejście eksperymentalne 

•

Podejście to bazuje na wykonaniu kilku krótkich 
przebiegów testowych w celu uzyskania wartości 
względnych funkcji ważności dla różnych 
możliwych wyborów parametrów źródłowych. 

•

Każdy test powinien posiadać ograniczony 
obszar (domenę) zmienności poszczególnych 
parametrów źródłowych, a po wykonaniu 
obliczeń wykorzystujemy uzyskaną funkcji 
ważności I(x) do utworzenia skorygowanego 
rozkładu:

•

Możemy też zaangażować niepewność obliczonej 
funkcji:

 

   

x

I

x

p

x

p

~

~

 

   

 

 

 

x

x

I

x

p

x

I

x

p

x

p

2

2

~

~

~





 

87

Optymalna modyfikacja -  

wyprowadzenie 

wzorów

•

Funkcja ważności I(P) jest wykorzystywana do modyfikacji 

źródła w następujący sposób:

gdzie:   
optymalna (nieznormalizowana dystrybucja 

do 

wykorzystania w symulacji)
prawdziwa (naturalna) dystrybucja (w 

poprzednich 

wzorach – p(x) )

Oczekiwana kontrybucja przy wyborze naturalnym - P

 

)

(

)

(

~

~

P

I

P

S

P

p







 

 )

(

 )

(

 

~

P

I

P

S

P

p







 

88

Optymalna modyfikacja 2

Wyprowadzenie:

•

Użycie           zamiast             daje nam oczekiwaną 
wariancję w zmodyfikowanej symulacji optymalnej w 
postaci:

=>

gdzie          jest prawdziwą kontrybucją 
niezmodyfikowaną,

•

Nałożenie warunku optimum (najmniejszej wariancji)

prowadzi do następującego wyniku:

 

 



 

 

 

P

d

I

P

p

P

I

P

S

P

p

I

P

d

I

P

I

w

P

p

I

P

i

P

corr

i

















2

2

2

2

~

)

(

)

(

~

)

(

~





























 

P

p



~

)

(P

I



)

(P

S



 

0

~

~

2









p

p

 

89

Optymalna modyfikacja - 3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I

P

I

P

S

P

p

P

d

I

P

p

P

I

P

S

I

P

p

P

I

P

S

P

d

P

p

P

I

P

S

I

P

p

P

I

P

S

I

P

p

P

I

P

S

P

d

P

p

P

I

P

S

I

P

p

P

I

P

S

P

p

I

P

p

P

I

P

S

P

d

I

P

p

P

I

P

S

P

p

p

P

d

I

P

p

P

I

P

S

P

p

p

p

p

P

P

P

P

P

)

(

)

(

~

0

~

)

(

)

(

~

)

(

)

(

~

)

(

)

(

2

~

)

(

)

(

~

)

(

)

(

~

)

(

)

(

~

)

(

)

(

~

2

~

)

(

)

(

~

)

(

)

(

~

~

~

)

(

)

(

~

~

~

~

2

2

2

2

2

2











































































































































































































































































  

if

   

 

90

Optymalna modyfikacja (4)

•

Problem pojawić się może gdy nasz wybór       

nie gwarantuje kontrybucji do naszego zliczania

•

Natomiast gdy wybór       wywołuje  kontrybucje 

z niepewnością             , wtedy optymalna 

modyfikacja ma postać:

•

W przypadku poisson’owskim (głęboka 

penetracja analogowa) mamy:

Wtedy optimum jest

(oznacza to zanik optymalnej modyfikacji)

)

(

)

(

)

(

P

I

P

I

P













P



)

(P

I



P



)

(P

I



)

(P





 

   

2

2

)

(

)

(

)

(

~

~

P

P

I

P

S

P

p













 

)

(

)

(

~

~

P

I

P

S

P

p







 

91

Ważenie po celach 

(MCNP)

•

Rules:

1.

Gdy cząstka przechodzi z obszaru A do B o 

większej ważności wtedy rozszczepiamy w 

stosunku:

2.

Gdy cząstka przechodzi z obszaru A do B o 

mniejszej ważności wtedy gramy rosyjską 
ruletkę z szansą przeżycia:

     

old

new

I

I

 

old

new

I

I

 

 

92

Okna wagowe w MCNP

Reguły:

Przy przejściu cząstki z obszaru A do nowego B 

gdzie następuje zmiana wagi dostosowujemy 
wagę tak by mieściła się w wymaganym 
przedziale (oknie wagowym) 

Stosujemy RR lub rozszczepienie wykorzystując:

desired

current

w

w

 

high

desired

low

w

w

w



