MCMChybrid GP_5.4 R语言

已于 2022-07-09 15:40:48 修改
阅读量243 收藏
点赞数
分类专栏： MCMC HMC
于 2022-06-26 17:54:19 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/ygf666/article/details/125467646
版权

高斯过程 马尔科夫链 蒙特卡洛 探索阶段 采样阶段
关键词由CSDN通过智能技术生成
MCMC 同时被 2 个专栏收录
16 篇文章 2 订阅
订阅专栏
HMC
4 篇文章 0 订阅
订阅专栏

示例代码：

文章目录

library('MCMChybridGP')
mu1 <- c(-1, -1)
mu2 <- c(+1, +1)
sigma.sq <- 0.1225
ub <- c(1.5, 3)
X0 <- generateX0(lb=c(-2,-2), ub=ub)
f <- function(x) {
        px <- 1/4/pi/sqrt(sigma.sq) * exp(-1/2/sigma.sq *
            sum((x - mu1)^2)) + 1/4/pi/sqrt(sigma.sq) *
            exp(-1/2/sigma.sq * sum((x - mu2)^2))
        return(log(px))
    }
explore.out <- hybrid.explore(f, X0, ub=ub, n=150, graph=TRUE)
sample.out <- hybrid.sample(explore.out, n=500, graph=TRUE)

opar <- par(mfrow=c(2,1))
plot(density(sample.out$SAMP[,1]), xlab="x1", ylab="f(x)")
plot(density(sample.out$SAMP[,2]), xlab="x2", ylab="f(x)")
par(opar)

注：在运行之前，需要安装MCMChybridGP库,但是我在windows操作系统上装不了，在ubuntu上能装。

链接：https://pan.baidu.com/s/1wDBphwj_IOOb8GwWLHY-8A
提取码：00ri

接下来，对相关代码进行解释

1. 各种参数初始化

  • 导入库MCMChybridGP
  • 均值点（峰值点）1：（-1,-1） u 11 = − 1 , u 12 = − 1 u_{11}=-1,u_{12}=-1 u11​=−1,u12​=−1
  • 均值点（峰值点）2：（+1,+1） u 21 = + 1 , u 22 = + 1 u_{21}=+1,u_{22}=+1 u21​=+1,u22​=+1
  • 标准差：0.1225
  • 点的x坐标上限：1.5，y坐标上限：3
library('MCMChybridGP')
mu1 <- c(-1, -1)
mu2 <- c(+1, +1)
sigma.sq <- 0.1225
ub <- c(1.5, 3)

2.产生多个样本x

在上下限ub，lb的范围内产生若干样本

X0 <- generateX0(lb=c(-2,-2), ub=ub)

在这里插入图片描述
其中，generateX0函数如下所示

generateX0 <- function (lb, ub, n0 = 10, npool = 100)
{
  if(!is.null(dim(lb))) stop("lb in error")
  if(!is.null(dim(ub))) stop("ub in error")
  d <- length(lb)
  if(length(ub) != d) stop("Dimension mismatch for lb, ub")
  if(sum(abs(lb)+abs(ub)) == Inf) stop("Finite lb, ub required")
  ranX <- function() {
    X <- matrix(rep(NA, n0 * d), ncol = d)
    for (i in 1:d) {
      X[, i] <- sample(seq(lb[i], ub[i], length.out = n0))
    }
    return(X)
  }
  X <- Xbest <- as.matrix(ranX())
  if (n0 < 2) {
    return(Xbest)
  }
  dmax <- -Inf
  DIST <- function(x, y) sum(((x - y)/(ub - lb))^2)
  for (pool in 1:npool) {
    dmin <- Inf
    for (i in 1:(n0 - 1))
      for (j in (i + 1):n0) if (DIST(X[i, ], X[j, ]) < dmin) 
      dmin <- DIST(X[i, ], X[j, ])
    if (dmin > dmax) {
      dmax <- dmin
      Xbest <- X
    }
    if (pool < npool) 
      X <- as.matrix(ranX())
  }
  return(Xbest)
}

3. 定义密度函数（单项似然函数）

P ( x ) = P ( x 1 , x 2 ) = 1 4 π σ exp ⁡ ( − 1 2 σ ( x 1 − μ 11 ) 2 + ( x 1 − μ 12 ) 2 ) + 1 4 π σ exp ⁡ ( − 1 2 σ ( x 1 − μ 21 ) 2 + ( x 1 − μ 22 ) 2 ) P(x)=P(x_1,x_2)=\frac{1}{4\pi\sigma}\exp(-\frac{1}{2\sigma}(x_1-\mu_{11})^2+(x_1-\mu_{12})^2)+\frac{1}{4\pi\sigma}\exp(-\frac{1}{2\sigma}(x_1-\mu_{21})^2+(x_1-\mu_{22})^2) P(x)=P(x1​,x2​)=4πσ1​exp(−2σ1​(x1​−μ11​)2+(x1​−μ12​)2)+4πσ1​exp(−2σ1​(x1​−μ21​)2+(x1​−μ22​)2)
f ( x ) = log ⁡ P ( x ) f(x)=\log P(x) f(x)=logP(x)

f <- function(x) {
        px <- 1/4/pi/sqrt(sigma.sq) * exp(-1/2/sigma.sq *
            sum((x - mu1)^2)) + 1/4/pi/sqrt(sigma.sq) *
            exp(-1/2/sigma.sq * sum((x - mu2)^2))
        return(log(px))
    }

4. explore 探索阶段

explore.out <- hybrid.explore(f, X0, ub=ub, n=150, graph=TRUE)

接下来分析hybrid.explore 函数

4.1 函数参数

  • f ：密度函数 f ( x ) = log ⁡ P ( x ) f(x)=\log P(x) f(x)=logP(x)
  • X0: X 0 = x 1 , . . . , x 10 X_0={x_1,...,x_{10}} X0​=x1​,...,x10​ 10个观测输入
  • y0
  • n
  • L
  • delta
  • nchains
  • T.mult
  • lb
  • ub
  • maxleap
  • finetune
  • Tinclude
  • nswaps
  • graph
hybrid.explore <-
function (f, X0, ..., y0=NULL, n = 200, L = 1000, delta = 0.003, nchains = 5, 
    T.mult = 1.5, lb = NULL, ub = NULL, maxleap=0.95, finetune = FALSE,
    Tinclude = nchains, nswaps = choose(nchains, 2), graph = FALSE)

4.2 参数接收

  • x=X0
  • f n ( x ) = f ( x ) fn(x)=f(x) fn(x)=f(x)
# 随机矩阵传递
X <- as.matrix(X0)
# 函数传递
fn <- function(par) f(par, ...)

4.3 初始化

# 声明yfinal 和 Xfinal
yfinal <- Xfinal <- NULL
# X 的行数（观测点数）
n0 <- dim(X)[1]
# x的列数（观测点维度）
d <- dim(X)[2]
# Taccept：6 6 6 6 6 6 6 6 6 6 由n0个Tinclude+1组成
Taccept <- rep(Tinclude+1, n0)
# 声明接受
acceptance <- NULL

# 检查维度
if (d < 2) 
	stop("Points with dimension of at least 2 assumed for X")
if (d < 2) 
    graph <- FALSE
# 检查样本错误
if (n0 > 1) 
	for (i in 2:n0) 
		for (j in 1:(i - 1)) 
			if (sum(abs(X[i,] - X[j, ])) == 0) 
            stop(paste(sep = "", "Duplicate rows (", i, " & ", 
                j, ") in X not allowed"))
# 是否有界
BOUNDS <- FALSE
if ((!is.null(lb))  || (!is.null(ub))) BOUNDS <- TRUE
if(is.null(lb)) lb <- rep(-Inf, d)
# 下界传入
low <- lb
negInf <- function(X)
      as.double(apply(X, 2, function(x) 101*min(x)-100*max(x)))
low[!(low>-Inf)] <- negInf(X)[!(low>-Inf)]
if(is.null(ub)) ub <- rep(Inf, d)
# 上界引入
upp <- ub
posInf <- function(X)
      as.double(apply(X, 2, function(x) 101*max(x)-100*min(x)))
upp[!(upp<Inf)] <- posInf(X)[!(upp<Inf)]
# 界调整
if (BOUNDS) {
	if (!is.null(lb)) 
    	if (length(lb) == 1) 
        	lb <- rep(lb, d)
   	if (!is.null(ub)) 
       	if (length(ub) == 1) 
              ub <- rep(ub, d)
   	if ((length(lb) != d) || (length(lb) != d)) 
    	stop("lb, ub incorrect dimension")
   	if (is.na(sum(lb + ub))) 
      	stop("NaN lb, ub")
    }
  1. 没太看懂的

KEY.VARS <- NULL
# key.vars=1,2
key.vars <- 1:dim(X)[2]
for (i in 1:(length(key.vars) - 1)) 
	for (j in (i + 1):length(key.vars)) {
        pair <- c(key.vars[i], key.vars[j])
        KEY.VARS <- rbind(KEY.VARS, pair)
    }

在这里插入图片描述

  1. keyVars=1,2
if (graph) {
	if(dim(X)[2] <= 2) 
		nrefresh.graph <- Inf 
	else
		nrefresh.graph <- ceiling(10/nchains)#向上取整
	graph.template <- .hybrid.template# 实例化函数
    key.indx <- 1
    keyVars <- KEY.VARS[key.indx, ]
}

其中：.hybrid.template为画图模板函数如下

.hybrid.template <-
function (X, key.vars, lb, ub) 
{
    if(is.null(lb)) lb <- rep(-Inf, length(x1))
    if(is.null(ub)) ub <- rep(Inf, length(x1))
    # 按列分离
    # 第一列x1变量
    x1 <- X[, key.vars[1]]
    # 第二列x2变量
    x2 <- X[, key.vars[2]]
    #  获取边界
    low1 <- min(x1)
    upp1 <- max(x1)
    low2 <- min(x2)
    upp2 <- max(x2)
    extra1 <- extra2 <- 0.25
    if(sum(abs(lb)) < Inf) {
      low1 <- lb[key.vars[1]]
      low2 <- lb[key.vars[2]]
      extra1 <- 0.05
    }
    if(sum(abs(ub)) < Inf) {
      upp1 <- ub[key.vars[1]]
      upp2 <- ub[key.vars[2]]
      extra2 <- 0.05
    }
    l1 <- low1 - extra1 * (upp1 - low1)
    l2 <- low2 - extra1 * (upp2 - low2)
    u1 <- upp1 + extra2 * (upp1 - low1)
    u2 <- upp2 + extra2 * (upp2 - low2)
    # 获取合适的坐标范围
    xlim <- c(l1, u1)
    ylim <- c(l2, u2)
    # 设定标签x1,x2,和范围，以及步长显示
    plot(x1, x2, xlim=xlim, type="n", ylim=ylim,
         xlab=paste(sep="", "x", key.vars[1]),
         ylab=paste(sep="", "x", key.vars[2]))
    # 标出边界
    abline(v=lb[key.vars[1]], lty=2)
    abline(v=ub[key.vars[1]], lty=2)
    abline(h=lb[key.vars[2]], lty=2)
    abline(h=ub[key.vars[2]], lty=2)
    #设置图例
    legend(xlim[2], ylim[1], c("Accepted points", "Rejected points"), 
        pch=c(20, 4), pt.cex=c(1.2, 1), 
        cex=0.7, col=c("black", "red"), xjust=1, 
        yjust=0)
}
  1. 回火
    temp : 1.0000 1.5000 2.2500 3.3750 5.0625
# 温度取值T.mult=1.5
# 1.0000 1.5000 2.2500 3.3750 5.0625
Temp <- T.mult^(1:nchains - 1)
# 初始化交换记录
swaps.record <- NULL
# f的样本数
fcalls <- 0
if(!is.null(y0)) {
      y <- y0 
} 
# 获取样本输出 y
# 以及统计样本输出数量 fcalls
else {
	if (T.mult < 1) 
    	stop("T.mult < 1 not allowed")
    # n0为样本数
    y <- rep(NA, n0)
    date1 <- date2 <- date()
    for (j in 1:n0) {
    	functionOK <- FALSE
        try({
        	Xj <- X[j,]
            OUTSIDE_BOUNDS <- FALSE
            if (BOUNDS) {
            	for (i in 1:d) {
                	if (Xj[i] < low[i]) {
                  		OUTSIDE_BOUNDS <- TRUE
                  		Xj[i] <- low[i]
   					}
                	if (Xj[i] > upp[i]) {
                  		OUTSIDE_BOUNDS <- TRUE
                  		Xj[i] <- upp[i]
               		}
              	}
            }
           if(OUTSIDE_BOUNDS) stop("X0 violates bounds supplied")
           date1 <- date()
           y[j] <- fn(X[j,])
           date2 <- date()
           functionOK <- TRUE
        })
        if(!functionOK) stop("f(x) failed")
        if (is.na(y[j])) stop("f(x) NaN")
        if (y[j] == Inf) stop("Infinite f(x) not allowed")
    }
    fcalls <- n0
}
 # 给每条链选择初始x值
select <- sample(1:n0, nchains)
x <- X[select, ]
if (is.null(dim(x))) 
	x <- t(x)
# 给每条链初始x值对应的y值
fx <- y[select]
if (graph) {
	#初始化图
	graph.template(X, key.vars = keyVars, lb, ub)
    #画初始点
    points(X[, keyVars[1]], X[, keyVars[2]], pch = 20, col = "grey")
    # Xfinal=NULL
    points(Xfinal[, keyVars[1]], Xfinal[, keyVars[2]], pch = 20, col = "black")
    #设置标题
    title(c("", "", "", "Exploratory phase"), adj = 0, cex.main = 0.8)
    # nchains条链，每条连按照"green"->"turquoise"->"orange", 
                      #   ->"yellow"->"red"-> "grey"->"pink"
    clrs <- c("blue", rep(c("green", "turquoise", "orange", 
            "yellow", "red", "grey", "pink"), nchains))
    reds <- c("blue", rep(c("green", "turquoise", "orange", 
            "yellow", "red", "grey", "pink"), nchains))
    }
# 初始化计数器
count10 <- 0
nacc <- 0
# 初始化恢复计数器
refresh_count <- 1
# 参数数量为x维度+1
params <- rep(1, d+1)
#初始化平均接受率
mean.acceptance <- 0

4.4 重复采样

repeat { 内容如下  }

4.4.1 gp过程拟合

# 高斯过程拟合
GPfit <- GProcess(X, y, params = params, request.functions = FALSE,
                          finetune = finetune)

具体内容如下：

  • X ：样本输入
  • y ：样本输出
  • param：params <- rep(1, d+1)（表示3个待优化的参数）
  • request.functions：例子中是False
  • finetune ：例子中是False
GProcess <-
function (X, y, params = NULL, request.functions = TRUE, finetune=FALSE)
{   
    # X转为样本输入矩阵
    X <- as.matrix(X)
    # X行数即数据量为n，维度为d
    n <- dim(X)[1]
    d <- dim(X)[2]
    if(d < 2) stop("Points with dimension of at least 2 assumed for X")
    # y中心化（分布在0的2侧）
    # 求出y的均值mY
    mY <- mean(y)
    # y×=y-mY
    ystar <- y - mY
    if (is.null(params)) {
        params <- rep(1, 1 + d)
        finetune <- TRUE
    } 
    #例子中走else
    else factr <- 1e10
    if(finetune) 
    	factr <- 1e15

    Lreject <- -1e13
    Lok <- TRUE
    #  初始化协方差矩阵Sigma
    Sigma <- matrix(rep(0, n*n), ncol=n)
	
	# 待优化函数L=likehood(params)
    L <- function(params) {
        # 第一个参数
        lsq <- params[1]
        # 除第一个外，所有其他参数
        eta <- params[-1]
        # 求方差元素函数
        c.scalar <- function(xi, xj) {
            h <- abs(xi - xj)
            return(lsq * prod((1 + eta * h) * exp(-eta * h)))
        }
        # 求 x与X 协方差函数
        c.vec <- function(x) {
            vec <- rep(NA, n)
            for (j in 1:n) vec[j] <- c.scalar(x, X[j, ])
            return(vec)
        }
        # 协方差矩阵？？
        Sigma <- .C("GPcovar", as.double(X), as.integer(n), as.integer(d),
                    as.double(params), Sigma = as.double(Sigma))$Sigma
        # 修正协方差矩阵的格式
        dim(Sigma) <- c(n, n)
        # 判断 Sigma 是否能够chol分解为 R
        cholOK <- FALSE
        try({
            # R表示分解后的Sigma
            R <- chol(Sigma)
            cholOK <- TRUE
        }, silent=TRUE)
        if (!cholOK) {
            val <- Lreject * runif(1, 1, 10)
            Lok <- FALSE
            return(val)
        }
		# log|Sigma|见公式6右边的第一项https://blog.csdn.net/ygf666/article/details/125393906
        ldet.Sigma <- 2 * sum(log(diag(R)))
        # 求逆
        Sigma.inv <- chol2inv(R)
        # 求似然函数
        val <- -1/2 * ldet.Sigma - 1/2 * t(ystar) %*% Sigma.inv %*% 
            ystar
        if (is.na(val)) {
            val <- Lreject * runif(1, 1, 10)
            Lok <- FALSE
            message(paste("L(lsq, eta) NaN set to", val))
            return(val)
        }
        val <- as.double(val)
        return(val)
    }

    try({
        # 优化参数
        solveL <- optim(params, L, control = list( fnscale = -1, factr=factr ),
          method = "L-BFGS-B", lower = 1e-06)
        if (Lok) 
            #返回优化之后的参数
            params <- solveL$par
    })
    
    #优化之后的参数
    lsq <- params[1]
    eta <- params[-1]
    Sigma <- .C("GPcovar", as.double(X), as.integer(n), as.integer(d),
                as.double(params), Sigma = as.double(Sigma))$Sigma
    dim(Sigma) <- c(n, n)
    # 定义能否求逆
    inverseOK <- FALSE
    try({
        Sigma.inv <- ginv(Sigma)
        inverseOK <- TRUE
    })
    
    if (!request.functions) 
        # 探索过程到此返回，不再往下计算gp的预测值
        return(list(Sigma = Sigma, Sigma.inv = Sigma.inv, inverseOK = inverseOK, 
            X = X, y = y, params = params))
    c.scalar <- function(xi, xj) {
        h <- abs(xi - xj)
        val <- lsq * prod((1 + eta * h) * exp(-eta * h))
        return(val)
    }
    c.vec <- function(x) {
        vec <- rep(NA, dim(X)[1])
        for (j in 1:dim(X)[1]) vec[j] <- c.scalar(x, X[j, ])
        return(vec)
    }
    #定义GP对y的预测均值
    if (!inverseOK) 
        Ef <- function(x) NA
    else Ef <- function(x) {
        # GP对y×的预测均值
        Ey_x <- as.double(t(c.vec(x)) %*% Sigma.inv %*% ystar)
        # GP对y的预测均值
        return(mY + Ey_x)
    }
    if (!inverseOK) 
        sigmaf <- function(x) NA
    else sigmaf <- function(x) {
        c.vec_x <- c.vec(x)
        val <- lsq - as.double(t(c.vec_x) %*% Sigma.inv %*% c.vec_x)
        if (is.na(val)) 
            return(Inf)
        if (val < 0) {
            if (val > -0.001) 
                val <- -val
            else {
                rprt <- paste("sigmaf(")
                for (i in 1:length(x)) rprt <- paste(rprt, x[i])
                rprt <- paste(rprt, ") =", val)
                stop(rprt)
            }
        }
        return(sqrt(val))
    }
    return(list(Sigma = Sigma, Sigma.inv = Sigma.inv, inverseOK = inverseOK, 
        X = X, y = y, params = params, Ef = Ef, sigmaf = sigmaf))
    # Sigma：   协方差矩阵
    # params：  优化后的参数
    # Ef：      预测输出均值
    # sigmaf：  预测输出协方差矩阵
}

4.4.2 返回拟合结果

if(!GPfit$inverseOK)
	stop("GProcess inverse covariance matrix failed")
# 协方差的逆
Sigma.inv <- GPfit$Sigma.inv
# 传回已优化的参数
params <- GPfit$params
lsq <- params[1]
eta <- params[-1]

4.4.3 定义相关函数

sigmaf : 求 σ \sigma σ

c.scalar <- function(xi, xj) {
            h <- abs(xi - xj)
            val <- lsq * prod((1 + eta * h) * exp(-eta * h))
            return(val)
        }
c.vec <- function(x) {
            vec <- rep(NA, dim(X)[1])
            for (j in 1:dim(X)[1]) vec[j] <- c.scalar(x, X[j,])
            return(vec)
        }
sigmaf <- function(x) {
	c.vec_x <- c.vec(x)
    val <- lsq - as.double(t(c.vec_x) %*% Sigma.inv %*% c.vec_x)
    if (is.na(val)) 
    	return(Inf)
    if (val < 0) {
    	if (val > -0.001) 
       		val <- -val
        else {
        	xrep <- (x)
            rprt <- paste("sigmaf(")
            for (i in 1:length(xrep)) 
            	rprt <- paste(rprt,xrep[i])
            rprt <- paste(rprt, ") =", val)
            stop(rprt)
         }
    }
    return(sqrt(val))
}

4.4.4 链信息交换

nswaps : C n c h a i n s 2 C_{nchains}^2 Cnchains2​

# 链信息交换
if (nchains >= 2) {
	# 迭代次数：C_chains^2, 平均每2个链需要进行一次按照条件交换
	for (reps in 1:nswaps) 
		# 分解每个链
	  	for (index1 in sample(1:nchains)) 
        	for (index2 in sample((1:nchains)[-index1])) #任意选择2个链，判断是否交换
            {
            	if (T.mult == 1) 
                  	Tindex <- 1:nchains
                else 
                	Tindex <- sort(Temp, index.return = TRUE)$ix

                swap <- Tindex[c(index1, index2)]
                T1 <- Temp[swap[1]]# 链1的温度
                T2 <- Temp[swap[2]]# 链2的温度
                sdf1 <- sigmaf(x[swap[1]])# 链1的标准差
                sdf2 <- sigmaf(x[swap[2]])# 链1的标准差
                corf1 <- sdf1/sqrt(lsq) 
                corf2 <- sdf2/sqrt(lsq)
                fswap1 <- fx[swap[1]] + sdf1
                fswap2 <- fx[swap[2]] + sdf2
                alpha <- (fswap2/T1 - fswap1/T1 + fswap1/T2 - fswap2/T2)
                if(is.na(alpha)) alpha <- -Inf
                if (log(runif(1)) < alpha) {
                  Temp[swap[1]] <- T2
                  Temp[swap[2]] <- T1
                  swaps.record <- c(swaps.record, 1)
                }
                else {
                  swaps.record <- c(swaps.record, 0)
                }
       }
}

4.4.5 对每个链采样

Xnext <- X
ynext <- y
if (T.mult == 1) 
	Rank <- 1:nchains
else 
	Rank <- rank(Temp)
# 对每个链进行采样
for (chain in 1:nchains) {
	path <- Rank[chain]
    # 链对应的当前输入
    x.C <- as.double(x[chain, ])
    # 声明链对应的建议输入
    p.C <- NULL
    # 标准差
    sdf_x.C <- sigmaf(x.C)
    # 相关性
    corf_x.C <- sdf_x.C/sqrt(lsq)
    # 实际当前输出
    fx.C <- fx[chain]
    xp_Early <- NULL
    
    count <- 0
    repeat {
    	# 建议
        p.C <- rnorm(d)
        # 组合：原始+建议+0
        xp_Early <- matrix(c(x.C, p.C, rep(0, d)), ncol = 3)
        xOK <- FALSE
        try({
        	out <- .C("Leap", as.double(X), as.integer(dim(X)[1]),
                    as.integer(d), as.double(y), as.double(x.C),
                    as.double(p.C), as.integer(L), as.double(delta),
                    as.double(GPfit$Sigma.inv), as.double(params),
                    as.double(Temp[chain]), as.double(sqrt(lsq)*maxleap),
                    as.double(low), as.double(upp),
                    as.double(rep(0,n)), as.double(rep(0,n)),
                    as.double(rep(0,d)), as.double(rep(0,d)),
                    xp_Early = as.double(xp_Early) )
             xp_Early <- out$xp_Early
             dim(xp_Early) <- c(d,3)
             xOK <- TRUE
        })
        if (is.na(sum(xp_Early))) 
        	xOK <- FALSE
        if (xOK) 
        	break
        message("hybrid.explore: Proposed x failed")
        count <- count + 1
        if (count == 10) {
        	for (remv in 1:nchains) {
                    X <- X[-length(y), ]
                    y <- y[-length(y)]
            }
            GPfit <- GProcess(X, y, params = params,
                    request.functions = TRUE)
            if(!GPfit$inverseOK)
            	stop("GProcess inverse covariance matrix failed")
            Sigma.inv <- GPfit$Sigma.inv
            params <- GPfit$params
            lsq <- params[1]
            eta <- params[-1]
            sigmaf <- GPfit$sigmaf
            Xnext <- X
           	ynext <- y
            count <- 0
      	}
	}#endrepeat
    out <- list(x = xp_Early[, 1], p = xp_Early[, 2], Early = xp_Early[1, 3])
    x.P <- out$x
    OUTSIDE_BOUNDS <- FALSE
    if (BOUNDS) {
    	for (i in 1:d) {
        	if (x.P[i] < lb[i]) {
           		OUTSIDE_BOUNDS <- TRUE
           	}
           	if (x.P[i] > ub[i]) {
            	OUTSIDE_BOUNDS <- TRUE
           	}
  		}
   	}
    p.P <- out$p
    sdf_x.P <- sigmaf(x.P)
    corf_x.P <- sdf_x.P/sqrt(lsq)
    if(as.integer(count10/10) == count10/10) 
    	MAX10 <- 0
    if(corf_x.P > MAX10) 
    	MAX10 <- corf_x.P
    if (out$Early < 0) {
    	message(paste(sep="","Early stop with sd(f(x))/lambda = ",
                  round(corf_x.P,5), " where maxleap = ", maxleap))
    }
    else {
		if (corf_x.P > maxleap) 
        	stop(paste("Caution: Leapfrog gave sigmaf =", 
       	out$Early/sqrt(lsq), "when sigmaf(x)/lambda =", corf_x.P))
    }
    if (graph) {
     	clr <- clrs[path]
      	clr_red <- reds[path]
	    if(!OUTSIDE_BOUNDS)
    		lines(c(x.C[keyVars[1]], x.P[keyVars[1]]), 
                	c(x.C[keyVars[2]], x.P[keyVars[2]]), lwd = 0.2, 
               		 lty = 2, col = clr)
   	}
	if(OUTSIDE_BOUNDS) {
   		message(paste(sep="","T = ", Temp[chain],": Ignore point out of bounds"))
              fx.P <- alpha <- -Inf
    } 
    else {
    	date1 <- date()
        fx.P <- fn(x.P)
        date2 <- date()
        if (is.na(fx.P)) 
        	stop("f(x) NaN")
        if (fx.P == Inf) 
            stop("Infinite f(x) not allowed")
        fcalls <- fcalls + 1
        H_xp.C <- -fx.C - sdf_x.C + W(p.C)
        H_xp.P <- -fx.P - sdf_x.P + W(p.P)
        # 接受率
        alpha <- min(((-H_xp.P + H_xp.C)/Temp[chain]), 0)
   	}
   if (log(runif(1)) < alpha) {
    	x[chain, ] <- x.P
        fx[chain] <- fx.P
        Xnext <- rbind(Xnext, x.P)
        ynext <- c(ynext, fx.P)
        if(path <= Tinclude) {
        	Xfinal <- rbind(Xfinal, x.P)
        	yfinal <- c(yfinal, fx.P)
         	Taccept <- c(Taccept, path)
        	acceptance <- c(acceptance, 1)
            if(dim(Xfinal)[1] %% floor(n/10) == 0)
             	message(paste(sep="","hybrid.explore: ",
            round(100*dim(Xfinal)[1]/n),"% Complete with ",
          	round(100*mean(acceptance),1),"% acceptance of proposals"))
        }
       	if (graph) 
        	points(x.P[keyVars[1]], x.P[keyVars[2]], 
                  pch = 20, col = clr)
 	}
    #用X表示被拒绝的样本
    else {
              
       	if(path == 1) if(!OUTSIDE_BOUNDS) {
          	Xnext <- rbind(Xnext, x.P)
          	ynext <- c(ynext, fx.P)
         	Taccept <- c(Taccept, Tinclude+1)
      	}
      	if(path <= Tinclude) 
      		acceptance <- c(acceptance, 0)
        if (graph) if(!OUTSIDE_BOUNDS)
       		points(x.P[keyVars[1]], x.P[keyVars[2]], 
                  pch = "X", col = clr_red)
 	}
# 当样本数量足够n时，计算接受率
if(dim(Xnext)[1] > n) {
	i <- 0
  	repeat {
    	i <- i+1
   		if(Taccept[i] > Tinclude) {
        	Xnext <- Xnext[-i,]
           	ynext <- ynext[-i]
            i <- i-1
        }
       	if(length(ynext) <= n) 
       		break
	}
}
if (is.null(acceptance)) 
	mean.acceptance <- 0 
else
	mean.acceptance <- mean(acceptance[ceiling(length(acceptance)/2):length(acceptance)])
    if(!is.null(Xfinal)) if(dim(Xfinal)[1] >= n) 
    	break
    count10 <- count10+1
}

4.4.6 处理样本

# 返回所有样本
X <- Xnext
y <- ynext
if (is.null(acceptance)) 
	mean.acceptance <- NA
else 
	mean.acceptance <- mean(acceptance)
if (graph) {
 	if (refresh_count < nrefresh.graph) 
    	refresh_count <- refresh_count + 1
   	else {
   		refresh_count <- 1
  		key.indx <- key.indx + 1
    	if (key.indx > dim(KEY.VARS)[1]) 
        	key.indx <- 1
      	keyVars <- KEY.VARS[key.indx, ]
     	graph.template(X, key.vars = keyVars, lb, ub)
       	title(c("", "", "", "Exploratory phase"), adj = 0, 
                  cex.main = 0.8)
     	points(X[, keyVars[1]], X[, keyVars[2]], pch = 20, col="grey")
    	points(Xfinal[, keyVars[1]], Xfinal[, keyVars[2]], pch = 20)
    }
}
# 当样本数量达到n时，停止重复采样
if(!is.null(Xfinal)) if (dim(Xfinal)[1] >= n) break

4.5 获取探索阶段成果

  • 样本x，y
  • 训练出的 gp
  • 训练细节：
    • 链数：nchains
    • 温度数：T.mute
    • 交换次数 ：nswaps
    • L：似然函数
    • 密度函数：f
    • fcall：函数数
X <- Xfinal
y <- yfinal
GPfit <- GProcess(X, y, params, request.functions = TRUE, finetune=TRUE)
if(!GPfit$inverseOK)
	message("hybrid.explore: GProcess inverse covariance matrix failed")
details <- list(nchains = nchains, T.mult = T.mult, nswaps = nswaps, L = L, delta = delta,
        lb = lb, ub = ub, KEY.VARS = KEY.VARS)
return(list(X = X, y = y, f = fn, maxleap=maxleap, function.calls = fcalls, details = details, GPfit = GPfit))

5. sample 采样阶段

sample.out <- hybrid.sample(explore.out, n=500, graph=TRUE)

接下来分析hybrid.explore 函数

5.1 函数参数

  • Explore：探索阶段的成果
  • n：采样数
hybrid.sample <-
function (Explore, n = 1000, replace.target = c(0, 1, 2),
    lb=NULL, ub=NULL, L = NULL, delta = NULL, nchains = NULL,
    T.mult = NULL, maxleap=NULL, r = 5, nswaps = NULL, graph = FALSE)
{

5.2 准备阶段（初始化）

5.2.1 参数接收与校准

negInf <- function(X)
	as.double(apply(X, 2, function(x) 101*min(x)-100*max(x)))
posInf <- function(X)
	as.double(apply(X, 2, function(x) 101*max(x)-100*min(x)))
#获取探索阶段的样本
X <- Explore$X
y <- Explore$y
f <- Explore$f
d <- dim(X)[2]
if (d < 2)
	graph <- FALSE
nexplore <- dim(X)[1]
if (is.null(L)) L <- Explore$details$L
if (is.null(L)) L <- 1000
if (is.null(delta)) Explore$details$delta
if (is.null(delta)) delta <- 0.003
if (is.null(nchains)) Explore$details$nchains
if (is.null(nchains)) nchains <- 5
if (is.null(nswaps)) nswaps <- Explore$details$nswaps
if (is.null(nswaps)) nswaps <- choose(nchains, 2)
if (is.null(T.mult)) T.mult <- Explore$details$nchains
if (is.null(T.mult)) T.mult <- 1.5
if (T.mult == 1) 
	if (replace.target[1] == 1) 
    	replace.target <- 0
Temp <- T.mult^(1:nchains - 1)
if(is.null(lb)) lb <- Explore$details$lb
if(is.null(ub)) ub <- Explore$details$ub
BOUNDS <- FALSE
if ((!is.null(lb)) || (!is.null(ub))) 
	BOUNDS <- TRUE
if(is.null(lb)) 
	lb <- rep(-Inf, d)
low <- lb
low[!(low>-Inf)] <- negInf(X)[!(low>-Inf)]
if(is.null(ub)) 
	ub <- rep(Inf, d)
upp <- ub
upp[!(upp<Inf)] <- posInf(X)[!(upp<Inf)]
if(is.null(lb)) 
	lb <- rep(-Inf,d)
if(is.null(ub)) 
	ub <- rep(Inf,d)
if (BOUNDS) {
	if (!is.null(lb)) 
    	if (length(lb) == 1) 
       		lb <- rep(lb, d)
	if (!is.null(ub)) 
       	if (length(ub) == 1) 
          	ub <- rep(ub, d)
 	if ((length(lb) != d) || (length(lb) != d)) 
      	stop("lb, ub incorrect dimension")
   	if (is.na(sum(abs(lb)+abs(ub))))
     	 stop("NaN for lb, ub")
}
W <- function(p) return(sum(p^2)/2)
if (!replace.target[1] %in% 0:2) 
        stop(paste("hybrid.sample: Choose replace.target from:\n", 
            "1: Use true target distribution in primary chain only\n", 
            "0: Use true target distribution in all chains\n", 
            "2: Replace target by Gaussian process target approximation\n"))
# 样本索引
select <- 1:nexplore
# 去除违规样本
for(i in nexplore:1) 
	for(j in 1:d) {
		if ((X[i,j] < lb[j]) || (X[i,j] > ub[j]) )
        	select <- select[-i]
    }
# 过滤后的样本数
nselect <- length(select)
if(nselect < nchains) 
	stop("Not enough points supplied within bounds")
y_select <- y[select]
X_select <- X[select,]
indx_best <- sort(y_select, decreasing=TRUE, index.return=TRUE)$ix[1:nchains]
x <- X_select[indx_best, ]
fx <- y_select[indx_best]
if(is.null(dim(x))) 
	x <- t(x)
fcalls <- 0
x1 <- x[1, ]
fx1 <- fx[1]
###本质上这一块没用
digits <- rep(NA, d)
for (j in 1:d) {
	rangeX <- IQR(X[, j])
	if (rangeX == 0) 
		digits[j] <- 4
	else 
		digits[j] <- max(1, round(3 - log10(rangeX)))
}
rnd <- function(x) return(round(10^digits * x)/10^digits)
###本质上这一块没用
KEY.VARS <- NULL
key.vars <- 1:d
for (i in 1:(length(key.vars) - 1)) 
	for (j in (i + 1):length(key.vars)) {
        pair <- c(key.vars[i], key.vars[j])
        KEY.VARS <- rbind(KEY.VARS, pair)
    }
if (graph) {
	if(dim(X)[2] <= 2) 
		nrefresh.graph <- Inf 
	else
		nrefresh.graph <- ceiling(10/nchains)
    graph.template <- .hybrid.template
    count <- key.indx <- 1
    if (is.null(graph.template)) 
        graph.template <- Explore$details$graph.template
    keyvars <- KEY.VARS[key.indx, ]
    graph.template(X, key.vars = keyvars, lb, ub)
    title(c("", "", "", "Sampling phase"), adj = 0, cex.main = 0.8)
}
clrs <- c("blue", rep(c("green", "turquoise", "orange", "yellow", 
        "red", "grey", "pink"), nchains))
reds <- c("blue", rep(c("green", "turquoise", "orange", "yellow", 
        "red", "grey", "pink"), nchains))

# 温度设置
if (T.mult == 1) 
	Rank <- 1:nchains
else Rank <- rank(Temp)
	if (graph) {
        keyvars <- KEY.VARS[key.indx, ]
        points(X[, keyvars[1]], X[, keyvars[2]], col = "grey", 
            pch = 20)
    }
if (T.mult == 1) 
	Rank <-1:nchains
else 
	Rank <- rank(Temp)
for (chain in 1:nchains) {
    path <- Rank[chain]
    if (graph) 
    	points(x[chain, keyvars[1]], x[chain, keyvars[2]], 
                pch = 19, col = clrs[path])
}

5.2.2.高斯过程初始化

  • 检查并校准
if(is.null(Explore$GPfit))
	Explore$GPfit <- GProcess(X, y, params = NULL, request.functions = FALSE)
GPfit <- Explore$GPfit
if (!GPfit$inverseOK) 
	stop("hybrid.sample: Gaussian Process inverse covariance matrix failed")
params <- GPfit$params
Sigma.inv <- GPfit$Sigma.inv
ystar <- y - mean(y)
lsq <- params[1]
eta <- params[-1]
  • 准备
c.scalar <- function(xi, xj) {
	h <- abs(xi - xj)
	val <- lsq * prod((1 + eta * h) * exp(-eta * h))
	return(val)
}
c.vec <- function(x) {
	vec <- rep(NA, nexplore)
	for (j in 1:nexplore) vec[j] <- c.scalar(x, X[j, ])
	return(vec)
}
Ef <- function(x) {
	Ey_x <- as.double(t(c.vec(x)) %*% Sigma.inv %*% ystar)
	return(Ey_x)
}
sigmaf <- function(x) {
	c.vec_x <- c.vec(x)
	val <- lsq - as.double(t(c.vec_x) %*% Sigma.inv %*% c.vec_x)
	if (is.na(val)) {
		warning("hybrid.sample: NaN sigmaf(", rnd(x), ")\n")
       	return(Inf)
    }
	if (val < 0) {
    	if (abs(val) < 0.001) 
       		val <- abs(val)
       	else {
       		rprt <- paste("sigmaf(")
          	for (i in 1:length(x)) rprt <- paste(rprt, x[i])
           	rprt <- paste(rprt, ") =", val)
           	stop(rprt)
       	}
   	}
	return(sqrt(val))
}
#声明
Ey <- cory <- NULL
for(i in 1:nexplore) {
	Ey <- c(Ey, Ef(X[i,]))
   	cory <- c(cory, sigmaf(X[i,])/sqrt(lsq))
}
if (is.null(maxleap)) 
	maxleap <- max(cory)
SAMP <- ySAMP <- NULL
its <- 1
acceptance <- NULL
swaps.record <- swaps.record1 <- NULL

5.3重复采样

repeat { 内容如下  }

5.3.1 链信息交换

case1：

if (nchains >= 3) 
	for (reps in 1:nswaps) 
		for (index1 in sample(2:(nchains - 1))) {
   			index2 <- index1 + 1
       		if (T.mult == 1) 
          		Tindex <- 1:nchains
            else 
            	Tindex <- sort(Temp, index.return = TRUE)$ix
            swap <- Tindex[c(index1, index2)]
            T1 <- Temp[swap[1]]
            T2 <- Temp[swap[2]]
            if (replace.target[1] == 0) {
             	Efswap1 <- fx[swap[1]]
               	Efswap2 <- fx[swap[2]]
            }
           	else {
               	corf1 <- sigmaf(x[swap[1], ])/sqrt(lsq)
              	corf2 <- sigmaf(x[swap[1], ])/sqrt(lsq)
               	Efswap1 <- Ef(x[swap[1], ]) - r * (corf1 - maxleap)^2/(1 - corf1)^2
                Efswap2 <- Ef(x[swap[2], ]) - r * (corf2 - maxleap)^2/(1 - corf2)^2
            }
            alpha <- (Efswap2/T1 - Efswap1/T1 + Efswap1/T2 - Efswap2/T2)
            if(is.na(alpha)) 
            	alpha <- -Inf
            if (log(runif(1)) < alpha) {
             	Temp[swap[1]] <- T2
              	Temp[swap[2]] <- T1
           		swaps.record <- c(swaps.record, 1)
                }
             else {
             	swaps.record <- c(swaps.record, 0)
             }
 		}

case2：

if (nchains >= 2) {
 	if (T.mult == 1) 
      	Tindex <- 1:nchains
  	else Tindex <- sort(Temp, index.return = TRUE)$ix
       	swap <- Tindex[c(1, 2)]
    T1 <- Temp[swap[1]]
    if (T1 != 1) 
    	stop(paste("hybrid.sample swaps: T1 =", T1))
  	T2 <- Temp[swap[2]]
   	if (replace.target[1] == 0) {
   		fswap1 <- Efswap1 <- fx[swap[1]]
		fswap2 <- Efswap2 <- fx[swap[2]]
	}
   	else if (replace.target[1] == 2) {
     	corf1 <- sigmaf(x[swap[1], ])/sqrt(lsq)
    	corf2 <- sigmaf(x[swap[1], ])/sqrt(lsq)
      	Efswap1 <- fswap1 <- Ef(x[swap[1], ]) - r * (corf1 - maxleap)^2/(1 - corf1)^2
        Efswap2 <- fswap2 <- Ef(x[swap[2], ]) - r * (corf2 - maxleap)^2/(1 - corf2)^2
    }
   	else if (replace.target[1] == 1) {
 		fswap1 <- fx1
      	xswap2 <- x[swap[2], ]
        corf1 <- sigmaf(x[swap[1], ])/sqrt(lsq)
    	corf2 <- sigmaf(x[swap[1], ])/sqrt(lsq)
        Efswap1 <- Ef(x[swap[1], ]) - r * (corf1 - maxleap)^2/(1 - corf1)^2
      	Efswap2 <- Ef(x[swap[2], ]) - r * (corf2 - maxleap)^2/(1 - corf2)^2
       	date1 <- date()
      	fswap2 <- f(x[swap[2], ])
       	date2 <- date()
      	mean.acceptance <- 0
        if (!is.null(acceptance))
          	mean.acceptance <- mean(acceptance)
        if (is.na(fswap2)) 
	        stop("f(x) NaN")
        if (fswap2 == Inf) 
           	stop(paste("Swaps: Infinite f(x) not allowed:\n", 
                    "f(", deparse(xswap2), ") =", fswap2))
       	fcalls <- fcalls + 1
     	if (graph) 
            .show.time(date1, date2)
 	}
  	alpha <- ((fswap2 - fswap1)/1 + (Efswap1 - Efswap2)/T2)
    if(is.na(alpha)) alpha <- -Inf
    if (log(runif(1)) < alpha) {
     	Temp[swap[1]] <- T2
       	Temp[swap[2]] <- 1
       	swaps.record1 <- c(swaps.record1, 1)
       	if (replace.target[1] == 1) {
         	x1 <- xswap2
          	fx1 <- fswap2
       	}
 	}
	else
		swaps.record1 <- c(swaps.record1, 0)
}

5.3.2 重排温度

if (T.mult == 1) 
	Rank <- 1:nchains
else 
	Rank <- rank(Temp)

5.3.3 对每个链采样

for (chain in 1:nchains) {
	path <- Rank[chain]
	x.C <- x[chain, ]
    p.C <- NULL
   	if (replace.target[1] == 0) 
  		fx.C <- fx[chain]
   	else if ((replace.target[1] == 1) && (path == 1)) 
        fx.C <- fx1
    else 
    	fx.C <- Ef(x[chain, ])
    xp_Early <- NULL
    repeat {
    	p.C <- rnorm(d)
       	xp_Early <- matrix(c(x.C, p.C, rep(0, d)), ncol = 3)
      	xOK <- FALSE
       	delta_adj <- delta
      	try({
           	out <- .C("Leap", as.double(X), as.integer(dim(X)[1]),
	            as.integer(d), as.double(y), as.double(x.C),
	            as.double(p.C), as.integer(L), as.double(delta),
	            as.double(GPfit$Sigma.inv), as.double(params),
	            as.double(Temp[chain]), as.double(sqrt(lsq)*maxleap),
	            as.double(low), as.double(upp),
	            as.double(rep(0,n)), as.double(rep(0,n)),
	            as.double(rep(0,d)), as.double(rep(0,d)),
            xp_Early = as.double(xp_Early) )
        	xp_Early <- out$xp_Early
            dim(xp_Early) <- c(d,3)
            xOK <- TRUE
           	xOK <- TRUE
   		})
        xout <- xp_Early[, 1]
        if (is.na(sum(xp_Early))) 
         	xOK <- FALSE
        if (xOK) 
           	break
	}
   	out <- list(x = xp_Early[, 1], p = xp_Early[, 2])
   	H_xp.C <- -fx.C + W(p.C)
	x.P <- out$x
   	p.P <- out$p
  	if (graph) 
  		lines(c(x.C[keyvars[1]], x.P[keyvars[1]]), 
            	c(x.C[keyvars[2]], x.P[keyvars[2]]), lwd = 0.5, 
                lty = 2, col = clrs[path])
   	if (replace.target[1] == 2) 
    	fx.P <- Ef(x.P)
	else if ((replace.target[1] == 1) && (path > 1)) 
     	fx.P <- Ef(x.P)
	else {
     	date1 <- date()
      	fx.P <- f(x.P)
       	date2 <- date()
        if (is.na(fx.P)) 
     		stop("f(x) NaN")
        if (fx.P == Inf) 
            stop(paste("Leapfrog: Infinite f(x) not allowed:\n", 
                  "f(", deparse(x.P), ") =", fx.P))
      	fcalls <- fcalls + 1
        if (graph) 
        if (path == 1) 
      		.show.time(date1, date2)
    }
  	H_xp.P <- -fx.P + W(p.P)
   	if ((replace.target[1] == 2) || ((replace.target == 1) && (path > 1))) {
       	corf_x.C <- sigmaf(x.C)/sqrt(lsq)
        corf_x.P <- sigmaf(x.P)/sqrt(lsq)
        H_xp.C <- H_xp.C + r * (corf_x.C - maxleap)^2/(1 - corf_x.C)^2
        H_xp.P <- H_xp.P + r * (corf_x.P - maxleap)^2/(1 - corf_x.P)^2
   	}
    alpha <- min(((-H_xp.P + H_xp.C)/Temp[chain]), 0)
    if(is.na(alpha)) alpha <- -Inf
   	if (log(runif(1)) < alpha) {
      	x[chain, ] <- x.P
   		if (replace.target[1] == 0) 
      		fx[chain] <- fx.P
   		else if (replace.target[1] == 1) 
   			if (path == 1) 
   				fx1 <- fx.P
		if (path == 1) {
       		SAMP <- rbind(SAMP, x.P)
   			ySAMP <- c(ySAMP, fx.P)
       		acceptance <- c(acceptance, 1)
     		if(dim(SAMP)[1] %% floor(n/10) == 0)
       		message(paste(sep="","hybrid.sample: ",
		              round(100*dim(SAMP)[1]/n),"% Complete with ",
		              round(100*mean(acceptance),1),"% acceptance of proposals"))
    }
	if (graph) 
	  	points(x.P[keyvars[1]], x.P[keyvars[2]], 
    				pch = 20, col = clrs[path])
   	}
   	else {
		if (graph)
       		points(x.P[keyvars[1]], x.P[keyvars[2]], 
                  pch = "X", col = reds[path])
      	if (path == 1) {
        	SAMP <- rbind(SAMP, x.C)
           	ySAMP <- c(ySAMP, fx.C)
         	acceptance <- c(acceptance, 0)
       	}
 	}
	mean.acceptance <- 0
  	if (!is.null(acceptance))
 		mean.acceptance <- mean(acceptance)
}

5.3.4 样本处理与检查

当采样数量一定时，退出循环

if (graph) {
 	if (count >= nrefresh.graph) {
    	count <- 1
   		key.indx <- key.indx + 1
  		if (key.indx > dim(KEY.VARS)[1]) 
       		key.indx <- 1
      	keyvars <- KEY.VARS[key.indx, ]
        graph.template(X, key.vars = keyvars, lb, ub)
        points(X[, keyvars[1]], X[, keyvars[2]], col = "grey", pch = 20)
        points(x.P[keyvars[1]], x.P[keyvars[2]], pch = 20, col = "grey")
       	title(c("", "", "", "Sampling phase"), adj = 0, cex.main = 0.8)
  	}
	points(x[, keyvars[1]], x[, keyvars[2]], pch = 20, col = "grey")
    points(SAMP[acceptance>0, keyvars[1]], SAMP[acceptance>0, keyvars[2]], pch = 20, col = "black")
    count <- count + 1
}
its <- its + 1
if (length(acceptance) > n) 
	break

5.4 抽样结束

当5.4样本检查跳出循环后，返回采样结果：

return(list(SAMP = SAMP, y = ySAMP, acceptance = acceptance, function.calls = fcalls))

6.画图

opar <- par(mfrow=c(2,1))
plot(density(sample.out$SAMP[,1]), xlab="x1", ylab="f(x)")
plot(density(sample.out$SAMP[,2]), xlab="x2", ylab="f(x)")
par(opar)

7.实验结果

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

风尘23187
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
专栏目录
23825736gp.rar_GP matlab_GP算法_gp_matlab GP_new chaos
07-14
【标题】"23825736gp.rar_GP matlab_GP算法_gp_matlab GP_new chaos" 提供的资源是一个MATLAB实现的GP（Genetic Programming）算法，它是一种基于进化计算的优化方法，用于解决各种复杂问题。"new chaos"可能指的是...
R语言--MCMC算法介绍以及例子
qq_45638191的博客
01-02 2614
MCMC算法介绍以及例子
MCMC_Gibbs/MH采样的简单案例R代码实现
qq_43331366的博客
10-29 6479
文章目录练习：用程序实现正态分布均值、方差的后验分布抽样。题目背景Gibbs抽样（Bayes公式推导）R代码实现μ\muμ的更新函数σ2\sigma^2σ2的更新函数主函数案例求解（补充）RCPP代码 练习：用程序实现正态分布均值、方差的后验分布抽样。 题目背景 记Xi,i=1...,n{X_i,i=1...,n}Xi​,i=1...,n来自一维正态N(μ,σ2)N(\mu,\sigma^2)N(μ,σ2)的独立同分布样本，μ\muμ具有先验分布N(a,b)N(a,b)N(a,b)，σ2\sigma^2σ
gpxe.rar_GPXE_gp_linux gp_pxe_无盘
09-23
总之，"gpxe.rar_GPXE_gp_linux gp_pxe_无盘"这个主题涉及到的是使用GPXE无盘启动技术在Linux环境下的应用，这包括了无盘工作站的设置、GPXE源码的分析以及可能的PXE协议增强。通过这种技术，可以实现高效、集中化的...
GP.zip_GP python_gaussian process_gp_zip
09-19
标题中的"GP.zip_GP python_gaussian process_gp_zip"表明这是一个关于使用Python实现高斯过程（Gaussian Process，简称GP）的压缩包文件。高斯过程是一种统计模型，常用于机器学习中的回归和分类任务，特别是在处理...
GP.rar_GP matlab_Guassian Process_gp_regression
09-20
【标题】"GP.rar_GP matlab_Guassian Process_gp_regression" 涉及的主要知识点是高斯过程（Gaussian Process, GP）在机器学习领域的应用，特别是针对分类和回归问题的解决方案。高斯过程是一种非参数概率模型，常...
GP_baozi9.rar_G-P_gp_gp algorithm_gp 算法_关联维
07-13
Matlab是一种广泛使用的数值计算和编程环境，而mex文件是Matlab与C、C++或Fortran等编译语言之间的接口，它可以提高代码运行效率。这里的mex版G-P算法程序意味着用户可以利用Matlab的便捷性，同时享受到接近原生编译...
略谈set与map的pair封装与进入哈希
creative_mind的博客
08-09 712
一，哈希是一种数学函数，它接受一个输入（或“消息”），然后返回一个通常更小的固定大小的输出，这个输出称为“哈希值”或“哈希码”。映射：哈希思想通过哈希函数将任意长度的数据映射到固定长度的哈希值。这个映射过程是单向的，即从数据到哈希值是容易的，但从哈希值回溯到原始数据几乎是不可能的。快速性：哈希函数的设计旨在快速计算，以便在大数据集中实现高效的数据访问。均匀分布：理想情况下，哈希函数应该能够将输入数据均匀地分布在哈希值空间中，以减少冲突并提高查找效率。
基于R语言APSIM模型进阶应用与参数优化、批量模拟技术教程
2301_78164062的博客
08-09 399
基于R语言APSIM模型进阶应用与参数优化、批量模拟技术教程
脱胎于 S 语言的R语言，Ross Ihaka 和 Robert Gentleman 和社区的力量让 R 在学术界与研究机构放光彩
2301_79425796的博客
08-09 369
R语言从一门用于统计学教学的编程语言，发展成为全球数据科学领域的重要工具，离不开其强大的功能、丰富的社区资源和开源精神。这些都离不开Ross Ihaka 和 Robert Gentleman 和 社区的力量。
MATLAB算法实战应用案例精讲-【数模应用】均值z 检验（附R语言、python和MATLAB代码实现）
qq_36130719的博客
08-09 36
Z检验（Z Test）又叫U检验。由于实际问题中大多数随机变量服从或近似服从正态分布，U作为检验统计量与X的均值是等价的，且计算U的分位数或查相应的分布表比较方便。通过比较由样本观测值得到的U的观测值，可以判断数学期望的显著性，我们把这种利用服从标准正态分布统计量的检验方法称为U检验（U-test）。z检验是指一种常用的统计学方法，用于检验样本均值与总体均值之间是否存在显著差异。z检验适用于样本数据服从正态分布的情况下，通过计算z值来判断样本均值是否显著偏离总体均值。
R语言医疗数据分析笔记
LJWWD的博客
08-06 1026
分组因子：用于标识和分离不同组的数据，用于统计分析（如性别、地区等）。数组：用于存储和处理数值数据的结构，可以用来分析或计算。分组因子一般是分类变量，帮助将数据分组。连续变量通常是在不同分组之间进行比较的对象，以分析差异和趋势。x是一个分组的数据列表，包含两个组的观测值（如 Age）。y是展平后的向量，包含了所有观测值。g是一个因子，表示每个观测值所属的组，用于计算 p 值的分组比较。详细分析和这两行代码。我们将通过实际的示例数据来说明它们的作用和含义。t.test。
基于R语言绘制GGE双标图3
maizeman126的博客
08-07 1365
把各环境产量均值第一的品种G6与高产与稳产品种G9进行比较，过双标图原点且与2个品种连线相垂直的一条直线将图分割为2个部分，E6、E3和E5位于分割线的下部分，意味着G9在这2个试验点的产量表现优于G6，特别是E6(抗旱鉴定环境)试验点。GGE 双标图可用于直接比较两个品种在各个试验点的表现，将要比较的两个品种用直线连接起来，过双标图圆点与两个品种连线相垂直的一条直线将图分割为2个部分。本图结论：品种G9在同心圆的内圆中，说明它们是青稞试验区的高产稳产品种，而G3距同心圆最远，是丰产性、稳产性最差的品种。
基于R语言绘制GGE双标图2
maizeman126的博客
08-06 597
在这个前提下,去掉一些在所有或多数环境下都低产的品种，可以减少双标图中G的分量，使较多的GE得以在双标图上表现，从而取得对环境更好的区分。SD-定标的GGE 双标图与多元分析中常用的主成分分析模型相一致， 因而在多元数据(如品种-性状两向表, 品种-遗传标记两向表, 等等)分析中有广泛应用。按照2-D双标图表现出的环境间关系把环境分组,然后按组分别研究,是解决数据较大、关系较复杂的有效途径。，故这里只讨论SD标定。在对双标图的解释中，有一个隐含的假设，就是所讨论的 2-D双标图可以充分近似所代表的两。
仿RabbiteMq简易消息队列基础篇（gtest的使用）
最新发布
m0_62812354的博客
08-10 719
gtest是google的一个开源框架，它主要用于写单元测试，检查自己的程序是否符合预期行为。可在多个平台上使用（包含Linux，MAC OC，Windows等）。它提供了丰富的断言，致命和非致命失败判断，能进行值参数化测试，类型参数化测试，“死亡测试”。
1章3节：R 语言的产生与发展轨迹
2301_79425796的博客
08-06 789
药物的研究阶段包括靶标确定、生物学模型建立、新化学实体的发现和优化新化学实体四个关键环节。每个环节都在药物研发过程中发挥着至关重要的作用。通过科学严谨的方法和技术手段，研究人员能够在药物研发过程中不断探索和优化，从而发现和开发出具有高效、安全性的候选药物，为疾病治疗提供更多的选择和希望。
SO_REUSEADDR 和 SO_REUSEPORT 的区别
u013669912的博客
08-09 794
…
代码+视频，R语言VRPM绘制多种模型的彩色列线图
dege857的博客
08-09 129
列线图，又称诺莫图（Nomogram），它是建立在回归分析的基础上，使用多个临床指标或者生物属性，然后采用带有分数高低的线段，从而达到设置的目的：基于多个变量的值预测一定的临床结局或者某类事件发生的概率。列线图(Nomogram)可以用于多指标联合诊断或预测疾病发病或进展。咱们既往已经多篇文章介绍绘制列线图，今天咱们来视频介绍一下VRPM包绘制彩色列线图，这个包可以绘制多个模型的列线图，咱们一一来介绍。近些年来在高质量SCI临床论文中用的越来越多。R语言VRPM绘制多种模型的彩色列线图。
R 语言学习教程，从入门到精通，R 基础运算（5）
qq_45746668的博客
08-05 463
这个赋值符号是 R 语言的一个形式上的优点和操作上的缺点：形式上更适合数学工作者，毕竟不是所有的数学工作者都习惯于使用 = 作为赋值符号。 操作上来讲，< 符号和 - 符号都不是很好打的字符，这会让很多程序员不适应。因此，R 语言的比较新的版本也支持 = 作为赋值符
OpenCASCADE基础教程：gp_Pnt与gp_Vec解析
"这篇内容主要介绍了Open CASCADE技术的基础知识，特别是gp_Pnt和gp_Vec类的使用，以及如何判断空间点的位置关系。" Open CASCADE（简称OCC）是一个开源的三维几何建模库，广泛应用于CAD、CAE和CAM等领域。它提供了...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

风尘23187

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值