Planning A New Cluster Fitting Protocol

By studying a wider array of literature on the subject, I have become aware of multiple systemic errors in my code that need modification and improvements that I could enact. What I had been doing, though not wrong in the simplest sense, was an oversimplification and not representing the full knowledge that we have on the subject.

Therefore I attempt to make these changes:
1. Change to writing in IDL

It can read in the bulk isochrone files so i don't have to input each age individually
2. Use the Cross Entropy and bootstrapping methods
This will allow me to perform 3
3. Take larger and finer time and metallicity grids

4. Use U-B information as well
5. Factor in the probability that each star is a binary
Accomplished by changing the star's distribution function to a smeared Gaussian
6. Weight each star based similar to Mendieros
7. Use one reference source for the data
8. Factor in the Initial Mass Function

Unsure of how to accomplish this... maybe weighting the isochrone points
9. Perform for multiple theoretical isochrone databases to analyze systematic error arising from theoretical underpinnings

After all these, I will test my program using simulated clusters before applying them to the real thing.

Fitting isochrones to open cluster photometric data

Monteiro et al
A&A 2010

Monteiro solved for the parameters of an open cluster using the Cross Entropy method and analyzed the errors using the bootstrap method, instead of direct computation of the direct tau-squared space. This allows for a much larger grid. Right now, I limit my space to those near the published values.

The cross entropy method first randomly draws a certain number of parameters, and evaluates some function S(x) to be minimized at each point The points are ordered and the best N are selected out. These N define a mean and standard deviation around which to normally draw the next set of points. Weighting the mean and standard deviation with those of the previous round allows the process to converge less quickly and avoid local minima.

The function S(x) used to evaluate the quality of a point was the -log likelihood. Like Flannery & Johnson 1982, they use a nearest point estimator, only paying attention to the point which is closest. I believe that the nearby points should be considered as well and the probabilities added. This makes it more favorable for a star to exist in a area with a larger possible mass range, an important factor.

When the compute the probability for one point, they use the same Gaussian distribution I have been using, but also factor in fitting U-B as well, using all available data.

One of the most important topics in the paper dealt with determining which stars to exclude and giving stars weighting factors. The triangle of cluttered background stars near the bottom was completely cut out. The exact point of cutting did not effect the end result significantly. A cluster radius was computed using density given by positions in the image. Background stars outside this radius were excluded.

To compute the weighting factor, each gained a box with 3 sigma on each side. A Gaussian was computed using the mean and standard deviation for B, B-V, and U-B of stars in that box. It also included a term for its radius from the cluster center. This put importance on stars in the densest of the diagram and in the cluster center. The weight of each star was a factor in its probability.

After setting up the process, they performed multiple checks to assess the quality of the program. They created multiple synthetic clusters and checked if they could regain the same values. They took well-studied clusters and cross-referenced their values with those published. They studied the effect that different IMF's and binary % had on the results. IMF did not have an effect, but binary ratio could effect the distance up to 10%, though age remained the same. The technique was largely successful.

They used standard photometric errors (not those intrinsic to the data) as published by Mointinho 2001. They worked of Padova isochrones like the ones I'm using. They also fitted for E(B-V) in the beginning using a color-color diagram and only allowed it to vary 10% from that value. They also isolated data from only one source, and chose those with U-B data as well. They took 100% binary ratio with companions pulled from the same IMF.

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)

Why still a delta function PDF?

So...
Since the mass is still way too precise, I decided to move backward and look at the cluster fitting parameters. They turn out to be extremely precise too. Way too precise.

So, I went back and looked at the fitting paper that I read and I section that didn't make since before now makes since and I may be greatly wrong as to my process. They take the original isochrone, populate the HR diagram with many,many stars from this isochrone according to the initial mass function, and make a certain percentage of them unseen binaries according to some other distribution function. This creates a density on the HR diagram that is used to compute probabilities of a data point from that isochrone. Though the code would not be incredibly hard to write, it would be computationally difficult and greatly increase the time necessary for one cluster. Is that necessary and where I am going wrong?

Before I get too overwhelmed and confused, I think I'll try this on another cluster to compare results.

Better than a delta function PDF?

So I found an error in my code (I was dividing by the standard deviation in two separate places). The peak now does not move with a changing RMS.

And the error bounds are slightly better being twice as large, .002 solar masses instead of .001 ...

Survey Update
Stars:70
Clusters:24

All clusters with NGC# < 2287 have been completely rechecked for any hint of viability

Delta Function PDF?

So I am improving my cluster parameters fitting function, and creating my mass fittingfunction, and have a feeling something is wrong. All the code seems right. The cluster parameters returned visually fit the data well. The mass returned is within .1 of the rough distance only mass approximation. But the PDF for the mass is a delta function. I'm pretty sure we can't know the mass to a greater precision than .001 solar masses.

Later on, I have discovered that the PDF and even the maximum is highly dependent on the error. The RMS quoted for the star is .003, but MATLAB won't handle the extremely small probabilities that result from that. And obviously 67% of the bright stars are not .003 within a predicted isochrone. So I used .03, which gives the delta at 2.227, but I move it to .3, and it changes and moves the maximum to 2.5 (part of the loop that it obviously doesn't belong to).

I believe astrophysical properties, and not instrumentation, determine the majority of the error. Thus, I will fit all the clusters, and determine the error from tabulating the distribution of stars around their isochrones.

Note: I just eliminated fainter data since it included many background stars not even close to being on the main sequence, and I used CCD data.