Wolfram Language Paclet Repository
Community-contributed installable additions to the Wolfram Language
MardiaDryden
|
Author Information
Wilfrid S . Kendall
Statistics, University of Warwick,
Coventry CV4 7 AL, UK.
Purpose
This is a Mathematica notebook that derives the Mardia-Dryden distribution for the shape of three points moving by independent Brownian motions started from different points on the plane.
IMPORTANT
This notebook assumes nothing has been previously defined. Quit Mathematica and reload if this is not the case!
Load package
In[1]:=
Needs["FernandoDuarte`Itovsn3`"]
Specification of problem parameters
We are dealing with pts=3 points moving in n=2 dimensions: the size of the configuration is measured by a Bessel process of dimension dim=n (pts-1) (allows for irrelevant translations). Note that the problem is well-formulated for larger numbers of points and higher dimensions: however we leave it as a (non-trivial!) exercise to generalize this approach accordingly.
In[2]:=
n=2;pts=3;dim=n(pts-1)
Out[2]=
4
Semimartingales and change of time-scale
First we initialize Itovsn3 with basic time semimartingale t and set up a Bessel size process (by modifying a Brownian motion - it would be a minor exercise in Mathematica programming to set up a Bessel procedure, but we choose not to do this here!). Notice the persistent need to declare constants using SetAttributes; a feature of Mathematica is that it chooses to assume implicit dependency unless otherwise stated (it is arguable whether this is (a) annoying or (b) a safety feature or (c) both!).
In[3]:=
 | ItoReset |
 | BrownSingle |
| AddDrift |
(dim-1)dt
2R
Out[3]//TableForm=
Itovsn3 resetting ... |
Itovsn3 initialized |
with time semimartingale t |
and time differential dt |
We need a new timescale dtau to allow for the way that shape changes faster when size is small. We introduce a new innovation Brownian motion W in this new timescale.
In[7]:=
dtau=
[W,dW];
[,dtau];
[dWdR,0];
[dW,0]
[]
dt
2
R
| Introduce |
| AddQuadVar |
2
dW
| AddQuadVar |
| AddDrift |
| ItoStatus |
Out[7]=
dt
2
R
Out[11]=
0
Out[12]=
Summary of current structure of stochastic differentials | ||||||||||||||||
Current second-order structure of semimartingale differentials: | ||||||||||||||||
| ||||||||||||||||
Current first-order structure of semimartingale differentials: | ||||||||||||||||
| ||||||||||||||||
Current initial values: | ||||||||||||||||
|
Definition of shape-space geometry concepts
Now we can introduce the shape-metric distance S from the shape-point representing all equilateral shapes. The argument of the trigonometric functions is 2 S rather than S to allow for the shape geometry, which is that of a sphere of radius 1/2. (Detailed justification of the stochastic differential equation for S follows from work cited in the references.)
In[13]:=
SetAttributes[S0,Constant];
[S,dSdW+Cot[2S]dtau,S0];
[]
| Itosde |
| ItoStatus |
Out[15]=
Summary of current structure of stochastic differentials | |||||||||||||||||||||||||
Current second-order structure of semimartingale differentials: | |||||||||||||||||||||||||
| |||||||||||||||||||||||||
Current first-order structure of semimartingale differentials: | |||||||||||||||||||||||||
| |||||||||||||||||||||||||
Current initial values: | |||||||||||||||||||||||||
|
Surrogate for time
The next step is to produce a surrogate for time, taking into account the Brownian scaling. This will be used as the dispersion parameter. Let t1 be the (real) time for which we allow the original Brownian points to move. The rule kappaSub allows us to rewrite formulae in a more compressed style if required.
In[16]:=
SetAttributes[t1,Constant];kappa=;SetAttributes[k,Constant];kappaSub=First[Solve[kappak,t1]]
2
R
4(t1-t)
Out[19]=
t1+4kt
2
R
4k
Backwards differential equation
Now we can derive the "backwards differential equation" for the Mardia-Dryden density:
In[20]:=
density=f[Cos[2S],kappa];bde=Simplify
[density]/.kappaSub0
 | Drift |
| ItoD |
dt
Out[21]=
1
2
R
(0,1)
f
2
k
(0,2)
f
(1,0)
f
2
Sin[2S]
(2,0)
f
Trial solutions motivated by perturbation theory
Instead of trying to solve this directly, we play around a little with trial solutions. Motivated by geometry, consider a "flat" trial solution and perturbations:
In[22]:=
p0=kappaExp[-kappa(1-Cos[2S])]SetAttributes[a0,Constant];SetAttributes[a1,Constant];SetAttributes[a2,Constant];p1=a0++a2(1-Cos[2S])p0
a1
kappa
Out[22]=
-(1-Cos[2S])
2
R
4(-t+t1)
2
R
4(-t+t1)
Constraints induced by small-time and small-scale
We are going to consider constraints on the constants a0, a1, a2 placed by requirements at small distance from equilateral shape:
Now study the constraint of requiring zero drift at the equilateral triangle shape:
Now consider the conservation of mass:
This must be constant in k so we can differentiate and then set k->0:
Now we can solve all these three constraints:
Our trial answer is therefore
Verification
By computation of drift and pole behaviour
By computation of drift and pole behaviour
Consider the drift (remembering we have to restore some time-dependence by resubstituting kappa for k):
Hence our trial density is correct as long as we check it has the correct behaviour for t near to t1 (its "pole"). This is what we get (a positive limit as t -> t1 for generic S):
By substitution in backwards differential equation
By substitution in backwards differential equation
Note that we can substitute the trial density in the "backwards differential equation" and see if we get the correct answer! (Alternative to computing the Drift.)
Contact Information
Email: w.s.kendall@warwick.ac.uk
URL: http://www.warwick.ac.uk/go/WSK
Acknowledgements
The research reported here was supported by EPSRC grants GR/71677 (Stochastic calculus in AXIOM using modules of stochastic differentials) and GR/L56831 (Perfect simulation in stochastic geometry), and a joint EPSRC/BBSRC research grant (Multi-strain species modelling and control via differential algebra reductions). This Mathematica notebook was constructed on a visit to MSRI Berkeley CA during its 1997-1998 program Stochastic Analysis. Finally, it is a pleasure to express my gratitude to my friends Suzanne Scotchmer and Joseph Farrell for the generous hospitality they showed to me during my visit to MSRI.
References
The original Mardia-Dryden derivation:
General introductions to shape:
◼
Chapter 8 of D. Stoyan, WSK and J. Mecke: Stochastic geometry and its applications (Second Edition), Wiley and Sons, Chichester (1995);
◼
WSK: "Geometry, statistics, and shape", in Geometry in Present Day Science, edited by O.E. Barndorff-Nielsen and E.B. Vedel-Jensen (1998).
Symbolic Ito calculus and Itovsn3 were developed to use stochastic calculus to investigate shape!
◼
WSK: "Symbolic computation and the diffusion of shapes of triads", Advances in Applied Probability, 20:775-797 (1988);
◼
WSK: "The Euclidean diffusion of shape", in "Disorder in Physical Systems", edited by D. Welsh and G. Grimmett, Oxford University Press, Oxford, pp 203-217 (1990);
◼
WSK: "A diffusion model for Bookstein triangle shape", Advances in Applied Probability, 30 number 2 (1998);
This particular application is also implemented in AXIOM:
◼
WSK: "Symbolic Ito calculus: an ongoing story", Research Report 327, Department of Statistics, University of Warwick. Submitted for publication. (1998).