Programming in Mathematical Notation

This example uses random sampling to approximate an integral.

Several new issues arise in this example, include handling Greek symbols, combing multiple sums into a single loop, and how to interface with a random number generator.

Random Sampling Integration Example

The input file (test10.xml)

<?xml version="1.0"?>

<program>
<comment>
   <center>
      <h1>Simple Monte Carlo Integration Example</h1>
   </center>
</comment>

<math>
  <comment>
     <p>Lower limit</p> 
  </comment>
  <declare type="fn" constant="true">
     <ci>A</ci>
     <lambda> 
       <cn>0</cn>
     </lambda> 
  </declare>

  <comment>
     <p>Upper limit</p> 
  </comment>
  <declare type="fn" constant="true">
     <ci>B</ci>
     <lambda> 
       <cn>1</cn>
     </lambda> 
  </declare>

  <comment>
     <p>Number of samples</p> 
  </comment>
  <declare type="fn" constant="true" output="int">
     <ci>N</ci>
     <lambda> 
       <cn>100</cn>
     </lambda> 
  </declare>

  <comment>
     <p>Function to integrate</p> 
  </comment>
  <declare type="fn" time_series="true">
     <ci>f</ci>
     <lambda> 
       <bvar><ci>x</ci></bvar>
       <apply><power/>
         <ci>x</ci>
         <cn>3</cn>
       </apply>
     </lambda> 
  </declare>


  <comment>
    <p>Foreign function interface</p>
  </comment>
  <declare type="ffi" lang="c" index="true">
    <ci>&#x3be;</ci>
    <ffi_name>drand48</ffi_name>
  </declare>


  <comment>
     <p>Sample value</p> 
  </comment>
  <declare type="fn" input="int" output="double" index="true">
     <ci>x</ci>
     <lambda> 
       <bvar><ci>i</ci></bvar>
       <apply><plus/>
         <ci>A</ci>
         <apply><times/>
           <apply><minus/>
             <ci>B</ci>
             <ci>A</ci>
           </apply>
           <apply> <ci>&#x3be;</ci>
           </apply>
         </apply>
       </apply>
     </lambda> 
  </declare>

  <comment>
     <p>Approximation to integral</p> 
  </comment>
  <declare type="fn" input="int">
    <ci>Integral</ci>
    <lambda> 
      <apply><times/>
        <apply><divide/>
          <cn>1</cn>
          <ci>N</ci>
        </apply>
        <apply><sum/>
          <bvar><ci>i</ci></bvar>
          <lowlimit><cn>1</cn></lowlimit>
          <uplimit><ci>N</ci></uplimit>
          <apply><ci>f</ci>
            <apply><ci>x</ci>
              <ci>i</ci>
            </apply>
          </apply>
        </apply>
      </apply>
    </lambda> 
  </declare>

  <comment>
     <p>Variance in approximation to integral</p> 
  </comment>
  <declare type="fn" input="int">
    <ci>Variance</ci>
    <lambda>
      <apply><minus/>
        <apply><times/>
          <apply><divide/>
            <cn>1</cn>
            <ci>N</ci>
          </apply>
          <apply><sum/>
            <bvar><ci>i</ci></bvar>
            <lowlimit><cn>1</cn></lowlimit>
            <uplimit><ci>N</ci></uplimit>
            <apply><power/>
              <apply><ci>f</ci>
                <apply><ci>x</ci>
                  <ci>i</ci>
                </apply>
              </apply>
              <cn>2</cn>
            </apply>
          </apply>
        </apply>
        <apply><power/>
          <apply><times/>
            <apply><divide/>
              <cn>1</cn>
              <ci>N</ci>
            </apply>
            <apply><sum/>
              <bvar><ci>i</ci></bvar>
              <lowlimit><cn>1</cn></lowlimit>
              <uplimit><ci>N</ci></uplimit>
              <apply><ci>f</ci>
                <apply><ci>x</ci>
                  <ci>i</ci>
                </apply>
              </apply>
            </apply>
          </apply>
          <cn>2</cn>
        </apply>
      </apply>
    </lambda>
  </declare>

  <comment>
     <p>Standard error in approximation to integral</p> 
  </comment>
  <declare type="fn" input="int">
    <ci>StandardError</ci>
    <lambda>
      <apply><root/>
        <apply><divide/>
          <ci>Variance</ci>
          <apply><minus/>
            <ci>N</ci>
            <cn>1</cn>
          </apply>
        </apply>
      </apply>
    </lambda>
  </declare>
</math>
<results>
 <ci>StandardError</ci>
 <ci>Integral</ci>
 <ci>Variance</ci>
</results>
</program>

Greek characters are represented by the construct &#x(unicode number). This is the general way of representing arbitrary Unicode characters in XML. (In the future, it could be possible to translate names, or use a DTD that defined entities such as &xi;, which is how the Greek characters are represented in the xhtml output ) For the text and C++ output, the names of the characters are spelled out.

To handle calling the random number generator, a C function, we define a declaration type of "ffi" (foreign function interface). The "lang" attribute is in the example, but is not used yet. The "ffi_name" tag describes the name of the function in the foreign language.

In this example, we want to compute the average, variance, and standard error of the function evaluated at a series of values. The problem is these quantities are defined in separate sums, but we only want to evaluate the summand function once. The solution here is to create a "time_series" attribute on the function being evaluated. This indicates that the function will return a different value every time it is evaluated, and consequently, it should only be called once per iteration.

Also, we add a "results" section at the end. If any functions in the results section depends on a time series function, the sums inside those functions must be combined.

The text output is

A = 0
B = 1
N = 100
f(x) = x^3
ffi:  xi  =  drand48
x(i) = A+(B-A)*xi()()
Integral = (1/N)*sum_{i=1}^N (f(x(i)))
Variance = (1/N)*sum_{i=1}^N (f(x(i))()^2)-((1/N)*sum_{i=1}^N (f(x(i))))^2
StandardError = sqrt (Variance/(N-1))

The output in MathML presentation form (test10.xhtml)

Simple Monte Carlo Integration Example

Lower limit

A = 0

Upper limit

B = 1

Number of samples

N = 100

Function to integrate

f x = x 3

Foreign function interface

ξ = drand48

Sample value

x i = A + ( B - A ) ξ

Approximation to integral

Integral = 1 N i = 1 N f ( x i )

Variance in approximation to integral

Variance = 1 N i = 1 N f ( x i ) 2 - ( 1 N i = 1 N f ( x i ) ) 2

Standard error in approximation to integral

StandardError = Variance ( N - 1 )

The C++ output (test10.cpp)

double A = 0.0;

double B = 1.0;

int N = 100;

double f(double x){
    return pow(x,3);
}

double x(int i){
    return A + (B - A) * drand48();
}



void StandardError_and_Integral_and_Variance(double &StandardError,double &Integral,double &Variance){
    double sum = 0.0;
    double sum2 = 0.0;
    int i;
    for (i = 1;i <= N;i++) {
        double tmp = f(x(i));
        sum += tmp;
        sum2 += tmp * tmp;
    }
    double tmp2 = (1.0 / N) * sum;
    Variance = (1.0 / N) * sum2 - tmp2 * tmp2;
    StandardError = sqrt(Variance / (N - 1.0));
    Integral = (1.0 / N) * sum;
}

To make use of the C++ output, use the following code (main10.cpp)

#include <stdio.h>
#include <math.h>
#include <stdlib.h>



#include "test10.cpp"

int main()
{
  double in,var,err;
  StandardError_and_Integral_and_Variance(err,in,var);
  printf("integral = %g Variance = %g Error = %g\n",in,var,err);
  return 0;
}

Previous page: An example using the trapezoidal rule to approximate an integral.

Next page: Discussion and the conversion program

Written by Mark Dewing on July 20, 2005. Last updated September 21, 2005.