rm(list=ls(all=TRUE))
library(LearnBayes)  ##multivariate M-H
library(SuppDists)##Sampling inverse Gaussian
library(bayesm)###截尾正态分布
#### real values
theta<-1.2 # 1.000125 10
sigma<-1.5# 0.019   5
mu<-1.81  #0.00395   5
T<-200 #外循环
###########################################
t<-c(0,100,400,1000,2000) #seq(0,1500,by=50)
lt<-length(t) ##
lamda<-rep(0,lt-1)
for(i in 2:length(t)){
  lamda[i-1]<-t-t[i-1]
}
S0<-50
S<-c(83,133,173) ##Stress
s<-numeric(3)
for(i in 1:3){   ##标准化stress
  s=(1/S-1/S0)/(1/S[3]-1/S0)
}
h<-s
h<-numeric(3)
for(i in 1:3){
  h=s/s[1]
}
sn<-c(50,50,50)#######每个压力水平下的样本个数
N<-sn*(lt-1)  ##每个压力水平下总共的样本点数
sN<-sum(N)

########## defined function #################
lamdai<-sum(lamda)*sn


#####定义后验所需函数######
fphi<-function(params){#(theta,sigma,mu)
  alpha<-sum(y1^2/(2*params[1]^h[1]*lamda))+sum(y2^2/(2*params[1]^h[2]*lamda))+sum(y3^2/(2*params[1]^h[3]*lamda))
  beta<-sum(y1/params[1]^h[1])+sum(y2/params[1]^h[2])+sum(y3/params[1]^h[3])
  gamma<-lamdai[1]/(2*params[1]^h[1])+lamdai[2]/(2*params[1]^h[2])+lamdai[3]/(2*params[1]^h[3])
  phi<-(alpha-beta*params[3]+params[3]^2*gamma)/(2*params[2]^2)                  
  return(phi)
}

######pij####
transpostpij<-function(theta0){
  theta<-exp(theta0[1])+1
  sigma<-exp(theta0[2])
  mu<-exp(theta0[3])
  logphi2<--sN*log(2*pi)/2-(sN+1)*log(sigma)-
    (h[1]*N[1]+h[2]*N[2]+h[3]*N[3]+2)*log(theta)/2-
    fphi(c(theta,sigma,mu))-lamdayijk+theta0[1]+theta0[2]+theta0[3]
  return(logphi2)
}


####################compute################
start=c(0,1,1)
iteration=2000
probs=seq(0,1,by=0.01)

paramspij<-matrix(0,nrow = T,ncol = 3)

CPsigma<-numeric(4)    ##用于累计sigma覆盖个数
CPtheta<-numeric(4)   ##用于累计theta覆盖个数
CPmu<-numeric(4)   ##用于累计mu覆盖个数

for (k in 1:T) {
  ##sampling delta y1
  y1<-matrix(0,ncol = sn[1],nrow = lt-1)
  for(j in 1:sn[1]){
    for (i in 1:(lt-1)) {
      y1[i,j]<-rnorm(1,mu*lamda, (sigma*theta^h[1]*lamda)^1/2)
    }
  }
  ##sampling delta y2
  y2<-matrix(0,ncol = sn[2],nrow = lt-1)
  for(j in 1:sn[2]){
    for (i in 1:(lt-1)) {
      y2[i,j]<-rnorm(1,mu*lamda, (sigma*theta^h[2]*lamda)^1/2)
    }
  }
  ##sampling delta y3
  y3<-matrix(0,ncol = sn[3],nrow = lt-1)
  for(j in 1:sn[3]){
    for (i in 1:(lt-1)) {
      y3[i,j]<-rnorm(1,mu*lamda, (sigma*theta^h[3]*lamda)^1/2)
    }
  }

  lamdayijk<-sum(sn)*(log(lamda[1]*lamda[2]*lamda[3]*lamda[4]))
  #############################################

  #######pij
  laplacefit=laplace(transpostpij,start)###这一步是套用别人的公式
  proposal=list(var=laplacefit$var,scale=1.5)###运行报错就是在两步
  sj=rwmetrop(transpostpij,proposal,laplacefit$mode,iteration)
  Btheta<-exp(sj$par[,1])+1
  Bsigma<-exp(sj$par[,2])
  Bmu<-exp(sj$par[,3])
  paramspij[k,]<-c(mean(Btheta),mean(Bsigma),mean(Bmu))
  #####
  cisigma<-c(quantile(Bsigma,0.025),quantile(Bsigma,0.975))
  citheta<-c(quantile(Btheta,0.025),quantile(Btheta,0.975))
  cimu<-c(quantile(Bmu,0.025),quantile(Bmu,0.975))
  if(sigma>cisigma[1]&&sigma<cisigma[2]){CPsigma[1]<-CPsigma[1]+1}
  if(theta>citheta[1]&&theta<citheta[2]){CPtheta[1]<-CPtheta[1]+1}
  if(mu>cimu[1]&&mu<cimu[2]){CPmu[1]<-CPmu[1]+1}
  plot(1,1,xlab=k)
}

######

这是报错
Error in chol.default(proposal$var) :
  the leading minor of order 3 is not positive definite
In addition: Warning message:
In log(det(h)) : NaNs produced