网站跳转怎么做360,做购物网站需要多少钱,教育政务网站建设,网站建设项目明细表背景介绍#xff1a;#xff08;了解采样的可以跳过#xff09;1#xff09;为什么需要采样#xff1a;简单的分布#xff0c;比如高斯、exponential、gamma等等的样本都可以直接用numpy.random生成#xff0c;但复杂的分布需要采样器生成。在贝叶斯、概率编程里面…背景介绍了解采样的可以跳过1为什么需要采样简单的分布比如高斯、exponential、gamma等等的样本都可以直接用numpy.random生成但复杂的分布需要采样器生成。在贝叶斯、概率编程里面有很多复杂的分布而贝叶斯更新需要对这些复杂分布再采样。2最简单的3种采样CDF采样把pdf转成cdf然后对0-1区间均匀采样拒绝采样蓝线是真实分布红线是一个处处y坐标值大于蓝线的函数正比于一个简单分布比如均匀分布或者高斯分布或者一个分段函数对红线的分布采样得到样本x然后计算蓝线(x)/红线(x)作为接受这个样本的概率。重要性采样类似拒绝采样不需要红线函数处处大于蓝线但每个样本x的权重是蓝线(x)/红线(x)3以上三种采样的问题cdf采样需要一个数值积分所以在高维空间采样的时候计算量太大拒绝采样需要先构造一个总是值大于自己的函数重要性采样在对分布情况不太清楚的时候采样方差可能很大即大部分样本的权重都很低都或多或少需要对分布本身有一定了解在高维空间不够高效4Metropolis-Hastings特点用一个多元正态分布随机游走不需要对分布本身有太多了解高维空间时有效分布本身不需要积分1基本只需要保证f(x)在-∞和∞处00f(x)upper bound就行了。方法初始点随机游走如果 则接受这个新的样本否则按照 的概率接受这个采样。问题随机游走的效率仍然不够高。代码def binormal_draw(xprev,beta1):mean [0, 0]cov [[beta**2,0],[0,(1.5*beta)**2]]binormal scipy.stats.multivariate_normal(mean,cov)return xprev binormal.rvs()def metropolis(F, qdraw, nsamp, x_init, burnin, thinning2, beta1):samplesnp.empty((nsamp,2))x_prev x_initaccepted 0j 0for i in range((nsampburnin)*thinning):x_star qdraw(x_prev, beta)logp_star np.log(F(x_star[0], x_star[1]))logp_prev np.log(F(x_prev[0], x_prev[1]))logpdfratio_p logp_star-logp_prevu np.random.uniform()if np.log(u) logpdfratio_p:x_prev x_starif i burnin*thinning and i%thinning 0:samples[j] x_starj 1accepted 1else:#we always get a sampleif i burnin*thinning and i%thinning 0:samples[j] x_prevj 1return samples, accepted可视化http://elevanth.org/blog/2017/11/28/build-a-better-markov-chain/elevanth.orgHMC1直观理解原理把MCMC的样本随机游走改成在一个势场中具有动能的质点势场由分布函数f(x)决定初始动能随机给定这个质点在运动时间t后的位置记录下来并且作为下次采样的初始点。2相比MH的优点更新样本点的时候使用了梯度所以对复杂分布采样比随机试的速度要快。2公式、可视化参考下列文章Markov Chains: Why Walk When You Can Flow?elevanth.orgXinyu Chen如何简单地理解「哈密尔顿蒙特卡洛 (HMC)」zhuanlan.zhihu.comhttps://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/theclevermachine.wordpress.comMCMC: Hamiltonian Monte Carlo (a.k.a. Hybrid Monte Carlo)theclevermachine.wordpress.comhttps://arxiv.org/pdf/1701.02434.pdfarxiv.orgHamiltonian Monte Carlo explained by Alex Rogozhnikovarogozhnikov.github.io3代码def HMC(F, u0, n_iter, N_iter, h0.01): # part 1: 初始化一个orbit的变量第一行是初始点u0的坐标orbit torch.zeros((n_iter1, 2))orbit[0] u0.detach()u orbit[0].unsqueeze(0)shape u.size()# part 2: 循环n_iter次每次初始化一个v0和u0然后让leapfrog走N_iter次每次时长h默认0.01for k in tqdm(range(n_iter)):v0 torch.randn(sizeu0.shape)u0 torch.randn(sizeu0.shape)*3u, v leapfrog(F, u0, v0, h, N_iter,shape)u0 orbit[k]a float(ratio(F, u0, v0, u, v, shape))r np.random.rand()if r a:orbit[k1] uelse:orbit[k1] u0return orbit #[10:, :]# part 3: leapfrog算法此外还有velocity verlet等等都可以尝试但不要用euler
def leapfrog(F, u, v, h, N_iter,shape):v v - h/2 * grad(F, u)for i in range(N_iter-1):u u h * vv v - h * grad(F, u)u u h * vv v - h/2 * grad(F, u)return u, v# part 4: 对函数F求梯度
def grad(F, u):u u.detach()if u.requires_grad False:u u.requires_grad_()output -torch.log(F(u))output.backward()ugrad u.grad.squeeze(0)u u.squeeze(0)return ugrad# part 5: 细节处理
def unsqueeze(tensor, shape):if len(list(tensor.size())) len(list(shape)):tensor tensor.unsqueeze(0)return tensor# part 6: 由于leapfrog计算结果仍然不是绝对的稳定所以接受概率是min{运行前后总能量的比值,1}
#大多数情况下都≈1
def ratio(F, u0, v0, u, v, shape):u0 unsqueeze(u0, shape)u unsqueeze(u, shape)v0 unsqueeze(v0, shape)v unsqueeze(v, shape)w0 - torch.log(F(u0)) 0.5*torch.mm(v0, torch.t(v0))w1 - torch.log(F(u)) 0.5*torch.mm(v, torch.t(v))return torch.exp(w0-w1)4调参细节一共有5组主要的超参要调v0的标准差u0的位置和标准差时间步长h总步长数N_iterburning。5效果此处暂时无图。。不过采样的速度总体比Metropolis Hastings要稳定样本质量较高使用了梯度所以效率高。6进阶NUTS采样The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carloarxiv.org
