metropolis是一種取樣方法，一般用於獲取某些擁有某些比較複雜的概率分佈的樣本。
1.取樣最基本的是隨機數的生成，一般是生成具有均勻分佈的隨機數，比如C++裡面的rand函式，可以直接採用。
2.取樣複雜一點是轉化發，用於生成普通常見的分佈，如高斯分佈，它的一般是通過對均勻分佈或者現有分佈的轉化形成，如高斯分佈就可以通過如下方法生成：
令:U,V為均勻分佈的隨即變數；
那麼：Gauss=sqrt(-2lnU)cos(2*pi*V)。
(詳細參照Box-Muller)
像這種高斯分佈matelab可以直接用randN生成。


3. 最後一種方法是類MCMC方法，其實MCMC是一種架構，這裡我們說一種最簡單的Metropolis方法，它的核心思想是構建馬爾剋夫鏈，該鏈上每一個節點都成為一個樣本，利用馬爾剋夫鏈的轉移和接受概率，以一定的方式對樣本進行生成，這個方法簡單有效，具體可以參考《LDA數學八卦》的，裡面都已經進行了詳細的介紹。本文給一個小demo，是利用C++寫的。例外程式碼的matlab版本也有，但不是我寫的，一併貼出來，matlab的畫圖功能可以視覺化的驗證取樣有效性，（好像matlab有個小bug），


待生成樣本的概率分佈為：pow(t,-3.5)*exp(-1/(2*t)) （csdn有沒有想latex那種編輯方式啊？？），這裡在metropolis中假設馬爾剋夫鏈是symmetric，就是p(x|y)=p(y|x)，整個程式碼如下：

#include <iostream>
#include <time.h>
#include <stdlib.h>
#include <cmath>
using namespace std;


#define debug(A) std::cout<<#A<<": "<<A<<std::endl;


class MCMCSimulation{

	public:
		int N;
		int nn;
		float* sample;

		MCMCSimulation(){

			this->N=1000;
			this->nn=0.1*N;
			this->sample=new float[N];
			this->sample[0]=0.3;

		}

		void run(){

			for (int i=1;i<this->N;i++){

				float p_sample=this->PDF_proposal(3);	

				float alpha=PDF(p_sample)/PDF(this->sample[i-1]);

				if(PDF_proposal(1) < this->min(alpha,1)){
				
					this->sample[i]=p_sample;

				}else{

					this->sample[i]=this->sample[i-1];
				
				} 
				//debug(p_sample);
				//debug(PDF_proposal(1));
				//debug(sample[i]);
				cout<<sample[i]<<endl;
			}

		}
	private:

		//the probability desity function
		float PDF(float t){

			return pow(t,-3.5)*exp(-1/(2*t));	

		} 

		//generate a rand number from 0-MAX, with a specified precision
		float PDF_proposal(int MAX){

			//srand
			return (float)(MAX)*(float)(rand())/(float)(RAND_MAX);

		}

		float min(float x, float y){
		
			if(x<y){
				
				return x; 

			}else{
			
				return y;
			}
		
		}

};

int main(){

	MCMCSimulation* sim =new MCMCSimulation();
	sim->run();
	delete sim;
	return 1;
}
 

matlab版本

%% Simple Metropolis algorithm example
%{
---------------------------------------------------------------------------
*Created by:                  Date:            Comment:
 Felipe Uribe-Castillo        July 2011        Final project algorithm  
*Mail: 
 [email protected]
*University:
 National university of Colombia at Manizales. Civil Engineering Dept.
---------------------------------------------------------------------------
The Metropolis algorithm:
First MCMC to obtain samples from some complex probability distribution
and to integrate very complex functions by random sampling. 
It considers only symmetric proposals (proposal pdf). 

Original contribution: 
-"The Monte Carlo method"
  Metropolis & Ulam (1949).
-"Equations of State Calculations by Fast Computing Machines"
  Metropolis, Rosenbluth, Rosenbluth, Teller & Teller (1953).
---------------------------------------------------------------------------
Based on:
1."Markov chain Monte Carlo and Gibbs sampling"
   B.Walsh ----- Lecture notes for EEB 581, version april 2004.
   http://web.mit.edu/~wingated/www/introductions/mcmc-gibbs-intro.pdf
%}
clear; close all; home; format long g;

%% Initial data
% Distributions 
nu           = 5;                                       % Target PDF parameter
p            = @(t) (t.^(-nu/2-1)).*(exp(-1./(2*t)));   % Target "PDF", Ref. 1 - Ex. 2
proposal_PDF = @(mu) unifrnd(0,3);                      % Proposal PDF
% Parameters
N        = 1000;              % Number of samples (iterations)
nn       = 0.1*(N);           % Number of samples for examine the AutoCorr
theta    = zeros(1,N);        % Samples drawn form the target "PDF"
theta(1) = 0.3;               % Initial state of the chain

%% Procedure
for i = 1:N
   theta_ast = proposal_PDF([]);          % Sampling from the proposal PDF
   alpha     = p(theta_ast)/p(theta(i));  % Ratio of the density at theta_ast and theta points
   if rand <= min(alpha,1)
      % Accept the sample with prob = min(alpha,1)
      theta(i+1) = theta_ast;
   else
      % Reject the sample with prob = 1 - min(alpha,1)
      theta(i+1) = theta(i);
   end
end

% Autocorrelation (AC)
pp = theta(1:nn);   pp2 = theta(end-nn:end);   % First ans Last nn samples
[r lags]   = xcorr(pp-mean(pp), 'coeff');
[r2 lags2] = xcorr(pp2-mean(pp2), 'coeff');

%% Plots
% Autocorrelation
figure;
subplot(2,1,1);   stem(lags,r);
title('Autocorrelation', 'FontSize', 14);
ylabel('AC (first samples)', 'FontSize', 12);
subplot(2,1,2);   stem(lags2,r2);
ylabel('AC (last samples)', 'FontSize', 12);

% Samples and distributions
xx = eps:0.01:10;   % x-axis (Graphs)
figure;
% Histogram and target dist
subplot(2,1,1); 
[n1 x1] = hist(theta, ceil(sqrt(N)));
bar(x1,n1/(N*(x1(2)-x1(1))));   colormap summer;   hold on;   % Normalized histogram
plot(xx, p(xx)/trapz(xx,p(xx)), 'r-', 'LineWidth', 3);        % Normalized function
grid on;   xlim([0 max(theta)]);   
title('Distribution of samples', 'FontSize', 14);
ylabel('Probability density function', 'FontSize', 12);
% Samples
subplot(2,1,2);     
plot(theta, 1:N+1, 'b-');   xlim([0 max(theta)]);   ylim([0 N]);   grid on;
xlabel('Location', 'FontSize', 12);
ylabel('Iterations, N', 'FontSize', 12); 

%%End
 

下面是樣本的分佈圖