Максимизация ожиданий (EM) - это своего рода вероятностный метод классификации данных. Пожалуйста, поправьте меня, если я ошибаюсь, если это не классификатор.
Что такое интуитивное объяснение этой техники EM? Что здесь expectation и что maximized?
Примечание: код этого ответа можно найти здесь.
Предположим, у нас есть некоторые данные, взятые из двух разных групп, красного и синего:
Здесь мы можем видеть, какая точка данных принадлежит красной или синей группе. Это позволяет легко найти параметры, которые характеризуют каждую группу. Например, среднее значение красной группы составляет около 3, среднее значение синей группы - около 7 (и мы могли бы найти точное среднее, если бы захотели).
Вообще говоря, это называется оценкой максимального правдоподобия. Учитывая некоторые данные, мы вычисляем значение параметра (или параметров), который лучше всего объясняет эти данные.
Теперь представьте, что мы не можем видеть, какое значение было выбрано из какой группы. Все выглядит фиолетовым для нас:
Здесь мы знаем, что существует две группы значений, но мы не знаем, к какой группе принадлежит какое-либо конкретное значение.
Можем ли мы оценить средства для красной и синей групп, которые лучше всего соответствуют этим данным?
Да, часто мы можем! Максимизация ожиданий дает нам способ сделать это. Общая идея алгоритма такова:
Эти шаги нуждаются в дополнительном объяснении, поэтому я рассмотрю проблему, описанную выше.
Я буду использовать Python в этом примере, но код должен быть довольно простым для понимания, если вы не знакомы с этим языком.
Предположим, у нас есть две группы, красная и синяя, со значениями, распределенными как на рисунке выше. В частности, каждая группа содержит значение, полученное из нормального распределения со следующими параметрами:
import numpy as np
from scipy import stats
np.random.seed(110) # for reproducible results
# set parameters
red_mean = 3
red_std = 0.8
blue_mean = 7
blue_std = 2
# draw 20 samples from normal distributions with red/blue parameters
red = np.random.normal(red_mean, red_std, size=20)
blue = np.random.normal(blue_mean, blue_std, size=20)
both_colours = np.sort(np.concatenate((red, blue))) # for later use...
Вот изображение этих красных и синих групп снова (чтобы избавить вас от необходимости прокручивать вверх):
Когда мы видим цвет каждой точки (то есть, к какой группе она принадлежит), очень легко оценить среднее значение и стандартное отклонение для каждой группы. Мы просто передаем значения красного и синего во встроенные функции в NumPy. Например:
>>> np.mean(red)
2.802
>>> np.std(red)
0.871
>>> np.mean(blue)
6.932
>>> np.std(blue)
2.195
Но что, если мы не можем видеть цвета точек? То есть вместо красного или синего каждая точка была окрашена в фиолетовый цвет.
Чтобы попытаться восстановить параметры среднего и стандартного отклонения для красной и синей групп, мы можем использовать максимизацию ожидания.
Наш первый шаг (шаг 1 выше) состоит в том, чтобы угадать значения параметров для каждой средней группы и стандартное отклонение. Нам не нужно догадываться разумно; мы можем выбрать любые числа, которые нам нравятся:
# estimates for the mean
red_mean_guess = 1.1
blue_mean_guess = 9
# estimates for the standard deviation
red_std_guess = 2
blue_std_guess = 1.7
Эти оценки параметров дают кривые колокола, которые выглядят так:
Это плохие оценки. Оба средства (вертикальные пунктирные линии) выглядят, например, далеко от любого вида "середины" для разумных групп точек. Мы хотим улучшить эти оценки.
Следующий шаг (шаг 2) - это вычисление вероятности появления каждой точки данных в соответствии с предположениями текущего параметра:
likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)
likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)
Здесь мы просто поместили каждую точку данных в функцию плотности вероятности для нормального распределения, используя наши текущие предположения о среднем и стандартном отклонении для красного и синего. Это говорит нам, например, что с нашими текущими догадками точка данных в 1.761, скорее всего, будет красной (0.189), чем синей (0.00003).
Для каждой точки данных мы можем превратить эти два значения правдоподобия в веса (шаг 3), чтобы они суммировали с 1 следующим образом:
likelihood_total = likelihood_of_red + likelihood_of_blue
red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total
С нашими текущими оценками и нашими недавно вычисленными весами мы можем теперь вычислить новые оценки для среднего и стандартного отклонения красной и синей групп (шаг 4).
Мы дважды вычисляем среднее и стандартное отклонение, используя все точки данных, но с разными весами: один раз для красных весов и один раз для синих.
Ключевой момент интуиции заключается в том, что чем больше вес цвета в точке данных, тем больше точка данных влияет на следующие оценки параметров этого цвета. Это имеет эффект "вытягивания" параметров в правильном направлении.
def estimate_mean(data, weight):
"""
For each data point, multiply the point by the probability it
was drawn from the colour distribution (its "weight").
Divide by the total weight: essentially, we're finding where
the weight is centred among our data points.
"""
return np.sum(data * weight) / np.sum(weight)
def estimate_std(data, weight, mean):
"""
For each data point, multiply the point squared difference
from a mean value by the probability it was drawn from
that distribution (its "weight").
Divide by the total weight: essentially, we're finding where
the weight is centred among the values for the difference of
each data point from the mean.
This is the estimate of the variance, take the positive square
root to find the standard deviation.
"""
variance = np.sum(weight * (data - mean)**2) / np.sum(weight)
return np.sqrt(variance)
# new estimates for standard deviation
blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)
red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)
# new estimates for mean
red_mean_guess = estimate_mean(both_colours, red_weight)
blue_mean_guess = estimate_mean(both_colours, blue_weight)
У нас есть новые оценки для параметров. Чтобы улучшить их снова, мы можем вернуться к шагу 2 и повторить процесс. Мы делаем это до тех пор, пока оценки не сойдутся или не будет выполнено некоторое количество итераций (шаг 5).
Для наших данных первые пять итераций этого процесса выглядят так (последние итерации выглядят сильнее):
Мы видим, что средние значения уже сходятся по некоторым значениям, и формы кривых (определяемых стандартным отклонением) также становятся более стабильными.
Если мы продолжим 20 итераций, мы получим следующее:
Процесс EM сходится к следующим значениям, которые оказываются очень близкими к фактическим значениям (где мы можем видеть цвета - без скрытых переменных):
| EM guess | Actual | Delta
----------+----------+--------+-------
Red mean | 2.910 | 2.802 | 0.108
Red std | 0.854 | 0.871 | -0.017
Blue mean | 6.838 | 6.932 | -0.094
Blue std | 2.227 | 2.195 | 0.032
В приведенном выше коде вы, возможно, заметили, что новая оценка для стандартного отклонения была вычислена с использованием предыдущей итерационной оценки для среднего значения. В конечном счете, не имеет значения, вычисляем ли мы сначала новое значение для среднего значения, поскольку мы просто находим (взвешенную) дисперсию значений вокруг некоторой центральной точки. Мы еще увидим, что оценки параметров сходятся.
EM - это алгоритм максимизации функции правдоподобия, когда некоторые из переменных в вашей модели ненаблюдаются (т.е. когда у вас есть скрытые переменные).
Вы можете спросить, если мы просто пытаемся максимизировать функцию, почему бы нам просто не использовать существующие механизмы для максимизации функции. Ну, если вы попытаетесь максимизировать это, взяв производные и установив их на ноль, вы обнаружите, что во многих случаях условия первого порядка не имеют решения. Там проблема с курицей и яйцом, чтобы решить ваши параметры модели, вам нужно знать распределение ваших ненаблюдаемых данных; но распределение ваших ненаблюдаемых данных является функцией ваших параметров модели.
E-M пытается обойти это, итеративно угадывая распределение для ненаблюдаемых данных, а затем оценивая параметры модели, максимизируя то, что является нижней границей функции фактического правдоподобия, и повторяется до сближения:
Алгоритм ЭМ
Начните с предположения о значениях параметров вашей модели
E-step: для каждого дататота, у которого отсутствуют значения, используйте уравнение модели для решения для распределения недостающих данных, учитывая текущее предположение параметров модели и учитывая наблюдаемые данные (обратите внимание, что вы решаете для распределения для каждого отсутствующего значения, а не для ожидаемого значения). Теперь, когда у нас есть распределение для каждого отсутствующего значения, мы можем вычислить математическое ожидание функции правдоподобия относительно ненаблюдаемых переменных. Если наше предположение для параметра модели было правильным, это ожидаемое значение будет фактической вероятностью наших наблюдаемых данных; если параметры были неправильными, это будет только нижняя граница.
М-шаг: теперь, когда у нас есть ожидаемая функция правдоподобия без каких-либо незаметных переменных в ней, максимизируйте функцию, как и в полностью наблюдаемом случае, чтобы получить новую оценку параметров вашей модели.
Повторяйте до сходимости.
Вот простой рецепт, чтобы понять алгоритм максимизации ожиданий:
1- Прочтите этот учебный материал по EM от Do и Batzoglou.
2- У вас могут быть вопросительные знаки в вашей голове, посмотрите объяснения на этой странице обмена математическим стеком.
3- Посмотрите на этот код, который я написал на Python, который объясняет пример в учебном материале по EM из пункта 1:
Предупреждение: код может быть грязным/неоптимальным, так как я не являюсь разработчиком Python. Но это делает работу.
import numpy as np
import math
#### E-M Coin Toss Example as given in the EM tutorial paper by Do and Batzoglou* ####
def get_mn_log_likelihood(obs,probs):
""" Return the (log)likelihood of obs, given the probs"""
# Multinomial Distribution Log PMF
# ln (pdf) = multinomial coeff * product of probabilities
# ln[f(x|n, p)] = [ln(n!) - (ln(x1!)+ln(x2!)+...+ln(xk!))] + [x1*ln(p1)+x2*ln(p2)+...+xk*ln(pk)]
multinomial_coeff_denom= 0
prod_probs = 0
for x in range(0,len(obs)): # loop through state counts in each observation
multinomial_coeff_denom = multinomial_coeff_denom + math.log(math.factorial(obs[x]))
prod_probs = prod_probs + obs[x]*math.log(probs[x])
multinomial_coeff = math.log(math.factorial(sum(obs))) - multinomial_coeff_denom
likelihood = multinomial_coeff + prod_probs
return likelihood
# 1st: Coin B, {HTTTHHTHTH}, 5H,5T
# 2nd: Coin A, {HHHHTHHHHH}, 9H,1T
# 3rd: Coin A, {HTHHHHHTHH}, 8H,2T
# 4th: Coin B, {HTHTTTHHTT}, 4H,6T
# 5th: Coin A, {THHHTHHHTH}, 7H,3T
# so, from MLE: pA(heads) = 0.80 and pB(heads)=0.45
# represent the experiments
head_counts = np.array([5,9,8,4,7])
tail_counts = 10-head_counts
experiments = zip(head_counts,tail_counts)
# initialise the pA(heads) and pB(heads)
pA_heads = np.zeros(100); pA_heads[0] = 0.60
pB_heads = np.zeros(100); pB_heads[0] = 0.50
# E-M begins!
delta = 0.001
j = 0 # iteration counter
improvement = float('inf')
while (improvement>delta):
expectation_A = np.zeros((5,2), dtype=float)
expectation_B = np.zeros((5,2), dtype=float)
for i in range(0,len(experiments)):
e = experiments[i] # i'th experiment
ll_A = get_mn_log_likelihood(e,np.array([pA_heads[j],1-pA_heads[j]])) # loglikelihood of e given coin A
ll_B = get_mn_log_likelihood(e,np.array([pB_heads[j],1-pB_heads[j]])) # loglikelihood of e given coin B
weightA = math.exp(ll_A) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of A proportional to likelihood of A
weightB = math.exp(ll_B) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of B proportional to likelihood of B
expectation_A[i] = np.dot(weightA, e)
expectation_B[i] = np.dot(weightB, e)
pA_heads[j+1] = sum(expectation_A)[0] / sum(sum(expectation_A));
pB_heads[j+1] = sum(expectation_B)[0] / sum(sum(expectation_B));
improvement = max( abs(np.array([pA_heads[j+1],pB_heads[j+1]]) - np.array([pA_heads[j],pB_heads[j]]) ))
j = j+1
Технически термин "EM" немного недосказан, но я предполагаю, что вы ссылаетесь на метод анализа кластеров Gaussian Mixture Modeling, который является примером общего принципа EM.
Собственно, EM-кластерный анализ не является классификатором. Я знаю, что некоторые люди считают, что кластеризация является "неконтролируемой классификацией", но на самом деле кластерный анализ - это совсем другое.
Ключевое различие и большая классификация недоразумений, которые всегда имеют кластерный анализ, таковы: в кластерном анализе, нет "правильного решения" . Это метод обнаружения знаний, на самом деле это означает найти что-то новое! Это делает оценку очень сложной. Его часто оценивают с использованием известной классификации в качестве ссылки, но это не всегда уместно: классификация, которую вы имеете, может или не может отражать то, что находится в данных.
Позвольте мне привести пример: у вас большой набор данных клиентов, включая гендерные данные. Метод, который разбивает этот набор данных на "мужские" и "женские", оптимален, когда вы сравниваете его с существующими классами. В "предсказании" способ мышления это хорошо, как и для новых пользователей, теперь вы можете предсказать их пол. В способе "открытия знаний" это действительно плохо, потому что вы хотели обнаружить новую структуру данных. Способ, который, например, разделить данные на пожилых людей и детей, однако будет забивать так же хуже, как и в отношении класса мужчины/женщины. Однако это было бы отличным результатом кластеризации (если возраст не был указан).
Теперь вернемся к EM. По существу, предполагается, что ваши данные состоят из множества многомерных нормальных распределений (обратите внимание, что это очень сильное предположение, в частности, когда вы фиксируете количество кластеров!). Затем он пытается найти локальную оптимальную модель для этого с помощью , чередуя улучшающую модель и назначение объекта модели.
Для получения наилучших результатов в контексте классификации выберите количество кластеров, большее, чем число классов, или даже примените кластеризацию только к отдельным классам (чтобы выяснить, существует ли какая-либо структура внутри класса!).
Скажите, что вы хотите обучить классификатор, чтобы рассказать обо всех "машинах", "велосипедах" и "грузовиках". Существует мало пользы в предположении, что данные состоят из ровно 3 нормальных распределений. Однако вы можете предположить, что существует более одного типа автомобилей (и грузовиков и велосипедов). Поэтому вместо обучения классификатора для этих трех классов вы кластерируете автомобили, грузовики и велосипеды в 10 кластеров каждый (или, может быть, 10 автомобилей, 3 грузовика и 3 велосипеда, независимо от того), затем тренируйте классификатор, чтобы рассказать обо всех этих 30 классах, а затем объединить результат класса с исходными классами. Вы также можете обнаружить, что есть один кластер, который особенно сложно классифицировать, например Trikes. Это несколько автомобилей и несколько мотоциклов. Или грузовые автомобили, которые больше похожи на негабаритные автомобили, чем грузовики.
Другие ответы хороши, я постараюсь предоставить еще одну перспективу и решить интуитивную часть вопроса.
Алгоритм EM (Ожидание-Максимизация) является вариантом класса итеративных алгоритмов с использованием duality
Выдержка (акцент мой):
В математике двойственность, вообще говоря, переводит понятия, теорем или математических структур в другие понятия, теоремы или структуры, индивидуально, часто (но не всегда) с помощью средств операции инволюции: если двойственный от А есть В, то двойственный от В A. Такие инволюции иногда имеют фиксированные точки, так что двойственные из A есть сам A
Обычно двойной B объекта A связан с A каким-то образом, что сохраняет некоторую симметрию или совместимость. Например, AB = const
Примеры итеративных алгоритмов, использующих двойственность (в предыдущем смысле):
Аналогично, алгоритм ЭМ также можно рассматривать как два двойных шага максимизации:
.. [EM] рассматривается как максимизирующая совместную функцию параметров и распределение по ненаблюдаемым переменным. Е-шаг максимизирует эта функция относительно распределения по ненаблюдаемому переменные; М-шаг по параметрам.
В итеративном алгоритме, использующем двойственность, существует явное (или неявное) предположение о равновесной (или фиксированной) точке сходимости (для ЭМ это доказывается с использованием неравенства Йенсена)
Таким образом, схема таких алгоритмов:
Примечание, что, когда такой алгоритм сходится к (глобальному) оптимуму, он нашел конфигурацию, которая лучше в обоих смыслах (т.е. как в x домена и параметров и домена y). Однако алгоритм может просто найти локальный оптимум, а не глобальный оптимум.
я бы сказал, что это интуитивное описание схемы алгоритма
Для статистических аргументов и приложений другие ответы дали хорошие объяснения (проверьте также ссылки в этом ответе)
Принятый ответ ссылается на Chuong EM Paper, который делает достойную работу, объясняя EM. Существует также видео youtube, которое более подробно объясняет эту статью.
Напомним, вот сценарий:
1st: {H,T,T,T,H,H,T,H,T,H} 5 Heads, 5 Tails; Did coin A or B generate me?
2nd: {H,H,H,H,T,H,H,H,H,H} 9 Heads, 1 Tails
3rd: {H,T,H,H,H,H,H,T,H,H} 8 Heads, 2 Tails
4th: {H,T,H,T,T,T,H,H,T,T} 4 Heads, 6 Tails
5th: {T,H,H,H,T,H,H,H,T,H} 7 Heads, 3 Tails
Two possible coins, A & B are used to generate these distributions.
A & B have an unknown parameter: their bias towards heads.
We don't know the biases, but we can simply start with a guess: A=60% heads, B=50% heads.
В случае первого пробного вопроса, интуитивно мы думаем, что B сгенерировал его, так как пропорция совпадений с вершинами B очень хорошо... но это значение было всего лишь догадкой, поэтому мы не можем быть уверены.
Имея это в виду, мне нравится думать о решении EM следующим образом:
Это может быть упрощение (или даже принципиально неправильное на некоторых уровнях), но я надеюсь, что это поможет на интуитивном уровне!
EM используется для максимизации вероятности модели Q со скрытыми переменными Z.
Это итеративная оптимизация.
theta <- initial guess for hidden parameters
while not converged:
#e-step
Q(theta'|theta) = E[log L(theta|Z)]
#m-step
theta <- argmax_theta' Q(theta'|theta)
е-шаг: при текущей оценке Z рассчитывает ожидаемую функцию правдоподобия
т-шаг: найти theta, который максимизирует это Q
Пример GMM:
e-step: оценивать присвоения меток для каждого datapoint, учитывая текущую оценку параметра gmm
m-step: максимизируйте новую тету с учетом новых настроек метки
K-средство также является алгоритмом EM, и на K-средствах много объясняющих анимаций.
Используя ту же статью Do и Batzoglou, приведенную в ответе Zhubarb, я применил EM для этой проблемы в Java. Комментарии к его ответу показывают, что алгоритм застревает в локальном оптимуме, что также происходит с моей реализацией, если параметры thetaA и thetaB одинаковы.
Ниже приведен стандартный вывод моего кода, показывающий сходимость параметров.
thetaA = 0.71301, thetaB = 0.58134
thetaA = 0.74529, thetaB = 0.56926
thetaA = 0.76810, thetaB = 0.54954
thetaA = 0.78316, thetaB = 0.53462
thetaA = 0.79106, thetaB = 0.52628
thetaA = 0.79453, thetaB = 0.52239
thetaA = 0.79593, thetaB = 0.52073
thetaA = 0.79647, thetaB = 0.52005
thetaA = 0.79667, thetaB = 0.51977
thetaA = 0.79674, thetaB = 0.51966
thetaA = 0.79677, thetaB = 0.51961
thetaA = 0.79678, thetaB = 0.51960
thetaA = 0.79679, thetaB = 0.51959
Final result:
thetaA = 0.79678, thetaB = 0.51960
Ниже приведена моя Java-реализация EM для решения проблемы в (Do and Batzoglou, 2008). Основной частью реализации является цикл для запуска EM до тех пор, пока параметры не сближаются.
private Parameters _parameters;
public Parameters run()
{
while (true)
{
expectation();
Parameters estimatedParameters = maximization();
if (_parameters.converged(estimatedParameters)) {
break;
}
_parameters = estimatedParameters;
}
return _parameters;
}
Ниже приведен весь код.
import java.util.*;
/*****************************************************************************
This class encapsulates the parameters of the problem. For this problem posed
in the article by (Do and Batzoglou, 2008), the parameters are thetaA and
thetaB, the probability of a coin coming up heads for the two coins A and B,
respectively.
*****************************************************************************/
class Parameters
{
double _thetaA = 0.0; // Probability of heads for coin A.
double _thetaB = 0.0; // Probability of heads for coin B.
double _delta = 0.00001;
public Parameters(double thetaA, double thetaB)
{
_thetaA = thetaA;
_thetaB = thetaB;
}
/*************************************************************************
Returns true if this parameter is close enough to another parameter
(typically the estimated parameter coming from the maximization step).
*************************************************************************/
public boolean converged(Parameters other)
{
if (Math.abs(_thetaA - other._thetaA) < _delta &&
Math.abs(_thetaB - other._thetaB) < _delta)
{
return true;
}
return false;
}
public double getThetaA()
{
return _thetaA;
}
public double getThetaB()
{
return _thetaB;
}
public String toString()
{
return String.format("thetaA = %.5f, thetaB = %.5f", _thetaA, _thetaB);
}
}
/*****************************************************************************
This class encapsulates an observation, that is the number of heads
and tails in a trial. The observation can be either (1) one of the
experimental observations, or (2) an estimated observation resulting from
the expectation step.
*****************************************************************************/
class Observation
{
double _numHeads = 0;
double _numTails = 0;
public Observation(String s)
{
for (int i = 0; i < s.length(); i++)
{
char c = s.charAt(i);
if (c == 'H')
{
_numHeads++;
}
else if (c == 'T')
{
_numTails++;
}
else
{
throw new RuntimeException("Unknown character: " + c);
}
}
}
public Observation(double numHeads, double numTails)
{
_numHeads = numHeads;
_numTails = numTails;
}
public double getNumHeads()
{
return _numHeads;
}
public double getNumTails()
{
return _numTails;
}
public String toString()
{
return String.format("heads: %.1f, tails: %.1f", _numHeads, _numTails);
}
}
/*****************************************************************************
This class runs expectation-maximization for the problem posed by the article
from (Do and Batzoglou, 2008).
*****************************************************************************/
public class EM
{
// Current estimated parameters.
private Parameters _parameters;
// Observations from the trials. These observations are set once.
private final List<Observation> _observations;
// Estimated observations per coin. These observations are the output
// of the expectation step.
private List<Observation> _expectedObservationsForCoinA;
private List<Observation> _expectedObservationsForCoinB;
private static java.io.PrintStream o = System.out;
/*************************************************************************
Principal constructor.
@param observations The observations from the trial.
@param parameters The initial guessed parameters.
*************************************************************************/
public EM(List<Observation> observations, Parameters parameters)
{
_observations = observations;
_parameters = parameters;
}
/*************************************************************************
Run EM until parameters converge.
*************************************************************************/
public Parameters run()
{
while (true)
{
expectation();
Parameters estimatedParameters = maximization();
o.printf("%s\n", estimatedParameters);
if (_parameters.converged(estimatedParameters)) {
break;
}
_parameters = estimatedParameters;
}
return _parameters;
}
/*************************************************************************
Given the observations and current estimated parameters, compute new
estimated completions (distribution over the classes) and observations.
*************************************************************************/
private void expectation()
{
_expectedObservationsForCoinA = new ArrayList<Observation>();
_expectedObservationsForCoinB = new ArrayList<Observation>();
for (Observation observation : _observations)
{
int numHeads = (int)observation.getNumHeads();
int numTails = (int)observation.getNumTails();
double probabilityOfObservationForCoinA=
binomialProbability(10, numHeads, _parameters.getThetaA());
double probabilityOfObservationForCoinB=
binomialProbability(10, numHeads, _parameters.getThetaB());
double normalizer = probabilityOfObservationForCoinA +
probabilityOfObservationForCoinB;
// Compute the completions for coin A and B (i.e. the probability
// distribution of the two classes, summed to 1.0).
double completionCoinA = probabilityOfObservationForCoinA /
normalizer;
double completionCoinB = probabilityOfObservationForCoinB /
normalizer;
// Compute new expected observations for the two coins.
Observation expectedObservationForCoinA =
new Observation(numHeads * completionCoinA,
numTails * completionCoinA);
Observation expectedObservationForCoinB =
new Observation(numHeads * completionCoinB,
numTails * completionCoinB);
_expectedObservationsForCoinA.add(expectedObservationForCoinA);
_expectedObservationsForCoinB.add(expectedObservationForCoinB);
}
}
/*************************************************************************
Given new estimated observations, compute new estimated parameters.
*************************************************************************/
private Parameters maximization()
{
double sumCoinAHeads = 0.0;
double sumCoinATails = 0.0;
double sumCoinBHeads = 0.0;
double sumCoinBTails = 0.0;
for (Observation observation : _expectedObservationsForCoinA)
{
sumCoinAHeads += observation.getNumHeads();
sumCoinATails += observation.getNumTails();
}
for (Observation observation : _expectedObservationsForCoinB)
{
sumCoinBHeads += observation.getNumHeads();
sumCoinBTails += observation.getNumTails();
}
return new Parameters(sumCoinAHeads / (sumCoinAHeads + sumCoinATails),
sumCoinBHeads / (sumCoinBHeads + sumCoinBTails));
//o.printf("parameters: %s\n", _parameters);
}
/*************************************************************************
Since the coin-toss experiment posed in this article is a Bernoulli trial,
use a binomial probability Pr(X=k; n,p) = (n choose k) * p^k * (1-p)^(n-k).
*************************************************************************/
private static double binomialProbability(int n, int k, double p)
{
double q = 1.0 - p;
return nChooseK(n, k) * Math.pow(p, k) * Math.pow(q, n-k);
}
private static long nChooseK(int n, int k)
{
long numerator = 1;
for (int i = 0; i < k; i++)
{
numerator = numerator * n;
n--;
}
long denominator = factorial(k);
return (long)(numerator / denominator);
}
private static long factorial(int n)
{
long result = 1;
for (; n >0; n--)
{
result = result * n;
}
return result;
}
/*************************************************************************
Entry point into the program.
*************************************************************************/
public static void main(String argv[])
{
// Create the observations and initial parameter guess
// from the (Do and Batzoglou, 2008) article.
List<Observation> observations = new ArrayList<Observation>();
observations.add(new Observation("HTTTHHTHTH"));
observations.add(new Observation("HHHHTHHHHH"));
observations.add(new Observation("HTHHHHHTHH"));
observations.add(new Observation("HTHTTTHHTT"));
observations.add(new Observation("THHHTHHHTH"));
Parameters initialParameters = new Parameters(0.6, 0.5);
EM em = new EM(observations, initialParameters);
Parameters finalParameters = em.run();
o.printf("Final result:\n%s\n", finalParameters);
}
}