slide1 l.
Skip this Video
Loading SlideShow in 5 Seconds..
Cycle Romand de Statistique, 2009 September 2009 Ovronnaz, Switzerland Random trajectories: some theory and appl PowerPoint Presentation
Download Presentation
Cycle Romand de Statistique, 2009 September 2009 Ovronnaz, Switzerland Random trajectories: some theory and appl

Loading in 2 Seconds...

play fullscreen
1 / 76

Cycle Romand de Statistique, 2009 September 2009 Ovronnaz, Switzerland Random trajectories: some theory and appl - PowerPoint PPT Presentation

  • Uploaded on

Cycle Romand de Statistique, 2009 September 2009 Ovronnaz, Switzerland Random trajectories: some theory and applications Lecture 2 David R. Brillinger University of California, Berkeley 2   1. Lecture 2: Inference methods and some results

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'Cycle Romand de Statistique, 2009 September 2009 Ovronnaz, Switzerland Random trajectories: some theory and appl' - sampson

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

Cycle Romand de Statistique, 2009

September 2009

Ovronnaz, Switzerland

Random trajectories: some theory and applications

Lecture 2

David R. Brillinger

University of California, Berkeley

2 1


Lecture 2: Inference methods and some results

Lecture 1 provided motivating examples

This lecture presents analyses

EDA and CDA (Stefan)


The Chandler wobble.

Chandler inferred the presence of 12 and approx 14 months components in the wobble.

Serious concern to scientists and at the end of the 1800s

Network of stations set up to collect North Star coordinates

Data would provide information on the interior structure of the Earth


Monthly data, t = 1 month.

Work with complex-values, Z(t) = X(t) + iY(t).

Compute the location differences, Z(t), and then the finite FT

dZT() = t=0T-1 exp {-it}[Z(t+1)-Z(t)]


IZZT() = (2T)-1|dZT()|2



Arato, Kolmogorov, Sinai, (1962) set down the SDE

dX = - Xdt -  Ydt +  dB

dY =  Xdt -  Ydt +  dC

Z = X + iY  = B + iC

General stimulus

dZ = -  Zdt +  d =  - i 

Adding measurement noise, the power spectrum is

|i + |-2f()+2|1-exp{-i}|2/2

But what is the source of ? Source of 12 mo, 14 mo


If series stationary, mixing periodograms, Is at  = 2s/T approximately independent exponentials parameter fs

Suggesting estimation criterion (quasi-likelihood)

L = s fs-1 exp{-Is/fs}

and approximate standard errors

Gaussian estimation, Whittle method



Perhaps nonlinearity

Looked for association with earthquakes, atmospheric pressure by filtering at Chandler frequency.

None apparent

Mystery "solved" by modern data and models.

Using 1985 to 1996 data, R. S. Gross (NASA) concluded two thirds of wobble caused by changes in ocean-bottom water pressure, one-third by changes in atmospheric pressure.

NASA interested. One of the biggest sources of uncertainty in navigating interplanetary spacecraft is not knowing Earth's rotation changes.


"Brownian-like" data. Perrin's mastic grain particles

Viscosity, so can't be exactly Brownian

Perrin checking on Einstein and Smoluchowsky

n = 48, t = 30 sec


Potential function.

Quadratic. H(x,y) = γ1x + γ2y + γ11x2 + γ12xy + γ22y2



μ = - grad H = - (γ1 + 2γ11x + γ12y , γ2 + γ12x + 2γ22y )


(r(ti+1)-r(ti))/ (ti+1-ti) =μ(r(ti)) + σZi+1/√(ti+1-ti)


martingale differences

asymptotic normality +, Lai and Wei (1982)



Ornstein-Uhlenbeck like

Potential function for O-U

H = (a - r)'A(a-r)/2 A symmetric  0



H(r) =  log |r| +  |r| + γ1x + γ2y + γ11x2 + γ12xy + γ22y2

r = (x,y)

attraction (goalmouth) plus smooth

|r – a0|, a0 closest point of goalmouth

(r(ti+1)-r(ti))/ (ti+1-ti) =μ(r(ti)) + σZi+1/√(ti+1-ti)

μ = -grad H, stack, WLS


Discussion. Modelled path, not score

Asymmetry, down one side of the field

Ball speed, slow, then quick


Hawaiian Monk seal.

2.2 m, 250 kg, life span 30 yr

Endangered – environmental change, habitat modification, reduction in prey, humans, random fluctuations

~ 1200 remain


Location data.

Satellite-linked time depth recorder + radio transmitter

Argos Data Collection & Location Service

Location estimate + index

(Location class (LC) = 3,2,1,0,A,B,Z)

UTM coordinates – better projection, euclidian geometry, km


Female, 4-5 years old

Released La’au Point 13 April 2004

Study period til 27 July

n = 573 over 87.4 days

(ti ,r(ti), LCi), i=1,…,I unequally spaced in time

well-determined: LC = 3, 2, 1

I = 189

Spatial feature: Molokai



Multi-d generalization of boxplot

Center is multi-d median

Bag contains 50% of observations with greatest depth (based on halfspaces)

Fence separates inliers from outliers – inflates bag by factor of 3

Equivariant under affine transforms




- distance from La’au Point

- foraging?



H(r,t) - two points of attraction, one offshore, one atshore

Potential function

½σ2log |r-a| - δ|r-a|

a(t) changes


Parametric μ= -grad H

Approximate likelihood from

(r(ti+1)-r(ti))/ (ti+1-ti) =μ(r(ti)) + σZi+1/√(ti+1-ti)

Robust/resistant WLS

Estimate σ2 from mean squared error


Discussion and summary.

Time spent foraging in Penguin Bank appeared constrained by a terrestrial atractor (haulout spot – safety?).

Seal spent more time offshore than thought previously


robust/resistant methods basic


Brownian motor. Kinesin

A two-headed motor protein that powers organelle transport along microtubules.

Biophycist's question. "Do motor proteins actually make steps?"

Hunt for the periodic positions at which a motor might dwell

Data via optical instrumentation


Kinesin motor attached to microtubule

Malik, Brillinger and Vale (1994)


Location (X(t),Y(t))

Rotate via svd to get parallel displacement, Z(t)

2 D becomes 1 D


Step function, N(t)?

Z(t) =  + N(t) + E(t)


As stationary increment process fZZ = 2 fNN + fEE

If N(t) renewal

fNN = p(1 - ||2) / (2  |1 - |2), p rate,  characteristic function

Interjump, time j+1 - j constant, v velocity of movement

power spectrum

j(/v - 2j/)

periodic spikes


Prewhitened for greater sensitivity.

Robust line fitted to Z(t)

Periodogram of residuals

Robust line fit to log(periodogram) at low frequencies and subtracted

Averaged results for several microtubules

To assess simulated various gamma distributions


For some l set

Y(t) = t + klklk (t/T) + noise


lk(x) = 2l/2(2lx - k)

Haar scaling wavelet

(x) = 1 0  x < 1

= 0 otherwise

Fit by least squares

Shrink: replace estimate alk by w(|alk|/slk)alk

w(u) = (1-u2)+


Discussion and summary.

"Malik et al (1994) were able to rule out large, regular (that is, strictly periodic) steps for kinesin movement along microtubules (for  > 12 nm) and argued they would not have been able to detect smaller steps (say 8 nm or less) unless the motions were highly regular (quasi periodic), with the step-to-step interval varying by less than 20%"

Since then Brownian motion has been reduced revealing steps



Kernel density estimate

based on the locations r(tj)

Relation with potential function in stationary case

(r) = c exp{-2H(r)/2} H(r) = -½ 2 log (r)/c


Vector field model.

(r(ti+1)-r(ti))/ (ti+1-ti) =μ(r(ti),ti) + σZi+1/√(ti+1-ti)

Robust fit via generalized additive model for smooth μ



The elk appear to be more active at 0600 and 1800 hours, but staying iin a local area at 1200 and 2400 hours.

There are hot spots

The time of day effect does not appear additive.


Elephant seal journey.

Were endangered

Formulas for: great circle, SDE


"Particle" heading towards North Pole

speed 

bivariate Brownian disturbance (U,V).  s.d.

(,): longitude, colatitude

dt =  dUt + (2/2tan t - )dt

dt = (/sin t) dVt

Brownian with drift on a sphere


Parameter estimation.

Discrete approximation

t+1 - t = 2/(2 tan t) -  + t+1

t+1 - t =  / (sin t )t+1

Measurement error

t' = t + t'

t' = t +  / (sin t')t


Results. Likelihood by simulation

No measurement error

parameter estimate s.e.

 .0112 rad .001

 .00805

Measurement error

out  .0126 .0001

in  .0109 .0001

 .000489 .0000004

 .0175 .0011


Discussion and summary.

Mostly location measurement error (done by light levels)


Model appraisal.

Synthetic plots (Neyman and Scott)

Simulate realization of fitted model

Put real and synthetic side by side


Turing test

Compute same quantity for each?


Ringed seal.

Might wish to simulate other such paths

Suppose seal travels in segments with a constant velocity then,

dr = vdt

i.e. the segments are straight.

It may be that the speed |v| is approximately constant for all


Discussion and summary.

Acoustic tracking - attached pinger

Dives forages and surfaces

Finds its way back to breathing hole - need to navigate back to and between holes

straight line segments?

running biweight?


Some formal aspects.

Consider,for example,

H(x,y) = ∑ βkl gk(x)hl(y) with gk , hl differentiable

Hx(x,y) = ∑ βkl g(1)k(x)hl(y), Hy(x,y) = ∑ βkl gk(x)h(1)l(y)

with  = -(Hx,Hy)

(r(ti+1) – r(ti))/(ti+1 – ti) = (r(ti)) + σZi+1

stack the data

linear model


Martingale difference.

E{Xn+1|{X0,...,Xn}) = 0

E{Xn} = 0



E{Sn+1|{S0,...,Sn}) = Sn



Theorem A.1. [Lai and Wei (1982)]. Suppose yi = xiTβ + εi, i=1,2,…with {εi} martingale differences wrt increasing sequence of σ-fields {Fn}. Suppose further that

supn E(||εn||α|Fn-1) < ∞ a.s.

for some α > 2 and that limn→∞ var{εn|Fn-1) = σ2 a.s.for some nonstochastic σ.

Assume that xnFn-1-measurable r.v. and existance of non-random positive definite symmetric L by L matrix Bn for which Bn-1(XnTXn)1/2→I and sup1≤i≤n|| Bn-1xn||→0 in probability.Then

(XnTXn)1/2(b-β) →N(0, σ2I), in distribution as n→∞.


Note that 0 mean independent observations like the σZi+1 of the basic model

(r(ti+1)-r(ti))/ (ti+1-ti) =μ(r(ti),ti) + σZi+1/√(ti+1-ti)

form a martingale difference sequence with respect to the σ-field Fi generated by {r(t1),…,r(ti)}.


CI forφ(r)T β

Theorem A.2. Under the assumptions of Theorem A.1 and lim log λmax(XnTXn)/n→0 almost surely, one has

((φ(r)(XnTXn)-1φ(r)T)-1/2φ(r)T (b-β)/sn → N(0,1)

in distribution as n→∞.

sn = ((n-1)p)-1RSS


Chang and Chin (1995) least squares when

var{εn|Fn-1) = g(zn; )

zn, is observable and Fn-1 measureable

The estimate of  from

min t=1n (êt2 - g(zn; ))2

the êt residuals from the OLS fit

Asymptotic normality results


Summary and discussion.

Array of biological and physical mechanisms control how animals move, particularly on large landscapes.

Models of movement is one useful tool to study ecology of animal behavior and to test ideas concerning foraging strategies, habitat preferences, and dynamics of population densities

Cleaning the data robust/resistant


Constrained trajectories by shrinking

Potential function - effective approach

SDE motivated parameter estimate

New stochastic models result


SDE benefits.

conceptual, extendable, simulation, analytic results, prediction, effective

Potential function benefit

Motivates parametric and nonparametric estimates

difficulties: enforcing boundary


Acknowledgements. Data/background providers and collaborators

Aager, Guckenheimer, Guttorp, Kie, Oster, Preisler, Stewart, Wisdom, Littnan, Roy Mendolssohn, Dave Foley, ?

Lovett, Spector