• MCMC and Bayesian Data Analysis(PPT在文件模块)


    How to generate a sample from $p(x)$?

    Let's first see how Matlab samples from a $p(x)$. In Matlab, there are several common probability distributions.

    Try univariate Gaussian distribution

    p= normpdf(linspace(xmin , xmax , num_of_points) , mean, standard_deviation);%PDF
    c= normcdf(linspace(xmin , xmax , num_of_points) , mean, standard_deviation);%CDF
    y=normrnd(mean, standard_deviation,number_of_samples, 1);%Random Number Generating Method
    

    Try PDF:

    x=linespace(-1,1,1000); 
    p=normpdf(x ,0, 1);
    plot(x,p);
    

    Note: linespace returns a vector which is usually accessed like this

    x(1)%the first elem, not x(0)
    x(:)
    x(1,:)

    Try RSM:

    y=normrnd(0, 1,100, 1);%试试采10000个样本
    hist( y , 20 );%20 bars
    

    Try univariate uniform distribution

    p= unifpdf(linspace(xmin , xmax , num_of_points) , a,b);%PDF
    c= unifcdf(linspace(xmin , xmax , num_of_points) , a,b);%CDF
    y=unifrnd(a,b,number_of_samples, 1);%RNG
    

    Try PDF:

    x=linespace(-10,10,1000); 
    p= unifpdf(x ,-5,5);
    plot(x,p);
    

    Try RSM:

    y=unifrnd(-5, 5,100, 1);%试试采10000个样本
    hist( y , 20 );%20 bars

    Matlab provides random number generating functions for some standard $p(x)$, it doesn't provide us sampling functions for a general $p(x)$. Here I show some common sampling methods.

    Inverse Transform Sampling(ITS)

    with descret variables

    This method generates random numbers from any probability distribution given the inverse of its cumulative distribution function. The idea is to sample uniformly distributed random numbers (between 0 and 1) and then transform these values using the inverse cumulative distribution function(InvCDF)(which can be descret or continous). If the InvCDF is descrete, then the ITS method just requires a table lookup, like shown in Table 1.  

    Table 1. Probability of digits observed in human random digit generation experiment

    There is a method called randsample in Matlab that can implement the sampling process using the Table 1. See the code below. 

    %Note: The randsample doesn't defaultly exist in Octave-core package, install statistic package from http://octave.sourceforge.net/statistics/ before using randsample.
    
    % probabilities for each digit
    theta=[0.000; ... % digit 0
    0.100; ... % digit 1
    0.090; ... % digit 2
    0.095; ... % digit 3
    0.200; ... % digit 4
    0.175; ... % digit 5
    0.190; ... % digit 6
    0.050; ... % digit 7
    0.100; ... % digit 8
    0.000];
    
    seed = 1; rand( 'state' , seed );% fix the random number generator
    K = 10000;% let's say we draw K random values
    digitset = 0:9;
    Y = randsample(digitset,K,true,theta);
    figure( 1 ); clf;
    counts = hist( Y , digitset );
    bar( digitset , counts , 'k' );
    xlim([-0.5 9.5]);
    xlabel( 'Digit' );
    ylabel( 'Frequency' );
    title( 'Distribution of simulated draws of human digit generator' );
    pause;

    Instead of using the built-in functions such as randsample or mnrnd, it is helpful to consider how to implement the underlying sampling algorithm using the inverse transform method which is:

    (1) Calculate $F(X)$. 

    (2) Sample u from Uniform(0,1).

    (3) Get a sample $x^{i}$ of $P(X)$, which is $F(u)^{-1}$.

    (4) Repeat (2) and (3) until we get enough samples.

    Note: For discrete distributions, $F(X)^{-1}$ is discrete, the way to get a sample $x^{i}$ is illustrated below where $u=0.8,~x^{i}=6$ .

    with continuous variables

    This can be done with the following procedure:

    (1) Draw U ∼ Uniform(0, 1).

    (2) Set $X=F(U)^{-1}$

    (3) Repeat

    For example, we want to sample random numbers from the exponential distribution where  its CDF is F (x|λ) = 1 − exp(−x/λ) . Then $F(u|gamma)^{-1}=-log(1-u)gamma$. Therefore replace $F(U)^{-1}$ with $F(u|gamma)^{-1}$.

    p=-log(1-unifrnd(0,1,10000,1))*2;
    hist(p,30);

    Reject Sampling

    Applied situation: impossible/difficult to compute CDF of $P(X)$.

    Advantage: unlike MCMC, it doesn't require of any “burn-in” period, i.e., all samples obtained during sampling can immediately be used as samples from the target distribution $p( heta)$.

    Based on the Figure above, the method is:

    (1) Choose a proposal distribution q(θ) that is easy to sample from.

    (2) Find a constant c such that cq(θ) ≥ p(θ) for all θ.

    (3) Draw a proposal θ from q(θ).

    (4) Draw a u from Uniform[0, cq(θ)].

    (5) Reject the proposal if u > p(θ), accept otherwise. Actually, since u is sampled from Uniform[0, cq(θ)], it is equal to state like this " Reject if $uin[p( heta),cq( heta)]$, accept otherwise".

    (6) Repeat steps 3, 4, and 5 until desired number of samples is reached; each accepted sample $ heta$ is a draw from p(θ).

    For example

     

    then the code is 

    k=100000;%draw k samples
    c=2;
    theta_vec=unifrnd(0,1,k,1)%gen a proposal vector from q($	heta$)
    cq_vec=c*unifpdf(theta_vec);%cq(theta) vector
    p_vec=2*theta_vec;%p(theta) vector
    u_vec=[];
    for cq=cq_vec
    u_vec=[u_vec;unifrnd(0,cq)];
    end
    r=theta_vec.*(u_vec<p_vec);
    r(r==0)=[];%remove the “0” elements
    hist(r,20);
    

    MCMC Sampling

    Before getting to know MCMC sampling, we first get to know Monte Carlo Integration and Markov Chain.

     

      For example:

    %Implement the Markov Chain involving x under Beta(200(0.9x^((t-1))+0.05),200(1-0.9x^((t-1)-0.05))
    
    fa=inline('','x')%parameter a for beta
    fb=inline('200*(1-0.9*x-0.05)','x');%parameter b for beta
    no4mc=4;%4 markove chains
    states=unifrnd(0,1,1,no4mc);%initial states
    N=1000;%200 samples drawn from 4 chains
    X=states;
    for i=1:N
    states=betarnd(fa(states),fb(states));
    X=[X;states];
    end;
    plot(X);
    pause;
    

    Metroplis Sampling

  • 相关阅读:
    Hibernate---对象的三种状态
    grunt+bower依赖管理
    grunt 的安装和简单使用
    sqlserver dmv 动态管理视图
    ado.net 数据库连接池
    桥接模式
    .net MVP
    主定理(分治算法)
    图中环的判断
    选举协议paxos 协议 理解
  • 原文地址:https://www.cnblogs.com/chaseblack/p/6001931.html
Copyright © 2020-2023  润新知