Importance Sampling[1] 

                                 

Introduce a positive definite sampling function g(x) and let

          

With a positive definite g, G is single valued function that starts at 0 for x = a and ends at 1 for x =b.  Thus t is a number between 0 and 1.  Notice that x has move the the upper limit of the integral.   

With this definition for the function x(t)

                             

So that

          

This is discussed at length for analytically integral functions for which an immediate solution for x(t) is possible in Laurent.doc.  the Laurent transform is for the 0 to range while an arc tangent transform is used for the - to range. 

Semi-infinite range

                Write the integral

 

Let

 

 

The region (a,b) in becomes

  

The value of the integral given in becomes

            

The value of x(t) is needed.  Equation becomes

 

               

This solves as

          

Numerically

For t=1, x=ln(1-1)/α=∞, but computers do not handle ln(1-1) very well.    In order to leave a few digits for the last term the ending point for the integral in tmax should be (1-10-13).  This means that the integration in t extends only to 13×2.3/α. 

Figure 1  r2(t)exp(-2r(t)) versus t
100, 200, 400 pts are used toevaluatethe integral

  AN     2.500000151518428E-01  100 pts

 AN2     2.500000011811078E-01 200 PTS

 AN4     2.500000000742313E-01 400 PTS

 ANR     2.500000000566618E-01 Richardson’s extrapolation

ANRF   2.500000005080060E-01  ± 4.443478404E-10

Fit of two points to AN, AN2, AN4  [..\Fittery\nlfit-r\StdDev\3ptLinFit.docx]  Answer is ¼.  The code is in SampleInt.zip

Solving for x(t)with an integrable function

The value of x(t) can be found by solving the equation

          

Newton’s method, ..\optimization\solving\Newton.doc .htm, involves expanding the as a function of x about  x0, setting the result equal to zero and solving for the next value of x

  

So that

                                 

Note that G’ is g so that the sequence

                  

can be iterated to find x(t).  Note that this could take a lot of computer time if every value of G(x0;a) requires a revaluation of the integral in the numerator of .  Normally, there would be a Lagrange interpolation of a single set of points, but this can introduce errors.  An exact method involves making g(x) explicitly the straight line connecting a set of values g(xi).

 

Line connecting the points modification.

                A very simple g(x) is the line connecting a group of points

           

Let y= x-xi-1 so that

         

For y = xi-xi-1 the positive first term cancels half of the negative part of the second term leading to the seemingly linear

               

The values in are the exact values of G(xi) for the g(x) that is the line connecting the points g(xi).  The values between these points are given exactly by .

Figure 2  Integral of line connecting points (black).  Linear interpolation of this same line.  The g values are g(1)=3, g(2)=6,g(3)=5.  The integral is below the linear interpolation for a positive g’ and below it for a negative g’.  Code is in infosamp.zip.

Solving for x(t)in line connecting the points modification

The value of G(xi) form an ascending series of values starting at 0 and ending at 1. 

      

The subroutine LOCATE(tG(NMAX),G,NMAX,J) ..\interpolation\Locate.doc returns a value J such that G(J)£t£G(J+1).  In the region X(J) £ x £ X(J+1)

         

 Equation becomes

               

Define

       

So that

                        

This is a quadratic equation with a general solution given by ( ..\solving\Quadratic.doc). 

           

C is greater than or equal to zero, since locate returned a J such that tG(a,b) > G(J).  B is always less than 0, while the sign of A is unknown.

 Equation is numerically unsuitable since it will involve large cancellations.   Following ( ..\solving\Quadratic.doc) rewrite as

          

The + sign yields y > 0.  So

       

The largest value of C is G(J+1)-G(j) =( g(J+1)+g(J))(X(J+1)-X(J))/2 so that the largest value of 4AC/B2 is

This means that the most negative value of the argument of the square root in is

        

Thus the value of y is never imaginary.

Summary

                Find an arrangement of g(J)=g(xJ).  Use equation to find G(J).  Then for values of t between 0 and 1, use locate(tG(NMAX),G,NMAX.J) to find the relevant J.  Use to define the terms in which yields y(t).  Finally

       


[1] J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods, Methune & Co. Ltd, London, John Wiley & Sons Inc, New York.  pp. 57 -59