Работа с нулями при преобразовании журнала (с учетом других ограничений) в R - PullRequest
Новая
       19

Работа с нулями при преобразовании журнала (с учетом других ограничений) в R

0 голосов
/ 12 апреля 2020

У меня есть ситуация, когда мне нужно преобразовать свои данные для работы с ним, но в моей матрице есть нули. В дополнение к нулям, матрица, которую я имею, также извлечена из распределения Дирихле, что означает, что матрица имеет ограничение, что все суммы по столбцам должны складываться до 1. Вот данные: 

> q[1:10, 1:5]
            V1          V2          V3           V4           V5
1  0.534410243 0.009358740 0.011295181 0.2141751740 0.0030129254
2  0.026653603 0.372426720 0.447847534 0.0179177507 0.4072904477
3  0.193317915 0.003605024 0.003186611 0.4832114736 0.0007095471
4  0.111881585 0.000000000 0.000000000 0.2296213741 0.0119233461
5  0.089696570 0.591163629 0.509774416 0.0032542030 0.5535847030
6  0.007543558 0.000000000 0.000000000 0.0364907757 0.0013148362
7  0.004862942 0.000000000 0.002123909 0.0146682272 0.0004053690
8  0.009276195 0.011710457 0.014367894 0.0000000000 0.0000000000
9  0.006903171 0.004314528 0.011404455 0.0000000000 0.0126889937
10 0.015454219 0.007420903 0.000000000 0.0006610215 0.0090698319

Обратите внимание, что все столбцы q складываются в один 

> colSums(q)[1:5]
V1 V2 V3 V4 V5 
 1  1  1  1  1

Мне нужно взять log (q) примерно так: 

> log(q)[1:10, 1:5]
           V1         V2         V3         V4         V5
1  -0.6265915 -4.6714446 -4.4833791 -1.5409610 -5.8048438
2  -3.6248309 -0.9877150 -0.8033024 -4.0219634 -0.8982287
3  -1.6434192 -5.6254270 -5.7487974 -0.7273009 -7.2508837
4  -2.1903142       -Inf       -Inf -1.4713235 -4.4292569
5  -2.4113228 -0.5256624 -0.6737870 -5.7278079 -0.5913405
6  -4.8870614       -Inf       -Inf -3.3106958 -6.6340431
7  -5.3261117       -Inf -6.1544972 -4.2220715 -7.8107129
8  -4.6803038 -4.4472730 -4.2427592       -Inf       -Inf
9  -4.9757744 -5.4457675 -4.4737512       -Inf -4.3670203
10 -4.1698733 -4.9034546       -Inf -7.3217241 -4.7028016

Как видите, есть тонна - Значения Inf, которые портят мои расчеты. Я думал о замене нуля очень маленькими числами, но тогда сумма больше не равна 1 по столбцам. Как мне написать код для построения альтернативной матрицы q, что 1) не имеет нулевых значений и поэтому обходит проблему log (0), а 2) все еще имеет столбцы, которые складываются в единицу, и не меняет базовое распределение данные по строкам?

Большое спасибо!

Редактировать: Чтобы обеспечить немного более широкий контекст: мне нужно взять преобразование журнала, так как я передаю вывод в расчет журнала функция правдоподобия. В моем приложении я перенастраиваем параметры логарифмической вероятности распределения Дирихле, поэтому я не вызываю обобщенную c функцию логарифмической вероятности из пакета.

Вот как выглядит моя общая функция: 

llikelihood = function(alpha0, beta, q, d, n) {
  llike = n*(lgamma(alpha0) - sum_a(alpha0, beta, d) + sum_b (alpha0, beta, q, d, n))
  return(llike)
}

sum_a = function(alpha0, beta, d) {
  sum_a = 0
  for (i in 1:d) {
    sum_a = sum_a + lgamma(alpha0*beta[i]) 
  }
  return(sum_a)
}

# returns the output to summation from 1 to k of (alpha0*beta[i] - 1)*log(x_i)
sum_b = function(alpha0, beta, q, d, n) { 
  # replace zero values
  sum_b = 0
  # find the log q
  logq = log(q)
  qlog = apply(logq, 1, sum)
  #  for each column, sum up the draws
  for (i in 1:d) {
    sum_b = sum_b + (alpha0*beta[i] - 1)*1/n*qlog[i]
  }

  # apply(log(q), 2, sum)
  return(sum_b)
}

Здесь sum_b - это место, где я вычисляю логарифм (q), как указано выше. Как видите, моя проблема в том, что мне нужно избавиться от нулей, нормализовать данные до одного, а затем вести журнал этого. Как я могу написать код, который делает это эффективно? Я предполагаю, что это будет похоже на Laplace Smoothing, но я мало что знаю об этом, и я новичок в программировании на R. Большое спасибо за комментарии!

1 Ответ

1 голос
/ 12 апреля 2020

1) Вы можете попробовать другие преобразования, которые не возвращают -Inf с нулями, например квадрат root или куб root.

2) Нормализовать результат из 1) путем деления всех элементов на суммы их столбцов. 

set.seed(123)
X <- t(rdirichlet(4, alpha=c(1,0,2,1)))
X
           [,1]      [,2]      [,3]       [,4]
[1,] 0.03562445 0.3384606 0.5700819 0.01357789
[2,] 0.00000000 0.0000000 0.0000000 0.00000000
[3,] 0.64748450 0.2927702 0.3297736 0.88378152
[4,] 0.31689105 0.3687692 0.1001445 0.10264059

colSums(X)
# [1] 1 1 1 1

Шаг 1) Квадрат root. 

X2 <- sqrt(X); X2
          [,1]      [,2]      [,3]      [,4]
[1,] 0.1887444 0.5817737 0.7550377 0.1165242
[2,] 0.0000000 0.0000000 0.0000000 0.0000000
[3,] 0.8046642 0.5410824 0.5742592 0.9400966
[4,] 0.5629308 0.6072637 0.3164561 0.3203757

Шаг 2) Нормализовать 

X3 <- sweep(X2, 2, colSums(X2), FUN="/"); X3

          [,1]      [,2]      [,3]       [,4]
[1,] 0.1212746 0.3362621 0.4587794 0.08462201
[2,] 0.0000000 0.0000000 0.0000000 0.00000000
[3,] 0.5170236 0.3127428 0.3489340 0.68271531
[4,] 0.3617018 0.3509952 0.1922865 0.23266269

> colSums(X3)
[1] 1 1 1 1
...