Artificial Data / Random - Department of Physics



Artificial data

[pic](1.1)

[pic](1.2)

[pic](1.3)

for\PEAK.FOR

IMPLICIT REAL*8 (A-H,O-Z)

PARAMETER (N=5000)

DATA CJ/1D12/,EJ/511D0/WJ/10D0/

OPEN(1,FILE='PEAK.OUT')

H=600D0/N

DO I=1,N

E=(I-.5D0)*H

ARG=(E-EJ)/WJ

PEAK=CJ*EXP(-ARG*ARG)

WRITE(1,'(2G15.6)')E,PEAK

ENDDO

END

[pic]

Changing the d12 to d6 shows a more typical experiment. In general the constants a,b,cj,wj,Ej,dj, and dj+1 need to be determined along with their standard deviations.

[pic]

The standard deviation is the difference between the value determined in a single experiment, and that that would be determined by repeating the experiment an infinite number of times. It is not the difference between the correct answer and that given by the experiment.

The data observed is not the number given by the two codes above, but differs by a “random amount”. In many cases it differs by a Gaussian centered at the data value so that

[pic](1.4)

D:\public_html\class2K\Progdet\sorting\Random.htm

Psuedo Random numbers

Random numbers are not wanted. What is desired is a set of completely reproducible numbers that are the same on all computers and in both C and Fortran. In addition, it must be true that the numbers give all averages the same as would be given by actual random numbers. In practice the numbers should be re-initialized with the numbers needed to recreate the desired set written out as

17894561 3678991

18345899 234567899

...

In this way, when the code bombs after 23 hours or running, the situation just before the bomb can be recreated by simply entering these numbers.

With this many random numbers, the "randomness" of the set can be a problem for some calculations[1]. The method used here, first suggested by MacLaren and Marsaglia[2], is to insert numbers into a table with a multiplicative congruential method with one multiplier and seed and to extract them from the table with a multiplicative congruential method with a different multiplier and seed. random.for contains the common

COMMON/RAN/IX,IY,ITAB(128)

RANDOM.CPP RANDOM.H

The file random.h contains

struct {int ix,iy,itab[128];} ran;

int rseed(int ix,int iy);

double rndmf();

The seeds are ix and iy. The random number set is initialized by calling rseed(ix,iy). This call causes 128 numbers to be fed into the table itab[128] with one random number generator. The call x=rndmf() gives a random number base on ix,iy and itab, between 0 and 1. It selects this by advancing ix and iy, pulling a value from a randomly selected table entry and replacing that with a new value.

This can be re-initialized by printing out ix and iy and then calling rseed. This wastes 128 random numbers, but gives a completely determined new set of random numbers.

Simple Monte Carlo Integration

[pic]

Change variables such that x=a+(b-a)*x'. dx = (b-a)dx'

[pic]

now select x positions uniformly and randomly between 0 and 1 and average them so that

[pic]

"obvious" error analysis

The expected error in [pic] can be estimated by considering each random selection as an independent estimator of [pic] and then averaging the squares of the differences to yield

[pic]

as the rms difference between [pic] and any single estimator. The expectation value of the sum is then N times (N so that the deviation in the sum is

[pic]

and the expected error in the sum is

[pic]

while the expected error in the average is

[pic]

so that

[pic]

cpp\trandom.c cpp\random.c cpp\random.h

#include

#include

#include "random.h"

void main(void)

{int n,ir,ix=37,iy=72;

double sum,sum2,err,xtemp;

rseed(ix,iy);

for (ir=0;ir ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download