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.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>ξ</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>ξ</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 ξ, 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)
Lower limit
Upper limit
Number of samples
Function to integrate
Foreign function interface
Sample value
Approximation to integral
Variance in approximation to integral
Standard error in approximation to integral
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.