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

1 / 76

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

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.

## Cycle Romand de Statistique, 2009 September 2009 Ovronnaz, Switzerland Random trajectories: some theory and appl

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)]

Periodogram

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

Model.

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

Discussion.

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

real-valued

drift.

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

stack

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

WLS

martingale differences

asymptotic normality +, Lai and Wei (1982)

Discussion.

Ornstein-Uhlenbeck like

Potential function for O-U

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

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.

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

Bagplot.

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

Robust/resistant

Journeys?

- distance from La’au Point

- foraging?

Modelling.

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

Potential function

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

a(t) changes

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

EDA

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

Model

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

with

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

Starkey.

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 μ

Discussion.

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

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.

 .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

Assessment

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

CLT

Martingale.

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

CLT

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→∞.

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

simulation

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