0

I am trying to write a matlab function to solve a test integral using the Metropolis Method. My function is listed below.

The integral is from 0 to infinity of x*e(-x^2), divided by the integral from 0 to infinity of e(-x^2)

This function converges to ~0.5 (notably, it does fluctuate about this answer a little) however analytically the solution is ~0.5642 or 1/sqrt(pi).

The code I use to run the function is also below.

What I have done wrong? How do I use metropolis to correct solve this test function?

% Metropolis Method for Integration
% Written by John Furness - Computational Physics, KTH

function [I S1 S2] = metropolis(f,a,b,n,sig)
% This function calculates an integral using Metropolis Mathod
% Only takes input as a function, f, on an interval between a and b,
% where n is the number of points.

%Defining burnin
%burnin = n/20;
burnin = 0;

% Finding maximum point
x = linspace(a,b,1000);
f1 = f(x);
max1 = max(f1);

%Setting Up x-vector and mu
x(1) = rand(1);
mu=0;

% Generating Random Poins for x with Gaussian distribution.
% Proposal Distribution will be the normal distribution
strg = 'exp(-1*((x-mu)/sig).^2)';
norm = inline(strg,'x','mu','sig');

for i = 2:n
% This loop generates a new state from the proposal distribution.
  y = x(i-1) + sig*randn(1);
  % generate a uniform for comparison
  u = rand(1);
  % Alpha is the acceptance probability
  alpha = min([1, (f(y))/((f(x(i-1))))]);

    if u <= alpha
       x(i) = y;
    else
       x(i) = x(i-1);
    end

end


%Discarding Burnin
%x(1:burnin) = [];

%I = ((inside)/length(x))*max1*(b-a);
I = (1/length(f(x)))*((sum(f(x))))/sum(norm(x,mu,sig));

%My investigation variables to see what's happening
%S1 = sum(f(x));
%S2 = sum(norm1(x,mu,sig));
S1 = min(x);
S2 = max(x);

end

Code used to run the above function:

% Code for Running Metropolis Method
% Written by John Furness - Computational Physics

% Clearing Workspace
clear all
close all
clc

% Equation 1
% Changing Parameters for Equation 1
a1 = 0;
b1 = 10;
n1 = 10000;
sig = 2;

N1 = @(x)(x.*exp(-x.^2));
D1 = @(x)(exp(-x.^2));

denom = metropolis(D1,a1,b1,n1,sig);
numer = metropolis(N1,a1,b1,n1,sig);

solI1 = numer/denom
John
  • 1
  • 1
  • Is your primary issue a statistical one or a coding one? – Glen_b Jan 24 '15 at 02:20
  • Did you try increasing the number of points to see if the answer converges to the true solution? You didn't tell us what you have tried other than asking you what's wrong with your code. – rayryeng Jan 24 '15 at 23:15
  • @Glen_b I think it may be my understanding of the mathematics. The code does not produce any errors and I have some experience with MATLAB. I run the function with the following code and it converges to a value of 0.5 (It does fluctuate about this point). The true value for this integral is 1/sqrt(pi). The numerator (numer) and denominator (denom) are both integrals from 0 to infinity. – John Jan 30 '15 at 16:19

0 Answers0