接受拒绝采样matlab代码

接受拒绝采样是一种概率分布采样方法，常用于生成服从非常规分布的随机变量。这种方法的基本思想是：首先生成一个服从某种比较容易处理的分布的随机变量，然后按照特定的条件接受或拒绝这个随机变量，最终得到所需的随机变量。Matlab代码实现如下：

function x = acceptReject(N, p, f, f_max)
%N为采样数，p为目标分布，f为仿真分布，f_max为仿真分布最大值
x = zeros(N,1); %初始化结果向量
for i = 1:N
accept = false; %初始化接受标志
while ~accept
x_prop = rand*f_max; %生成随机变量
u = rand; %生成接受概率
if u <= p(x_prop)/(f(x_prop)*f_max) %按照条件计算接受概率
x(i) = x_prop; %接受样本
accept = true; %修改接受标志
end
end
end

其中，N表示采样数，p为目标分布，f为仿真分布，f_max为仿真分布的最大值。首先，初始化结果向量x，然后对每个采样数执行循环。在循环内部，生成一个在仿真分布中的随机变量x_prop和一个介于0和1之间的接受概率u。如果接受概率满足p(x_prop)/(f(x_prop)*f_max)，则接受该样本，并将接受标志设为真。最终，返回结果向量x。

可以使用该函数生成服从各种非常规分布的随机变量，如beta分布、伽马分布、指数分布等。

相关问题

M-H采样代码 MATLAB实现

### 回答1：
M-H采样是一种基于马尔科夫链蒙特卡罗方法的采样算法，可以用于从某个目标分布中采样。下面是MATLAB实现的M-H采样代码：

function [samples, accept_rate] = metropolis_hastings(target_dist, proposal_dist, init_state, num_samples, burn_in)
% target_dist: 目标分布
% proposal_dist: 提议分布
% init_state: 初始状态
% num_samples: 采样数
% burn_in: 热身期

samples = zeros(1, num_samples);
accept_rate = 0;
state = init_state;
for i = 1 : burn_in + num_samples
    % 从提议分布中采样得到候选状态
    candidate_state = proposal_dist(state);
    % 计算接受率
    acceptance_prob = min(1, target_dist(candidate_state) / target_dist(state) * proposal_dist(state, candidate_state) / proposal_dist(candidate_state, state));
    % 决定是否接受候选状态
    if rand() < acceptance_prob
        state = candidate_state;
        if i > burn_in
            samples(i - burn_in) = state;
        end
        accept_rate = accept_rate + 1;
    else
        if i > burn_in
            samples(i - burn_in) = state;
        end
    end
end
accept_rate = accept_rate / num_samples;
end

其中，`target_dist`表示目标分布的概率密度函数，`proposal_dist`表示提议分布的概率密度函数；`init_state`表示初始状态，`num_samples`表示采样数，`burn_in`表示热身期。`proposal_dist`可以是对称分布，如均匀分布或正态分布，也可以是非对称分布，如指数分布或Gamma分布。

使用时，需要自己定义目标分布和提议分布，并对其进行封装成函数。例如，假设我们要从一个标准正态分布中采样，使用均值为当前状态、方差为1的正态分布作为提议分布，则可以这样调用函数：

% 目标分布：标准正态分布
target_dist = @(x) exp(-x^2/2) / sqrt(2 * pi);
% 提议分布：均值为当前状态、方差为1的正态分布
proposal_dist = @(x) normrnd(x, 1);
% 初始状态：0
init_state = 0;
% 采样数：1000
num_samples = 1000;
% 热身期：100
burn_in = 100;

% 进行M-H采样
[samples, accept_rate] = metropolis_hastings(target_dist, proposal_dist, init_state, num_samples, burn_in);

% 输出采样结果和接受率
disp(['采样结果：', num2str(samples)]);
disp(['接受率：', num2str(accept_rate)]);

注：这里的代码仅作示例，实际应用时需要注意提议分布的选择和参数的调整，以得到较高的采样效率。

### 回答2：
M-H采样是一种用于从目标分布中采样的马尔可夫链蒙特卡罗方法。这种方法在统计学和机器学习中得到广泛应用。下面是一段用MATLAB实现的M-H采样代码：

function [samples, acceptance_ratio] = metropolis_hastings(target_distribution, initial_sample, num_samples)
    % 初始化采样样本
    samples = zeros(1, num_samples);
    samples(1) = initial_sample;
    
    % 初始化接受率
    acceptance_count = 0;
    
    for i = 2:num_samples
        % 提议新的样本
        proposal = normrnd(samples(i-1), 1);
        
        % 计算接受概率
        accept_prob = min(1, target_distribution(proposal) / target_distribution(samples(i-1)));
        
        % 接受或拒绝样本
        if rand() < accept_prob
            samples(i) = proposal;
            acceptance_count = acceptance_count + 1;
        else
            samples(i) = samples(i-1);
        end
    end
    
    % 计算接受率
    acceptance_ratio = acceptance_count / (num_samples - 1);
end

上述代码中，`target_distribution`是目标分布的概率密度函数，`initial_sample`是初始样本，`num_samples`是需要采样的样本总数。代码中使用了正态分布作为提议分布来生成新的样本，根据接受概率来决定是否接受该样本。最后，计算并返回了接受率。

### 回答3：
M-H采样（Metropolis-Hastings Sampling）是一种马尔科夫链蒙特卡洛（MCMC）方法，用于从一个复杂的分布中采样。

在MATLAB中，我们可以使用以下代码实现M-H采样：

1. 定义目标分布函数：

function p = target_dist(x)
    p = 0.3*normpdf(x, -2, 1) + 0.7*normpdf(x, 2, 1);
end

2. 设置初始状态和采样次数：

n_samples = 1000;
x = zeros(1, n_samples);
x(1) = 0; % 初始状态

3. 使用M-H采样算法进行采样：

for i = 2:n_samples
    % 从当前状态生成候选样本
    x_cand = normrnd(x(i-1), 0.5);

    % 计算接受率
    accept_prob = min(1, target_dist(x_cand) / target_dist(x(i-1)));

    % 判断是否接受候选样本
    if rand < accept_prob
        x(i) = x_cand;
    else
        x(i) = x(i-1);
    end
end

4. 绘制采样结果的概率密度函数估计：

histogram(x, 'Normalization', 'pdf');
hold on;
x_range = linspace(min(x), max(x), 100);
y = target_dist(x_range);
plot(x_range, y, 'LineWidth', 2);
legend('采样结果概率密度', '目标分布概率密度');

以上代码将从一个由两个正态分布混合而成的目标分布中进行采样。在M-H采样算法中，通过在当前状态附近生成候选样本，并计算接受率来决定是否接受该候选样本。最终得到的样本序列将逼近目标分布。使用histogram和plot函数可以对采样结果进行可视化。

需要注意的是，M-H采样的效率与候选样本的生成方式和接受率计算方法有关，因此在实际应用中可能需要根据具体问题进行调整和优化。

用matlab模拟马尔可夫链蒙特卡罗拒绝采样

马尔可夫链蒙特卡罗拒绝采样（Markov Chain Monte Carlo Rejection Sampling）是一种基于马尔可夫链的概率采样方法。它是通过建立一个马尔可夫链，并在该链上进行随机游走，最终得到符合所需概率分布的样本。

在Matlab中模拟马尔可夫链蒙特卡罗拒绝采样的步骤如下：
1. 确定目标概率分布P(x)，其中x是要采样的变量。
2. 选择一个合适的马尔可夫链转移矩阵T(x, x')。该矩阵定义了在当前状态x下，转移到下一个状态x'的概率。
3. 选择一个起始状态x0作为马尔可夫链的初始状态。
4. 对于每个迭代步骤n=1,2,3,...，
a. 从当前状态x(n-1)开始，按照转移矩阵T(x, x')的概率选择下一个状态x'(n)。
b. 计算接受率r=min(1, P(x'(n))/P(x(n-1)))。
c. 生成一个均匀分布的随机数u~U(0, 1)。
d. 如果u≤r，则接受状态x'(n)，否则保持当前状态x(n-1)。
5. 重复步骤4直到达到预定的采样次数。

需要注意的是，对于某些复杂的概率分布，上述方法可能无法高效收敛。此时，可以尝试使用其他改进的马尔可夫链蒙特卡罗方法，如Metropolis-Hastings方法或Gibbs采样等。此外，在编写Matlab代码时，应注意使用合适的数据结构和算法，以提高计算效率。