这是一个基于蒙特卡洛算法的matlaab程序，有错误实在是修改不出来，希望明白的大神给些修改意见。

运行顺序

/gf_mmc_driver.m

/mcmc1.m

/mcmc2.m

            

gf_mmc_driver.m

s

%  "This is the driver script to be run from the command line. Also %  requires functions mcmc1 and mcmc2. This code is referred to by an %  upcoming article by T. A. Driscoll and K. Maki in the Education %  section of SIAM Review, to appear 2006-2007."%%  MMC iteration applied to the growth factors of random matrices.  %  Runs 'numberchain' instances of MMC simultaneously using the Distributed%  Computing Toolbox.% N = 8;                              % size of nxn matrixbinedge = -.5:1:(2^(N-1)+.5);       % edges of target binsNbin = length(binedge) - 1;         % number of binsM = 5e5;                            % number of realizations per iterationNiter = 70;                         % number of iterationsburnin = 5e5;                       % burn-in period for MCMCnumberchain = 12;                   % number of independent MMC runs % ***** Initialization ***** randn('state',sum(100*clock)), rand('state', sum(100*clock));  Accept = zeros(Niter,numberchain);  % acceptance rate per iteration per runsigma = (1/N)*ones(numberchain,1);  % sigma parameter per runH = zeros(Nbin,Niter,numberchain);  % hits/bin per iteration per runP = ones(Nbin,numberchain)/Nbin;    % computed histogram per rung = zeros(Nbin-1,Niter,numberchain);% parameter needed for Berg's iterationPnew = zeros(Nbin,1);               % updated inverse weights GFinitial = zeros(numberchain,1);   % computed initial growth fac for MCX = zeros(N,N,numberchain);         % initial matrices for MC runschaingrowth = 10*(1:numberchain-2); % desired initial growth fac for MC % Establishing communication for distributed computingjm = findResource('jobmanager','Name','woprdce');   % Calculating one initial matrix for Markov ChainX(:,:,1) = randn(N,N);  % ***** Start MMC iteration ***** for i=1:Niter        % Calculating remaining initial matrices for Markov Chain    X(:,:,2)=randn(N);        % Implementing Higham and Higham's Theorem    for z=3:numberchain        while( GFiter(z)>chaingrowth(z-2)+.5 || GFiter(z) < chaingrowth(z-2)-.5 )            Mstart = toeplitz([1 -.99999*(-1+chaingrowth(z-2)^(1/(N-1)))*ones(1,N-1)], [1 zeros(1,N-1)]);            T = triu( (2*rand(N-1)-1)/N );            T = [T;zeros(1,N-1)];            Tlc = (chaingrowth(z-2).^((0:N-1)/(N-1)));            D = diag( sign(randn(N,1)) );            X(:,:,z) = D*Mstart*[T Tlc'];            [L,U]=lu(X(:,:,z));            GFiter(z) = (max(abs(U(:))))/(max(max(abs(X(:,:,z)))));        end    end     % ***** Running MCMC for burn-in *****        % Telling the Distributed Computing Toolbox to complete one job with     % 'numberchain' tasks.  Each task is comprised of running a MCMC     % for the burnin time (mcmc1.m) with a different initial matrix.     j = createJob(jm);    set(j,'Timeout',1200);    set(j,'FileDependencies',{ 'mcmc1.m'});    for z = 1:numberchain        createTask( j, @mcmc1, 2, {X(:,:,z), burnin, sigma(z), P(:,z), binedge});    end    submit(j);    waitForState(j,'finished');    data = getAllOutputArguments(j);        % Received computed data and recording results     for z = 1:numberchain        Accept(i,z) = data{z,1};        X(:,:,z) = data{z,2};    end    clear j     % Adjusting the parameter sigma based on burnin results    for z = 1:numberchain        if( Accept(i,z)/burnin < .1 )            sigma(z) = max( sigma(z)/2, eps);        elseif( Accept(i,z)/burnin > .4 )            sigma(z) = 2*sigma(z);        else            sigma(z) = sigma(z);        end    end     % ***** Running MCMC for M *****        % Again telling the Distributed Computing Toolbox to complete one job     % composed of 'numberchain' tasks.  Each task continues running the     % independent MCMC iterations (script mcmc2.m) from above.  Now    % statistics are calculated.        j = createJob(jm);    set(j,'Timeout',1200);    set(j,'FileDependencies',{'mcmc2.m'});    for z=1:numberchain        createTask(j,@mcmc2, 3, {X(:,:,z), M, sigma(z), P(:,z), binedge});    end    submit(j);    waitForState(j,'finished');    data = getAllOutputArguments(j);        % Received computed data and recording results    for z = 1:numberchain        Accept(i,z) = data{z,1};        H(:,i,z) = data{z,2};        X(:,:,z) = data{z,3};    end    clear j     % Tuning the parameter sigma again    for z = 1:numberchain        if( Accept(i,z)/(M) < .1 )            sigma(z) = max( sigma(z)/2, eps);        elseif( Accept(i,z)/(M) > .4 )            sigma(z) = 2*sigma(z);        else            sigma(z) = sigma(z);        end    end        % *****Berg update*****    for z =1:numberchain        Pnew = P(:,z);           for k = 1:(Nbin-1)            if (H(k+1,i,z)*H(k,i,z) > 0)                g(k,i,z) = H(k+1,i,z)*H(k,i,z)/(H(k+1,i,z)+H(k,i,z));                ghat(k) = g(k,i,z)/(sum(g(k,1:i,z)));                Pnew(k+1) = Pnew(k)*(P(k+1,z)/P(k,z))*((H(k+1,i,z)/H(k,i,z))^ghat(k));            else                Pnew(k+1) = Pnew(k)*(P(k+1,z)/P(k,z));            end        end        P(:,z) = Pnew/sum(Pnew);    end    end

mcmc1( X, burnin, sigma, P,binedge

% Run MCMC for the burnin period.

function [ Accept, X ] = mcmc1( X,burnin, sigma, P, binedge )

   acceptval = rand(burnin,1);     %all acceptance tests

   N = length(X);                  %size of matrices

   Accept = 0;                     %number of matrices accept during run

  

   % Calculating growth factor of initial matrix

   [L,U] = lu(X);                    

   growth = max( abs(U(:)) )/max( abs(X(:)) );

   [junk, bin] = histc(growth,binedge);

   

   % *****MCMC starts*****

   for m = 1:burnin         

       Xpro = X + sigma*trnd(8,N,N);      % proposal matrix

      

       % calculating growth factor of proposal matrix

       [L,Upro] = lu(Xpro);      

       growthXpro = max( abs(Upro(:)))/max( abs(Xpro(:)) );

       [junk, newbin] = histc(growthXpro,binedge);

      

       % calculation pi(X) and pi(X_pro)

       pi_x = exp(-.5*sum(sum(reshape( X.^2, [N^2 1]))));   

       pi_y = exp(-.5*sum(sum(reshape( Xpro.^2, [N^2 1]))));

      

       % Accept proposal?

       alpha = min(1, ((pi_y*P(bin))/(pi_x*P(newbin))));

       if acceptval(m) < alpha

            X = Xpro;

            bin = newbin;

            growth = growthXpro;

            Accept = Accept + 1;

       end

   end

end

mcmc2( X,M, sigma, P, binedge )

% Run of MCMC --- statistics recorded

function [ Accept, Hits, X ] = mcmc2( X,M, sigma, P, binedge )

   acceptval = rand(M,1);                 % all acceptance test

   N = length(X);                         % size of matrix

   Accept = 0;                            % number of acceptances

   Hits = zeros( length(binedge)-1,1 );   % hits/bin for current iteration

   % Calculating initial growth factor

   [L,U] = lu(X);

   growth = max( abs(U(:)) )/max( abs(X(:)) );

   [junk, bin] = histc(growth,binedge);

   % *****MCMC Starts*****

   for m=1:M        

      

       Xpro = X + sigma*trnd(8,N,N);      % proposal matrix

      

       % Calculating growth factor of proposal

       [L,Upro] = lu(Xpro);

       growthXpro = max( abs(Upro(:)))/max( abs(Xpro(:)) );

       [junk, newbin] = histc(growthXpro,binedge);

      

       % Calculating pi(X) and pi(X_pro)

       pi_x = exp(-.5*sum(sum(reshape( X.^2, [N^2 1]))));

       pi_y = exp(-.5*sum(sum(reshape( Xpro.^2, [N^2 1]))));

      

       % Accept proposal?

       alpha = min(1, ((pi_y*P(bin))/(pi_x*P(newbin))));

       if acceptval(m) < alpha

            X = Xpro;

            bin = newbin;

            growth = growthXpro;

            Accept = Accept + 1;

       end

      

       % Record statistics

       Hits(bin) = Hits(bin)+1;

   end

end