首页 > 解决方案 > 评估观察向量的 logpdf,其中每个观察具有不同的平均参数

问题描述

Julia 新手,只是尝试实现一个基本的贝叶斯模型。我想评估每个数据点的对数似然,其中每个数据点根据其对应的协变量具有不同的平均参数,而不必在所有数据点上实现 for 循环。

using Distributions 

y = -50:1:49
a = 1
b = 1
N = 100
x = rand(Normal(0, 1), N)
mu = a .+ b.*x
sigma = 5
# Can we evaluate the logpdf of every point in one call to logpdf without doing a for loop
loglikelihood = logpdf(Normal(mu, sigma), y)

MethodError: no method matching Normal(::Vector{Float64}, ::Int64)

编辑:我想澄清一下,mu上面指定的是一个与 相同维度的向量y,而不是在迭代过程中logpdf使用该函数评估每个观察值Normal(::Real, ::Real),我想处理一些影响 logpdf(Normal(::Array, ::Real), ::Array). 我在以下块中提供的代码通过对观察结果求对数似然之和来满足我的要求,但我宁愿不必转换为多元分布。

using LinearAlgebra

logpdf(MvNormal(mu, diagm(repeat([sigma], outer=N))), y)

谢谢你的帮助。

标签: statisticsjuliabayesian

解决方案


您的代码实际上并未运行,因为存在未定义的变量 ( a, b, y)。但总的来说,您所要求的开箱即用:

julia> using Distributions

julia> μ = 2.0; σ = 3.0;

julia> logpdf(Normal(μ, σ), 0:0.5:4)
9-element Vector{Float64}:
 -2.2397730440950046
 -2.1425508218727822
 -2.073106377428338
 -2.0314397107616715
 -2.0175508218727827
 -2.0314397107616715
 -2.073106377428338
 -2.1425508218727822
 -2.2397730440950046

在这里,我在值 0、0.5、1、...、3.5、4 处获取日志 pdf。这是有效的,因为有一种方法logpdf需要AbstractArray第二个参数:

julia> @which logpdf(Normal(μ, σ), 0:0.5:4)
logpdf(d::UnivariateDistribution{S} where S<:ValueSupport, X::AbstractArray) in Distributions at deprecated.jl:70

julia> @which logpdf(Normal(μ, σ), 0.5)
logpdf(d::Normal, x::Real) in Distributions at ...\Distributions\bawf4\src\univariate\continuous\normal.jl:105

正如您在那里看到的那样,该方法签名实际上已被弃用。让我们从 Julia 开始depwarn=yes,看看弃用通知:

$> julia --depwarn=yes

julia> using Distributions

julia> logpdf(Normal(), 1:10)
┌ Warning: `logpdf(d::UnivariateDistribution, X::AbstractArray)` is deprecated, use `logpdf.(d, X)` instead.
│   caller = top-level scope at REPL[4]:1
└ @ Core REPL[4]:1

这告诉您的是,实际上您不需要接受数组的方法签名,因为 Julia 的内置广播语法 - 在函数调用上附加一个点 - 免费为您提供。回到第一个例子:

julia> logpdf.(Normal(μ, σ), 0:0.5:4)
9-element Vector{Float64}:
 -2.2397730440950046
 -2.1425508218727822
 -2.073106377428338
 -2.0314397107616715
 -2.0175508218727827
 -2.0314397107616715
 -2.073106377428338
 -2.1425508218727822
 -2.2397730440950046

在这里,我实际上是在调用该logpdf(d::Normal, x::Real)方法,但.之后logpdf将函数元素应用于 range 0:0.5:4

广播语法还扩展到构造函数,因此您可以使用它来构造具有不同均值的多个正态分布:

julia> μ = rand(3)
3-element Vector{Float64}:
 0.5341692431981215
 0.5696647074299088
 0.3021675356902611

julia> Normal.(μ, 5)
3-element Vector{Normal{Float64}}:
 Normal{Float64}(μ=0.5341692431981215, σ=5.0)
 Normal{Float64}(μ=0.5696647074299088, σ=5.0)
 Normal{Float64}(μ=0.3021675356902611, σ=5.0)

这就是上面的错误告诉你的——Normal构造函数不接受向量作为第一个元素,而是一个值。如果要将其应用于多个值,只需广播!


推荐阅读