Fortran - Add Noise

Author

Add Noise

Filippos Louis

2004-03-27, 12:16 am

Hi all
Does anyone know of a Fortran routine to add noise to data on a vector?
(Gaussian or Random)!
Thank you in advance

FLouis

David Frank

2004-03-27, 12:16 am

"Filippos Louis" <flouis@geol.uoa.gr> wrote in message
news:c3c3cm$33i$1@nic.grnet.gr...
> Hi all
> Does anyone know of a Fortran routine to add noise to data on a
> vector?
> (Gaussian or Random)!
> Thank you in advance
>
> FLouis
>
>
>

No routine needed...

call random_number(noise)
noisy_vector = vector + (noise - 0.5)*scale

tholen@antispam.ham

2004-03-27, 12:16 am

David Frank writes:

> Filippos Louis wrote:

[color=darkred]
> No routine needed...
>
> call random_number(noise)
> noisy_vector = vector + (noise - 0.5)*scale

That won't add Gaussian noise.

Richard Maine

2004-03-27, 12:16 am

tholen@antispam.ham writes:

> David Frank writes:
>
>
>
>
> That won't add Gaussian noise.

True. Though it will add random noise of some kind. I find it a
little hard to decipher "Gaussian or Random". Technically this
does fit the "or random" part, but I don't know what was really
meant there.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Karagiorgos

2004-03-27, 12:16 am

Thank you all!!

"Richard Maine" <nospam@see.signature> wrote in message
news:m1oequhz5y.fsf@macfortran.local...
> tholen@antispam.ham writes:
>
>
> True. Though it will add random noise of some kind. I find it a
> little hard to decipher "Gaussian or Random". Technically this
> does fit the "or random" part, but I don't know what was really
> meant there.
>
> --
> Richard Maine | Good judgment comes from
experience;
> email: my first.last at org.domain | experience comes from bad judgment.
> org: nasa, domain: gov | -- Mark Twain

Joost VandeVondele

2004-03-27, 12:17 am

tholen@antispam.ham wrote in message news:<dJg6c.7907$h85.3597@twister.socal.rr.com>...
> David Frank writes:
>
>
>
>
> That won't add Gaussian noise.
noisy_vector=vector
DO i=1,N
call random_number(noise)
noisy_vector=noisy_vector+(noise -0.5)*scale
ENDDO

Rich Townsend

2004-03-27, 12:17 am

Joost VandeVondele wrote:
> tholen@antispam.ham wrote in message news:<dJg6c.7907$h85.3597@twister.socal.rr.com>...
>
>
> noisy_vector=vector
> DO i=1,N
> call random_number(noise)
> noisy_vector=noisy_vector+(noise -0.5)*scale
> ENDDO

Still ain't Gaussian, if that was what you were going for. The
RANDOM_NUMBER() intrinsic returns pseudorandom numbers uniformly
distributed between 0 and 1. To convert this distribution into a
Guassian one, with a specific standard deviation, use something along
the lines of the gasdev() routine in Numerical Recipes:

http://lib-www.lanl.gov/numerical/bookfpdf/f7-2.pdf

cheers,

Rich

tholen@antispam.ham

2004-03-27, 12:17 am

Joost VandeVondele writes:

[color=darkred]
[color=darkred]
[color=darkred]
[color=darkred]
[color=darkred]
> noisy_vector=vector
> DO i=1,N
> call random_number(noise)
> noisy_vector=noisy_vector+(noise -0.5)*scale
> ENDDO

Inefficient, and it's not immediately obvious what the
distribution would look like, though it should at least
be symmetric.

Richard Maine

2004-03-27, 12:17 am

Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:

> Joost VandeVondele wrote:
[color=darkred]
>
> Still ain't Gaussian...

But by the central limit theorem, it approaches Gaussian as N
approaches infinity (assuming you normalize scale appropriately).
And in this case, numbers on the order of 12 or so are a decent
approximation to infinity. I've seen code that uses this kind
of method to approximate Gaussians (and I'm thinking I recall
N=12 being used; might have even been N=6).

Perhaps not the world's most efficient method, but it is hard
to get simpler, and for some applications it is fine.

Hey, if you needed to be really picky about things, you'd be using a
3rd party generator instead of whatever the compiler vendor happened
to throw in for random_number.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Rich Townsend

2004-03-27, 12:17 am

Richard Maine wrote:
> Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:
>
>
>
>
>
>
> But by the central limit theorem, it approaches Gaussian as N
> approaches infinity (assuming you normalize scale appropriately).

You're right, I don't think I read or understood the code properly. As
you point out, though, there is a gotcha in normalizing scale properly;
if scale is kept independent of N, then the standard deviation of the
resulting Gaussian distribution will grow as SQRT(N).

This is, of course, why a random walk can get a particle from A to B. It
is also why the so-called 'law of averages' -- that a heads is highly
probable after tossing 6 consecutive tails -- is total bunk. And also
why the phenomenon of 'gamblers ruin' will mean that Casinos will always
win (assuming they are richer than the punters), even when the odds are
completely fair.

cheers,

Rich

Gordon Sande

2004-03-27, 12:17 am

In article <m14qslj3f8.fsf@macfortran.local>,
Richard Maine <nospam@see.signature> wrote:

>Subject: Re: Add Noise
>From: Richard Maine <nospam@see.signature>
>Organization: NASA Dryden
>Date: 18 Mar 2004 12:39:07 -0800
>Newsgroups: comp.lang.fortran
>
>Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:
>
>
>
>But by the central limit theorem, it approaches Gaussian as N
>approaches infinity (assuming you normalize scale appropriately).
>And in this case, numbers on the order of 12 or so are a decent
>approximation to infinity. I've seen code that uses this kind
>of method to approximate Gaussians (and I'm thinking I recall
>N=12 being used; might have even been N=6).

The variance of the common uniform is 1/12 so the sum of 12 of them
has the right variance of 1.0. It is a workable approximation at 1:00
in the morning and you need a quick and dirty fix by 1:10, providing
you comment it as being not much better than nothing.

>Perhaps not the world's most efficient method, but it is hard
>to get simpler, and for some applications it is fine.
>
>Hey, if you needed to be really picky about things, you'd be using a
>3rd party generator instead of whatever the compiler vendor happened
>to throw in for random_number.
>
>--
>Richard Maine | Good judgment comes from experience;
>email: my first.last at org.domain | experience comes from bad judgment.
>org: nasa, domain: gov | -- Mark Twain

beliavsky@aol.com

2004-03-27, 12:17 am

"Filippos Louis" <flouis@geol.uoa.gr> wrote in message news:<c3c3cm$33i$1@nic.grnet.gr>...
> Hi all
> Does anyone know of a Fortran routine to add noise to data on a vector?
> (Gaussian or Random)!
> Thank you in advance
>
> FLouis

Here is a Fortran 90/95 function with driver to simulate Gaussian
variates using the Box-Muller algoritm.

module box_muller
implicit none
public :: normal_bm
contains
function normal_bm() result(zz)
! generate a random normal deviate using the Box-Muller method
real :: zz
real :: dist,xx(2)
real, parameter :: two = 2.0, minus_two = -2.0
do
call random_number(xx)
xx = two*xx - 1.0
dist = xx(1)**2 + xx(2)**2
if (dist < 1.0) exit
end do
zz = xx(1)*sqrt(minus_two*log(dist)/dist)
end function normal_bm
end module box_muller

program xnormal_bm
use box_muller, only: normal_bm
implicit none
integer :: i
integer, parameter :: n = 1000000
real :: xx(n)
call random_seed()
do i=1,n
xx(i) = normal_bm()
end do
write (*,"(' moments:',100f8.4)") &
(/ sum(xx),sum(xx**2),sum(xx**3),sum(xx**4)
/) / n
end program xnormal_bm

With Compaq Visual Fortran 6.6 it produce the result

moments: 0.0010 0.9970 0.0037 2.9859

which looks good. Can anyone recommend a Fortran code to test the
quality of a vector of variates that are supposed to be IID normal?

Some sources for Fortran codes for random number generation are listed
at http://www.dmoz.org/Computers/Progr...d_Econometrics/
..

Sponsored Links