评估观测向量的 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)

感谢您的帮助。

解决方法

您的代码实际上并未运行，因为存在未定义的变量（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

这里我得到的 log pdf 值为 0,0.5,1,...,3.5,4。这是有效的，因为 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

不过，正如您在那里看到的，该方法签名实际上已被弃用。让我们以 depwarn=yes 开始 Julia 以查看弃用通知：

$> 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 将函数元素应用于范围 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 构造函数不接受向量作为第一个元素，而是接受单个值。如果要将其应用于多个值，只需广播即可！

bayesian julia statistics