在提问之前请确定你已经努力阅读了文档,并且尝试自己在互联网上搜索。
请尽可能提供你的demo代码或者GitHub的gist地址。
function standarddeviation(p)
# 读取数据
data = load("data($(p))OLSoutput.jld") # 需要用到JLD包
# 读取数据中的数字和矩阵
T = size(data["startingvalues(p = $(p))"]["Mydata"], 1)
K = data["startingvalues(p = $(p))"]["Total Variables"]
TT = T - p # 仅出现在S的方差表达式中
G = data["startingvalues(p = $(p))"]["G"]
S = data["startingvalues(p = $(p))"]["Std-Dev vector"]
R = data["startingvalues(p = $(p))"]["Correlation Matrix"]
OLS_residuals = data["startingvalues(p = $(p))"]["OLS_residuals"]
Smean = data["prior(p = $(p))"]["Smean"] # 为0
Ssd = data["prior(p = $(p))"]["Ssd"] # 为2
restrictionsSD = diag(data["restrictions(p = $(p))"]["Σ"])
# restrictionsSD[1] = 0
# S的方差参考 (http://mathworld.wolfram.com/StandardDeviationDistribution.html) 的第10式
V = 1/TT * ( TT - 1 - 2 * exp(2*(loggamma(TT/2)-loggamma((TT-1)/2)))) .* S
# 迭代开始
for i in 1:length(restrictionsSD)
#######################
# 只有i = 1是正确的结果
#######################
if restrictionsSD[i] != 0
# ⛔默认为全部通过
lowerbound = S[i] - 2 * sqrt(V[i])
upperbound = S[i] + 2 * sqrt(V[i])
# 👆在5%的显著性水平下,每一个标准差S[i]均相应地产生G-1个置信区间分段。Sgrid是分段的端点集(G个)
# 计算S.grid
Sgrid = LinRange(lowerbound, upperbound, G)
# 👆因为标准差S的个数有K个,所以总共有K组Sgrid,可以在最后进行保存
GSi = vec(zeros(G, 1)) # 👈设置初始值,以便后续进行迭代
#👆一共有K组GSi,可以每次保存迭代后的GSi
#👇迭代开始
for j in 1:G
# 首先构造S的均值向量
# μ = repeat([0], 2, 1) #这是❌的
μ = vec(repeat([0.0], K, 1)) # μ必须是向量类型!且Float64!
# 再构造关于S的协方差矩阵Sigmagrid
Sigrid = S # 防止原标准差向量S被替换后不可恢复,故意多出Sigrid这个中间变量
Sigrid[i] = Sgrid[j] # 随着迭代过程的更新,标准差向量S的第i个数被Sgrid重新赋值!
##############################################################################
# 👆问题出自这里。当i = 2时,迭代对象是S的第二个元素。
# 但S的第一个元素并不是data(6)OLSoutput.jld内标准差Std-Dev vector的第一个元素,
# 而是为i = 1, G = 100时的元素。
# 我的目的是令i = 2时,S的第一个元素自动迭代成Std-Dev vector的第一个元素。
##############################################################################
# ⛔此时Sigrid便是重新赋值后的结果,它不是原标准差向量S!
# 👇协方差矩阵Sigmagrid是随着迭代过程的更新而变化的
Sigmagrid = diagm(Sigrid) * R * diagm(Sigrid)
# ⛔为了避免问题:PosDefException: matrix is not Hermitian; Cholesky factorization failed.
# 特将Sigmagrid进行如下两行代码的处理:
Ω = cholesky(Hermitian(Sigmagrid))
Sigmagrid = Ω.L * Ω.U
# 通过上述变换,MvNormal(μ, Sigmagrid)得以接受
# 于是残差序列OLS_residuals的似然函数值为:
pdfresiduals = pdf(MvNormal(μ, Sigmagrid), OLS_residuals)
# 对残差序列的似然值逐个取对数:
η = log.(pdfresiduals)
# 将对数化后的似然值进行加总,可得似然概率值
loglikeli = sum(η)
# 👆这仅是一个数!
logSiprior = logpdf(LogNormal(Smean, Ssd), Sgrid[j])
# 👆用到了先验的超参数Smean和Ssd
GSi[j] = loglikeli + logSiprior
end