# Generate N parameter vectors from the prior.

set.seed(1)           # Random number seed
N <- 1000000          # Size of sample from prior

p  <- runif(N)        # Mixing proportion 
v  <- rgamma(N,3)/20  # Reciprocal of variance 
u1 <- rnorm(N,20,10)  # Mean for component no. 1
u2 <- rnorm(N,20,10)  # Mean for component no. 2

# Find the likelihoods for these parameter vectors.

galaxies <- scan("galaxy.data")

l <- rep(1,N)         # Vector of likelihoods

for (x in galaxies)
{ l <- l * ( p * sqrt(v/(2*pi)) * exp(-v*(x-u1)^2/2)
       + (1-p) * sqrt(v/(2*pi)) * exp(-v*(x-u2)^2/2) );
}

# Find mean likelihood, and associated std. error.

m <- mean(l)          # Estimate for marg. likelihood
e <- sqrt(var(l)/N)   # Standard error for estimate

# Print estimate & std. error, in direct and log form.

cat (m, e, log(m), e/m, "\n")