初学R，想做单元素Gibbs，用的是股票收益率

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

有人有编程吗？望高手指导

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Gibbs 股票收益率 股票收益 bbs 收益率 收益率 元素

看看r的相关help吧

是Gibbs抽样吧

这里有一个例子，你看看有没有用

一个核电厂有十个水泵，已知数据各水泵出故障次数，以及观察到出故障的时间
      y <- c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)t <- c(94, 16, 63, 126, 5, 31, 1, 1, 2, 10)
建立关于水泵故障次数的泊松分布，λi 为单位时间每个水泵的故障次数。似然函数有：
                   ∏i=1~10 Poisson(λi *ti)
假设Gamma(α, β) 是 λ 的先验分布。α=1.8. 所有的 λi都来自于这个分布。
假设Gamma(γ , δ) 是 β 的先验分布。γ=0.01 , δ=1
所以有：
p(λ, β|y, t) ∝  (∏i=1~10 Poisson(λi *ti)*Gamma(α, β))*Gamma(γ , δ)
经过整理有：
p(λ, β|y, t) ∝ ( ∏i=1~10 λi^(yi +α−1)*e^(−(ti +β)λi))*(β^(10α+γ−1)e^(−δβ))
通过上述联合分布，可以求得满条件分布:
p(λi |λ−i , β, y, t) ∝ λi^(yi +α−1)*e^(−(ti +β)λi)
p(β|λ, y, t)∝ β^(10α+γ−1)e^−β(δ+Sumi=1~10 λi)

gibbs<-function(n.sims,beta.start,alpha,gamma,delta,y,t,burnin=0,thin=1){beta.draws<-c()lambda.draws<-matrix(NA,nrow=n.sims,ncol=length(y))beta.cur<-beta.startlambda.update<-function(alpha,beta,y,t){rgamma(length(y),y+alpha,t+beta)}beta.update<-function(alpha,gamma,delta,lambda,y){rgamma(1,length(y)*alpha+gamma,delta+sum(lambda))}for(i in 1:n.sims){lambda.cur<-lambda.update(alpha=alpha,beta=beta.cur,y=y,t=t)beta.cur<- beta.update(alpha=alpha,gamma=gamma,delta=delta,lambda=lambda.cur,y=y)if(i>burnin&&(i-burnin)%%thin==0){lambda.draws[(i-burnin)%/%thin,]<-lambda.curbeta.draws[(i-burnin)/thin]<-beta.cur}}return(list(lambda.draws=lambda.draws,beta.draws=beta.draws))}
其中lambda.update 方法和beta.update 方法是Gibbs Sampler过程。循环语句是对样本进行burn-in和thinning操作。

最后求取平均值，即所求参数λ, β的期望值。
posterior <- gibbs(n.sims = 10000, beta.start = 1, alpha = 1.8, gamma = 0.01, delta = 1, y = y, t = t)colMeans(posterior$lambda.draws)mean(posterior$beta.draws)