<jupyter_start><jupyter_text>Create a synthetic 2D stellar population with some noise<jupyter_code>samples = np.random.multivariate_normal([-0, -0], [[1, 0],[0, 1]], 600)
x = samples[:,0]
y = samples[:,1]
r = np.sqrt(x**2 + y **2)
plt.plot(x,y,'.',ms=2)
plt.xlim(-4,4)
plt.ylim(-4,4)
plt.xlabel('X')
plt.ylabel('Y')
# bin into radial bins of given step size
step = 0.3
bins = np.arange(0, np.max(r), step)
annuli_density = []

for ii in bins:
    cts = len(r[r < (ii+step)]) - len(r[r < ii])
    area = np.pi*(((ii+step)**2) - (ii**2))
    surf_dens = cts/area # counts divided by area
    annuli_density.append(surf_dens)
plt.plot(bins, annuli_density, 'o', label='Subcomponent #1 Only')
plt.xlabel('Radius (px)')
plt.ylabel('Surface Density (stars/px^2)')
plt.legend(loc=1, numpoints=1)<jupyter_output><empty_output><jupyter_text>Perform kernel density estimation on the single-component toy model:<jupyter_code>xmin = x.min()
xmax = x.max()
ymin = y.min()
ymax = y.max()
X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([x, y])
kernel = stats.gaussian_kde(values)
Z = np.reshape(kernel(positions).T, X.shape)
fig, ax = plt.subplots()
ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
          extent=[xmin, xmax, ymin, ymax])
ax.plot(x, y, 'k.', markersize=2)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
plt.show()<jupyter_output><empty_output><jupyter_text>Test of fit: Anderson-Darling for a normal distribution<jupyter_code>stats.anderson(annuli_density, dist='norm')<jupyter_output><empty_output><jupyter_text>Statistic is > 0.6, which is encouraging. This should ideally be larger than that for a synthetic population with a little subpopulation. That is, a larger value for a single population suggests that we are more likely to reject the hypothesis that there are two populations. Let's try injecting a subpopulation that has 40/70 = 0.57 the surface density of the core at roughly the same spatial scale of 0.25 of the cluster radial extent. (3 pc/12 pc, with the extent of the subpopulation about 3 pc in size too, so 0.25 of the full cluster extent.)<jupyter_code>sm_samples = np.random.multivariate_normal([1, -1], [[0.1, 0],[0, 0.1]], 100)
sm_x = sm_samples[:,0]
sm_y = sm_samples[:,1]
sm_r = np.sqrt(sm_x**2 + sm_y**2)
plt.plot(x, y, '.',ms=2)
plt.plot(sm_x, sm_y, 'r.', ms=2)
plt.xlim(-4,4)
plt.ylim(-4,4)
plt.xlabel('X')
plt.ylabel('Y')
# bin into radial bins of given a given step size, in this case scaled to match that of Holland+2000

step = 0.3
sm_bins = np.arange(np.sqrt(2.0), np.max(sm_r), step)
sm_annuli_density = []

for ii in sm_bins:
    cts = len(sm_r[sm_r < (ii+step)]) - len(sm_r[sm_r < ii])
    area = np.pi*(((ii+step)**2) - (ii**2))
    surf_dens = cts/area 
    sm_annuli_density.append(surf_dens)

plt.plot(sm_bins, sm_annuli_density, 'ro', label='Subcomponent #2 Only')
plt.xlabel('Radius (px)')
plt.ylabel('Surface Density (stars/px^2)')
plt.legend(loc=1, numpoints=1)
sm_bins<jupyter_output><empty_output><jupyter_text>Combine the two populations (in a really simple way) and re-attempt whether or not this is drawn from a Gaussian distribution:<jupyter_code>x_new = np.concatenate([x, sm_x])
y_new = np.concatenate([y, sm_y])
r_new = np.concatenate([r, sm_r])
# bin into radial bins of given step size
step = 0.3
bins_new = np.arange(0, np.max(r_new), step)
annuli_density_new = []

for ii in bins_new:
    cts = len(r_new[r_new < (ii+step)]) - len(r_new[r_new < ii])
    area = np.pi*(((ii+step)**2) - (ii**2))
    surf_dens = cts/area 
    annuli_density_new.append(surf_dens)

plt.plot(bins_new, annuli_density_new, 'ko', label='Combined Populations')
plt.xlabel('Radius (px)')
plt.ylabel('Surface Density (stars/px^2)')
plt.legend(loc=1, numpoints=1)<jupyter_output><empty_output><jupyter_text>Perform KDE on the two-component toy model:<jupyter_code>xmin_new = x_new.min()
xmax_new = x_new.max()
ymin_new = y_new.min()
ymax_new = y_new.max()
X_new, Y_new = np.mgrid[xmin_new:xmax_new:100j, ymin_new:ymax_new:100j]
positions_new = np.vstack([X_new.ravel(), Y_new.ravel()])
values_new = np.vstack([x_new, y_new])
kernel_new = stats.gaussian_kde(values_new)
Z_new = np.reshape(kernel_new(positions_new).T, X_new.shape)

fig, ax = plt.subplots(1,2, figsize=(10,10))
ax[0].imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
          extent=[xmin, xmax, ymin, ymax])
# ax.plot(x_new, y_new, 'k.', markersize=2)
ax[0].plot(x_new, y_new, 'k.', markersize=2)
# ax[0].plot(sm_x, sm_y, 'r.', markersize=2)
ax[0].set_xlim([xmin, xmax])
ax[0].set_ylim([ymin, ymax])

ax[1].imshow(np.rot90(Z_new), cmap=plt.cm.gist_earth_r,
          extent=[xmin_new, xmax_new, ymin_new, ymax_new])
# ax.plot(x_new, y_new, 'k.', markersize=2)
ax[1].plot(x, y, 'k.', markersize=2)
ax[1].plot(sm_x, sm_y, 'r.', markersize=2)
ax[1].set_xlim([xmin, xmax])
ax[1].set_ylim([ymin, ymax])

plt.show()<jupyter_output><empty_output><jupyter_text>Test Anderson-Darling again<jupyter_code>stats.anderson(annuli_density_new, dist='norm')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17096390321301072, 'fraction_numerical': 0.02598175327251091, 'mean_word_length': 5.5408560311284045, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '654890', 'license': 'MIT', 'n_tokens_mistral': 2221, 'n_tokens_neox': 2090, 'n_words': 454, 'path': 'toymodels_for_hypothesis_testing.ipynb', 'repo_name': 'stephtdouglas/M44_subcluster', 'size': '5042'}