从概率后缀树获得对数似然

Question

这是我的代码：

library(RCurl)
library(TraMineR)
library(PST)

x <- getURL("https://gist.githubusercontent.com/aronlindberg/08228977353bf6dc2edb3ec121f54a29/raw/c2539d06771317c5f4c8d3a2052a73fc485a09c6/challenge_level.csv")
data <- read.csv(text = x)

# Load and transform data
data <- read.table("thread_level.csv", sep = ",", header = F, stringsAsFactors = F)

data.seq <- seqdef(data[2:nrow(data),2:ncol(data)], missing = "NA", right = "*")

# Make a tree
S1 <- pstree(data.seq, ymin = 0.05, L = 6, lik = TRUE, with.missing = F)
logLik(S1)

由于某种原因，它拒绝返回对数似然值？为什么会这样？如何获得对数似然值？

Answer 1

seqdef 命令中的 missing 和 right 参数值错误，这会导致 pstree 出现错误。

有

data.seq <- seqdef(data[2:nrow(data),2:ncol(data)], missing = NA, right= NA, nr = "*")
# Make a tree
S1 <- pstree(data.seq, ymin = 0.05, L = 6, lik = TRUE, with.missing = TRUE)
logLik(S1)

我们得到

'log Lik.' -31011.32 (df=47179)

请注意，由于您缺少值，因此我在 pstree 命令中设置了 with.missing = TRUE。

================

要忽略正确的缺失，请将right='DEL' 设置为seqdef。

seq <- seqdef(data[2:nrow(data),2:ncol(data)], missing = NA, right= "DEL")
S2 <- pstree(seq, ymin = 0.05, L = 6, lik = TRUE, with.missing = F)
logLik(S2)

我不知道 PST 计算为 logLik(S2) 以及为什么我们会在这里得到 NA。可以通过predict 函数获得使用树S2 生成数据的可能性，该函数返回数据中每个序列的可能性。那么数据的对数似然应该是

sum(log(predict(S2, seq)))

给了

 [>] 984 sequence(s) - min/max length: 1/32
 [!] sequences have unequal lengths
 [>] max. context length: L=6
 [>] found 1020 distinct context(s)
 [>] total time: 0.588 secs
[1] -4925.79

Answer 2

确实，在计算模型拟合不等长度序列的可能性时存在问题。这是固定的。新版 PST 包 (0.94) 将在几个小时内在 R-Forge 上提供，安装：

install.packages("PST", repos="http://R-Forge.R-project.org") 

后来在 CRAN 上。

请注意，由于您的序列不包含任何缺失值但长度不等，因此在使用pstree 函数时您不必设置with.missing=TRUE，在使用seqdef 时也不必设置任何选项。

现在运行以下代码时：

library(RCurl)
library(TraMineR)
library(PST)

x <- getURL("https://gist.githubusercontent.com/aronlindberg/08228977353bf6dc2edb3ec121f54a29/raw/c2539d06771317c5f4c8d3a2052a73fc485a09c6/challenge_level.csv")
data <- read.csv(text = x)

data.seq <- seqdef(data[2:nrow(data),2:ncol(data)])

# Make a tree
S1 <- pstree(data.seq, ymin = 0.05, L = 6) 

我明白了：

> S1@logLik
[1] -4925.79