Домашнее задание № 3

Восстановление параметров статистического распределения

В этом задании вы будете восстанавливать параметры смешанных статистических распределений, используя два метода: метод максимального правдоподобия и EM-метод. Научным заданием будет выделение двух рассеянных скоплений h и χ Персея в звёздном поле.

Дедлайн 22 апреля в 23:55

Вы должны реализовать следующие алгоритмы в файле mixfit.py:

1. Метод максимального правдоподобия для смеси двух одномерных нормальных распределений. Напишите функцию max_likelihood(x, tau, mu1, sigma1, mu2, sigma2, rtol=1e-3), где x — массив данных, остальные позиционные аргументы — начальные приближения для искомых параметров распределения, rtol — относительная точность поиска параметров, функция должна возвращать кортеж из параметров распределения в том же порядке, что и в сигнатуре функции. Для оптимизации разрешается использовать scipy.optimize.

2. Expectation-maximization метод для той же задачи: em_double_gauss(x, tau, mu1, sigma1, mu2, sigma2, rtol=1e-3).

3. EM-метод для смеси трех нормальных распределений — τ1 N(µ1, σ) + τ2 N(µ2, σ) + (1-τ1-τ2) N(0, σ0). Напишите функцию em_double_cluster(x, tau1, tau2, muv, mu1, mu2, sigma02, sigmax2, sigmav2, rtol=1e-5), где

   • x — массив (N, 4) состоящий из двух столбцов координат звезд, и двух столбцов собственых движений звезд. Собственными движениями называется скорость движения звезды в картинной плоскости в угловых единицах;
   • τ1 (tau1), τ2 (tau2) — относительные количества звезд в первом и втором скоплениях. Обратите внимание, что относительное количество звезд поля может быть найдена из условия нормировки (1-τ1-τ2);
   • µ1 — четырехмерный параметр распределения звезд в первом скоплении. Состоит из среднего положения звезд скопления mu1 и среднего собственного движения скопления muv;
   • µ2 — четырехмерный параметр распределения звезд во втором скоплении. Состоит из среднего положения звезд скопления mu2 и среднего собственного движения скопления muv. Обратите внимание, что предполагается, что скорость движения скоплений — одинаковая;
   • σ — диагональная матрица 4x4. Первые два компоненты которой равны sigmax2 и отвечают за разброс положения звезд скоплений вокруг среднего, а последние два компонента sigmav2 отвечают за разброс собственных движений звезд скоплений относительно среднего. Иначе говоря, предполагается, что разброс по координатам имеет симметричный характер и одинаковую дисперсию по обоим направлениям. Аналогичное предположение действует для разброса проекций скоростей.
   • σ0 — диагональная матрица 2x2, элементы которой равны sigma02. Разброс собственных движений звезд поля. Обратите внимание, что предполагается, что среднее собственное движение звезд поля — нулевое.

Для отладки ваших алгоритмов и самопроверки, вы можете написать несколько модульных тестов в файле test.py. Для генерации случайных данных подойдут, например, методы scipy.stats.multivariate_normal.rvs и scipy.stats.uniform.rvs.

4. Примените последний EM-метод в файле per.py для решения задачи нахождения центров и относительного числа звёзд в скоплениях h и χ Персея. Вам понадобиться модуль astroquery.vizier для того, чтобы загрузить координаты звёзд поля. Координаты центра двойного скопления 02h21m00s +57d07m42s, используйте поле размером 1.0 x 1.0 градус. Пример запроса для получения данных:

   import numpy as np
   import astropy.units as u
   from astropy.coordinates import SkyCoord
   from astroquery.vizier import Vizier
   
   center_coord = SkyCoord('02h21m00s +57d07m42s')
   vizier = Vizier(
       columns=['RAJ2000', 'DEJ2000', 'pmRA', 'pmDE'],
       column_filters={'BPmag': '<16', 'pmRA': '!=', 'pmDE': '!='}, # число больше — звёзд больше
       row_limit=10000
   )
   stars = vizier.query_region(
       center_coord,
       width=1.0 * u.deg,
       height=1.0 * u.deg,
       catalog=['I/350'], # Gaia EDR3
   )[0]
   
   ra = stars['RAJ2000']._data   # прямое восхождение, аналог долготы
   dec = stars['DEJ2000']._data  # склонение, аналог широты
   x1 = (ra - ra.mean()) * np.cos(dec / 180 * np.pi) + ra.mean()
   x2 = dec
   v1 = stars['pmRA']._data
   v2 = stars['pmDE']._data

5. В файл per.json сохрнаите результаты:

   {
     "size_ratio": 1.2,
     "motion": {"ra": -0.63, "dec": -1.15},
     "clusters": [
       {
         "center": {"ra": 35.35, "dec": 57.07},
       },
       {
         "center": {"ra": 36.36, "dec": 57.57},
       }
     ]
   }

6. В файле per.png изобразите график рассеяния точек звёздного поля, и график рассеяния собственных движений. На обоих графиках для каждой точки отметьте цветом условную вероятность принадлежности к одному из скоплений (используйте вычисленные параметры T1 и T2). Обозначьте скопления окружностями с центром в центре распределения и радиусом равным стандартном отклонению (корень из дисперсии).