在R语言中从非标准函数生成样本的技术问询

作为刚接触R的新手，碰到这种没有现成采样函数的自定义分布确实容易卡壳，我来一步步帮你解决问题——先从修正你写的函数开始，再逐个分析你尝试过的方法，给出可直接运行的代码。

一、先把目标密度函数写对 #

你提供的Falp函数有几个明显的问题：没传入关键参数x和N、定义了没用的变量ti、初始化矩阵时依赖未定义的x。先把基础的未归一化目标密度函数写好（采样时大多不需要归一化常数，省掉计算麻烦）：

# 未归一化的目标密度函数，x必须>0，N为样本数量
target_density <- function(x, N) {
  if (any(x <= 0)) return(0)  # 过滤x<=0的无效输入
  1 / (gamma(1 + 1/x))^N
}

# 测试一下：x=1时gamma(2)=1，所以结果为1；x=2时gamma(1.5)=sqrt(pi)/2≈0.886，结果≈1/(0.886^10)≈0.63
target_density(x = c(1, 2), N = 10)

二、为什么分位数法行不通？ #

分位数法需要知道累积分布函数（CDF）的反函数，但这个分布的CDF是$\int_0^x f(t)dt$，这个积分没有解析解，根本没法写出反函数，所以这个方法直接放弃就好，不是你的问题~

三、拒绝采样的简便实现：用自适应拒绝采样包 #

手动找合适的提议分布和常数M对新手太麻烦，R里的ars包专门处理单峰分布的自适应拒绝采样，自动帮你搞定这些细节：

1. 先安装并加载包：

install.packages("ars")
library(ars)

2. 定义对数形式的目标密度（ars包要求输入对数密度，数值稳定性更好）：

log_target <- function(x) {
  if (x <= 0) return(-Inf)  # x<=0时密度为0，对应对数负无穷
  -N * log(gamma(1 + 1/x))
}

3. 生成样本并验证：

N <- 10  # 假设样本数量为10
# 生成1000个样本，初始点选两个能覆盖密度主要区域的点（比如0.5和2）
samples_ars <- ars(n = 1000, logf = log_target, x = c(0.5, 2))

# 对比采样结果和目标密度的形状
x_vals <- seq(0.1, 5, length.out = 1000)
dens_vals <- sapply(x_vals, target_density, N = N)

plot(density(samples_ars), main = "拒绝采样样本密度 vs 目标密度", col = "red", lwd=2)
# 把目标密度缩放后对齐绘制
lines(x_vals, dens_vals / max(dens_vals) * max(density(samples_ars)$y), col = "blue", lwd=2)
legend("topright", legend = c("采样样本密度", "目标密度"), col = c("red", "blue"), lty = 1)

四、Metropolis-Hastings(MH)算法的手动实现 #

如果不想依赖第三方包，自己写MH算法也不难。这里用对数变换来避免处理x<=0的问题（把x转成y=log(x)，y的范围是全体实数，提议分布用正态分布更方便）：

mh_sample <- function(n_samples, N, sigma = 0.5, burn_in = 2000) {
  # 初始化：从x=1开始，对应y=log(1)=0
  y_current <- 0
  samples <- numeric(n_samples + burn_in)
  
  for (i in 1:(n_samples + burn_in)) {
    # 生成候选y值，用正态提议分布
    y_proposal <- rnorm(1, mean = y_current, sd = sigma)
    x_proposal <- exp(y_proposal)
    x_current <- exp(y_current)
    
    # 计算接受率（因为提议分布对称，雅可比行列式抵消，直接用对数密度差）
    log_alpha <- log_target(x_proposal) - log_target(x_current)
    alpha <- min(1, exp(log_alpha))
    
    # 决定是否接受候选样本
    if (runif(1) <= alpha) {
      y_current <- y_proposal
    }
    
    samples[i] <- exp(y_current)
  }
  
  # 去掉前2000个燃烧期样本，返回有效样本
  samples[(burn_in + 1):(burn_in + n_samples)]
}

# 生成1000个样本
samples_mh <- mh_sample(n_samples = 1000, N = 10, sigma = 0.5)

# 验证采样结果
plot(density(samples_mh), main = "MH算法样本密度 vs 目标密度", col = "green", lwd=2)
lines(x_vals, dens_vals / max(dens_vals) * max(density(samples_mh)$y), col = "blue", lwd=2)
legend("topright", legend = c("MH样本密度", "目标密度"), col = c("green", "blue"), lty = 1)

注意：sigma是提议分布的标准差，需要调整——如果接受率太低（比如<0.2）就调小sigma，太高（比如>0.5）就调大sigma，一般维持在0.2-0.5之间采样效率最高。

五、你的原函数Falp修正版 #

如果只是想计算给定x和N的密度值，修正后的函数可以写成这样：

Falp <- function(x, N) { 
  if (any(x <= 0)) stop("x必须大于0哦！")
  1 / ((gamma(1 + 1/x))^N)
}

# 测试
Falp(x = c(0.5, 3), N = 5)

内容的提问来源于stack exchange，提问作者miki