[F77] Modyfikacja rozkładu Gaussa

[Vigl]

Otrzymuję standardowy rozkład Gaussa, korzystając z generatora liczb pseudolosowych. Jednak dla moich potrzeb ten rozkład musi być bardziej rozlazły (czyt. większa wariancja). Próbowałem to zrobić po ludzku, więc zwyczajnie matematycznie:
N(0,σ)=N(0,1)⋅1/σ⋅exp(x²/2-x²/2σ²),
ale to się sypie.
Kod mojego drivera dla r. standardowego:

Kod: Zaznacz cały

	PROGRAM
	INTEGER seed,N,Nmax,j,i
	PARAMETER(N=500,Nmax=501,NP=10000)
	REAL x,dist(Nmax)

	seed=-12 
	DO i=1,Nmax
	dist(j)=0.0
	END DO

	!rozklad gaussa na podstawie NP punktow
	DO i=1,NP
	j=nint(gasdev(seed)/0.1)	!randomizer gaussowy + rzutowanie
	if ((j.ge.1).and.(j.le.Nmax)) dist(j)=dist(j)+1 !tworzymy krzywa dzwonowa
	END DO

	!wypis pierwszych Nmax wartosci
	DO j=1,Nmax
	x=0.1*float(j)
	dist(j)=dist(j)/NP
	WRITE (*,*) x,dist(j)
	END DO

	END

Skorzystałem tu - jak widać - z funkcji gasdev(seed) (biblioteka num. recipes). Jej opis:

Kod: Zaznacz cały

http://www.nrbook.com/a/bookfpdf/f7-2.pdf

Będę wdzięczny za rady.

[bstq]

\(\displaystyle{ X\sim\mathcal{N}\left(0,1^{2}\right)\Rightarrow\sigma X\sim\mathcal{N}\left(0,\sigma^{2}\right)\text{, bo }Var\left(\sigma X\right)=\sigma^{2}Var(X)=\sigma^{2}}\)
dalej mamy:
\(\displaystyle{ F_{\sigma X}(t)=P\left(\sigma X<t\right)=P\left(X<\frac{t}{\sigma}\right)=F_{X}\left(\frac{t}{\sigma}\right)\Rightarrow f_{\sigma X}(t)=\frac{d}{dt}F_{\sigma X}(t)=\frac{d}{dt}F_{X}\left(\frac{t}{\sigma}\right)=f_{X}\left(\frac{t}{\sigma}\right)\cdot\frac{d}{dt}\left(\frac{t}{\sigma}\right)=\frac{1}{\sigma}f_{X}\left(\frac{t}{\sigma}\right),\;\sigma>0}\)
ostatecznie:
\(\displaystyle{ \frac{1}{\sigma}f_{X}\left(\frac{x}{\sigma}\right)=\frac{1}{\sigma}\cdot\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\frac{x^{2}}{\sigma^{2}}}{2}\right)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right)}\)
a ty chcesz narysować:
\(\displaystyle{ \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right)}\)

[Vigl]

Bajer! Dzięki. Wystarczyło przemnożyć gasdev-a przez wariancje - co jest dość oczywiste. Dla \(\displaystyle{ 10 ^{7}}\) punktów otrzymuję teraz piękną krzywą - dokładnie taką jaką potrzebuję. Nie jestem pewien czy nadal jest unormowana, ale to się jeszcze sprawdzi i najwyżej poprawi.