Monte Carlo Simulation Summary and Explanation

Conceptual: The Monte Carlo simulation I performed works on a format called simulated annealing. The process mirrors the metallurgy technique to improve the crystalline nature of a metal.

At high temperatures, metal atoms possess too much kinetic energy to settle down into any lattice framework. If they reach a potential energy minimum, they jump right back out again. The structure is completely random. Quickly cooling a metal forces atoms to stop where ever they are, leaving the mostly random state intact. But cooling the metal slowly enough allows atoms to move around to find the lowest potential minimum before they settle for good. After some cooling, when they find a potential minimum they could stay there, but temperature could allow them increase in potential energy. This increase could allow them to travel to a global minimum instead of a mere local minimum.

So our test point randomly explores the survey quality function (I'm not sure if what I used could be considered a likelihood function) acting like a atom exploring a potential energy function. Our test point starts hot, with a high probability of accepting even an unfavorable new point during its random walk. During the cooling timescale, it steadily calms its jumpy nature and hones in on a minimum.

Practical: To best understand the actual implementation, I will walk through the computing steps.

After receiving a starting point, it randomly permutes both the number of observations and the number of stars which would characterize a survey. Safety measures must be implemented to keep those parameters from exceeding possible bounds and ending the simulation.
In the beginning I imported an array of the magnitudes. Given the magnitudes of the brightest required number of stars, 60 hours of observing time, and each star needing to get the specific number of observations, the lowest instrument noise is calculated. Stellar jitter for each specific star is randomly pulled from a Gaussian centered at 15m/s with a standard deviation of 3 m/s. Stellar jitter and instrument error are added in quadrature.

Next, the parameters are sent into a survey function. 100 simulated surveys are run and the average number of planets discovered is returned. Within a survey, a star is randomly given a planet with a probability equal to .07 times its mass (the masses were also imported earlier on). The planet's properties are randomly selected from a distribution modeled on available data for intermediate mass stars. Each star receives n observations at times randomly pulled from the possible time vector. Three possible days exist per month over a 36 month period. "Observations" consist of calculating the radial velocity due to the planet and adding Gaussian jitter according to the noise. A planet is considered detected if the RMS of the observations is three times greater than the noise level.

Once back in the main simulation, the new set of parameters is kept if either

A. the number of planets produced is greater than that before
B. A uniformly randomly generated number is less than exp( - diff in planets/ Energy)

Hence at high energies, there exists a high chance of excepting the number, and at low energies the probability shrinks. I set the energy timescale to a decaying exponential according to advice on the internet. The energy starts at 20, and each step the energy decreases by 0.2%. I selected the numbers to start very jumpy and end very cold.

If the point is accepted, it is returned to a file.
This process is repeated 5000 times.

(And I just learned a good deal about making aesthetic graphs in MATLAB)