# 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>

[itex]
<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>
[/itex]
<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\left(x\right)={x}^{3}$

Foreign function interface

$\xi =\mathrm{drand48}$

Sample value

${x}_{i}=A+\left(B-A\right){\xi }_{}$

Approximation to integral

$\mathrm{Integral}=\frac{1}{N}\underset{i=1}{\overset{N}{\sum }}f\left({x}_{i}\right)$

Variance in approximation to integral

$\mathrm{Variance}=\frac{1}{N}\underset{i=1}{\overset{N}{\sum }}{f\left({x}_{i}\right)}^{2}-{\left(\frac{1}{N}\underset{i=1}{\overset{N}{\sum }}f\left({x}_{i}\right)\right)}^{2}$

Standard error in approximation to integral

$\mathrm{StandardError}=\sqrt{\frac{\mathrm{Variance}}{\left(N-1\right)}}$

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.