Я использую разложение Холецкого для моделирования коррелированных случайных величин с учетом матрицы корреляции. Дело в том, что результат никогда не воспроизводит структуру корреляции так, как он задан. Вот небольшой пример на Python, чтобы проиллюстрировать ситуацию.
import numpy as np
n_obs = 10000
means = [1, 2, 3]
sds = [1, 2, 3] # standard deviations
# generating random independent variables
observations = np.vstack([np.random.normal(loc=mean, scale=sd, size=n_obs)
for mean, sd in zip(means, sds)]) # observations, a row per variable
cor_matrix = np.array([[1.0, 0.6, 0.9],
[0.6, 1.0, 0.5],
[0.9, 0.5, 1.0]])
L = np.linalg.cholesky(cor_matrix)
print(np.corrcoef(L.dot(observations)))
Это печатает:
[[ 1. 0.34450587 0.57515737]
[ 0.34450587 1. 0.1488504 ]
[ 0.57515737 0.1488504 1. ]]
Как вы можете видеть, аппроксимируемая оценочная матрица корреляции резко отличается от предыдущей. Есть ли ошибка в моем коде, или есть какая-то альтернатива использованию декомпозиции Холецкого?
редактировать
Прошу прощения за этот беспорядок. Я не думал, что в коде и / или в том, как применялась декомпозиция Холецкого, произошла ошибка из-за неправильного понимания материала, который я изучал ранее. На самом деле я был уверен, что сам метод не должен быть точным, и я был в порядке с этим, пока ситуация не заставила меня задать этот вопрос. Спасибо, что указали на заблуждение, которое у меня было. Я отредактировал название, чтобы лучше отразить реальную ситуацию, предложенную @Silverfish.
источник
Ответы:
Подход, основанный на разложении Холецкого, должен работать, он описан здесь и показан в ответе Марка Л. Стоуна, опубликованном почти одновременно с этим ответом.
Пример в
R
(извините, я не использую то же программное обеспечение, которое вы использовали в вопросе):Вас также может заинтересовать этот пост и этот пост .
источник
Люди, вероятно, найдут вашу ошибку намного быстрее, если вы объясните, что вы сделали со словами и алгеброй, а не с кодом (или, по крайней мере, написали ее с использованием псевдокода).
Вы, кажется, делаете эквивалент этого (хотя, возможно, транспонирован):
высчитываетY= L X
Что вы должны сделать, это:
Есть много объяснений этого алгоритма на сайте. например
Как генерировать коррелированные случайные числа (с учетом средних, дисперсий и степени корреляции)?
Могу ли я использовать метод Холецкого для генерации коррелированных случайных величин с заданным средним значением?
Этот обсуждает это непосредственно в терминах желаемой ковариационной матрицы, а также дает алгоритм для получения желаемой выборки ковариации:
Генерация данных с заданной выборочной ковариационной матрицей
источник
Нет ничего плохого в факторизации Холецкого. В вашем коде есть ошибка. Смотрите редактирование ниже.
Вот код и результаты MATLAB, сначала для n_obs = 10000, как у вас, а затем для n_obs = 1e8. Для простоты, поскольку это не влияет на результаты, я не беспокоюсь о средствах, т. Е. Делаю их нулями. Обратите внимание, что chol MATLAB производит верхний треугольный фактор Холецки R матрицы M, так что R '* R = M. numpy.linalg.cholesky производит нижний треугольный фактор Холецки, поэтому необходима корректировка по сравнению с моим кодом; но я считаю, что ваш код в этом отношении хорош.
Изменить: я нашел вашу ошибку. Вы неправильно применили стандартное отклонение. Это эквивалент того, что вы сделали, что неправильно.
источник
Резюме не о коде, но я был заинтригован, увидев, как это выглядело бы после всех хороших ответов и, в частности, @Mark L. Stone. Фактический ответ на вопрос предоставляется на его посту (пожалуйста, укажите его в случае сомнений). Я перемещаю эту дополнительную информацию сюда, чтобы облегчить поиск этого поста в будущем. Не преуменьшая других превосходных ответов, после ответа Марка это завершает проблему, исправляя сообщение в ОП.
Источник
В питоне:
IN [R]:
источник
Как уже показали другие: холеский работает. Вот фрагмент кода, который очень короткий и очень близок к псевдокоду: кодовый фрагмент в MatMate:
источник