\documentstyle[a4,12pt,amstex,harvard,psfig,norsk]{article}

\setlength{\parindent}{0cm}
\setlength{\parskip}{\baselineskip}

\title{SIF 5064 Romlig Statistikk, \\{\O}ving 1: Gaussiske felt og Kriging}
\author{}
\date{}

\begin{document}

\maketitle
{\bf Innledning}
Denne {\o}vingen handler om Gaussiske felt, kriging,
kovariansfunksjoner (variogrammer) m.m.  \O vingen er laget slik at 
det den benytter spatial biblioteket i Splus
Kapittel 15 i Venables \& Ripley: Modern Applied Statistics with Splus, og notatet
``S-plus ved IMF, Del~3: Funksjoner, bibliotek og internett'', 
samt web linker fra hjemme siden kan v\ae re nyttige. 
''spatial''  laster du inn  inn i S-plus ved {\aa} skrive 
\texttt{library(spatial)}.

\subsection*{Oppgave 1: 1D Gaussisk felt}

\textbf{(a)} Lag plott som viser hvordan eksponensiell og 2.ordens
eksponensiell  kovariansfunksjon ser ut for 
forskjellige valg av ''sill'', og  ''range''. 
Gj\o r parameter valg som gj\o r det mulig \aa\ sammenligne de
forskjellige korrelasjons funksjonene. Gi en tolkning av parameterene.
 
\textbf{Hint:} Benytt funksjonene  \texttt{gaucov} og \texttt{expcov} 
fra biblioteket ``spatial'' i splus.

Hva er den teoretiske  sammenhengen mellom kovariansfunksjon,
og variogram  ? N\aa r vil du bruke den ene parametriseringen fremfor
den andre.  

\textbf{(b)} Trekk et normalisert 1D gaussisk felt diskretisert p\aa\
et grid av st\o rrelse 50,
ved \aa\ bruke multigaussisk fordeling. Bruk forskjellige covarians
funksjoner. Diskuter innvirkningen av parameter valgene p\aa\ simuleringen.

\textbf{Hint:} Under ``Examples''
p{\aa} hjelpeteksten for \texttt{rnorm} finner du en funksjon som du
kan benytte.  

\textbf{(c)} Velg \'{e}n av realisasjonene med 1. ordens eksponensiell
kovarians fra oppgave (b).  Bruk verdiene i node 10,25 og 30 som
observasjoner. Finn utrykk for betinget fordeling for nodene som ikke
er observert gitt de tre observasjonene. Plot betinget forventning og
betinget varians (eller standard avvik ) i dette eksempelet. Plot
ogs\aa\ realisasjonen du valgte. Diskuter resultatene.
 

\textbf{(d)} Lag 25 realisasjoner fra det Gaussiske feltet betinget p{\aa}
observasjonene i (c).  Plott de betingede realisasjonene.  Sammenlign
med resultatene i (b) og (c) og kommenter.

\textbf{(e)} Beregn gjennomsnitt og empirisk varians for de  25
realisasjoner punktvis dvs., for hver av grid nodene.
Plot og kommenter resultatet. 


\textbf{(f)} Definer en terskel verdi p\aa\ $1.0$. En \o nsker \aa\ vite
hvor stort arealet over terskel verdien og under den utvalgte
realisasjonen er. For \aa\ predikere dette arealet, betinget p\aa\ de tre
observasjonen kan en tenke seg to prediktorer.
\begin{enumerate}
\item  arealet mellom betinget forventning fra (c)
    og terskel verdien.
\item  gjennomsnittet av arealet over terskelverdien, og under  
 gjentatte simuleringer fra den betingede fordelingen i (d).
\end{enumerate}
Benytt begge prediktorene for \aa\ gi en prediksjon. 
Hvilken prediktor gir st\o rst verdi ?
Kan du gi en teoretisk begrunnelse for dette ?
Hvordan vil du presentere usikkerheten i prediktoren
for strategi 2 ?
Beregn det sanne arealet over terskel verdien og under  din utvalgte  
realisasjon fra (b). Sammenlign dette estimatet med prediksjonen og
usikkerheten i strategi 2.


\textbf{Tips:} Splus syntaxen \texttt{sum(y[y>0])} kan forenkle
programmet.


\textbf{(g)} Repeter oppgave (c)-(f) med 2.ordens eksponensiell
kovarians funksjon. Kan noen av forskjellene du finner tilskrives
forskjell i kovarians funksjoner.
 
\textbf{Merk:} Dette eksempelet er en skjematisk presentasjon av "oil
in place" problemet, vi estimerer berg volumet under en geologisk
horisont og over olje vann kontakten.


\begin{figure}[hbtp]
    \begin{center}
        \leavevmode
        \psfig{figure=oljevolum.ps,width=9.5cm,angle=270}
        \caption{Figuren fremstiller situasjonen over.}
        \label{Oljevolum}
    \end{center}
\end{figure}

\subsection*{Oppgave 2: Trendflater, kriging og plotting i S-plus}

Funksjoner og data som brukes i denne oppgaven finnes i biblioteket
\texttt{spatial}.  Data settet
\texttt{topo} inneholder 52 m{\aa}linger av h{\o}yden over havet (i fot)
innenfor et $310 \times 310$ fot kvadrat, omregnet til 50 fots enheter.

\textbf{(a)} Vi vil f{\o}rst se p{\aa} trendflater, polynomiske
regresjonsflater av formen
\[
    f((x,y)) = \sum_{r+s \leq p} a_{rs} x^r y^s
\] 
hvor parameteren $p$ er polynomets orden.  
Tilpass trendflater med orden  p=  2,3,4,6. Bruk splus funksjonen 
\texttt{surf.ls}. Plot de forskjellige flatene med plot funksjonene:
;\texttt{contour}, \texttt{image} and
\texttt{persp} inkluder data punkter i disse plottene.
Funksjonen \texttt{trmat} kan v\ae re nyttig.
Diskuter fordeler og ulemper med de forskjellige plottings rutinene.
Kommenter resultatene.
Hvor mange parametere m{\aa} estimeres for en 6. ordens flate?  
Hvilke svakheter eller ulemper medf{\o}rer dette?

\textbf{(b)} Ved universell kriging tilpasses en flate av
formen
\[
    Z(x) = f(x)^T \beta + \epsilon(x) 
\]
der $x$ er 2D romlig referanse, $f(x)$ er gitte funksjoner, $\beta$ er
ukjente koeffisienter og  $\epsilon(x)  $  er et Gaussisk felt med forventning 0
og romlig kovariansfunksjon:
\[
    \mbox{cov} \{Z(x), Z(y) \} = \mbox{cov} \{ \epsilon(x), \epsilon(y) \}
\]

Pr{\o}v universell kriging, benytt funksjonene \texttt{surf.gls} med
varians 1, 1.ordens eksponensiell kovarains med parameter $d=0.7$.
La $f(x)$ v\ae re polynomer som i (a), og benytt forskjellige verdier
for $p$,dvs polynomer av forskjellig orden. Benytt spesielt $p=0$.
Plot og kommenter resultatene.
Funksjonen \texttt{prmat} kan her v\ae re nyttig. 
Pr\o v ogs\aa\  forskjellige verdier av $d$.

Endre til 2. ordens eksponensiell kovarians funksjon og gj\o r
beregningene paa nytt. Plot og kommenter.

\textbf{(c)} Plot den estimerte trendflaten fra kriging modellene i
(b), benytt \texttt{trmat}.

Hva er forskjellen mellom dette estimatet og det vi fikk i a, dvs. hva er
forskjellen mellom funksjonene \texttt{surf.gls} og \texttt{surf.ls}.

Hva er forskjellen p{\aa} trendflatene og krigingflatene?

\textbf{(d)} Plot prediksjonsstandardavvik for krigings flaten
fra (b), Kommenter resultatene. Benytt splus funksjonen \texttt{semat}.

\subsection*{Oppgave 3: Simulering og  estimering}

Simulering av Gaussiske felt i dimensjon 2 eller h{\o}yere blir fort
ganske regnetungt.  I denne oppgaven skal vi bruke et ferdig program
for simulering av 2D Gaussiske felt.  Programmet simulerer et Gaussisk
felt med forventning null, varians 1 og eksponensiell kovariansfunksjon
$c(r) = \exp( -3 r/d)$ der $r$ er avstanden mellom punktene.
Input til
programmet er gridst{\o}rrelse og parameteren $d$ i
kovariansfunksjonen.  \textbf{NB:} Merk 3-tallet!  I S-plus er
eksponensiell kovariansfunksjon definert som $c(r) = \exp(-r/d)$.
 Programmet
startes ved {\aa} skrive \texttt{/local/kurs/romstat/simGauss
    \$RANDOM}.  
Det kryptiske
\texttt{\$RANDOM} gir fr{\o} (seed) til randomgeneratoren.  Hvis det
utelates, brukes en standard verdi.  
Programmet skriver ut resultatet til en fil med en matrise av verdiene
til det Gaussiske feltet i tilh{\o}rende gridnoder.

\textbf{(a)} Simuler et Gaussisk felt med
gridst{\o}rrelse: $50 \times 50$ og parameter $d=20$. 
Pr{\o}v  et par andre parametervalg ogs{\aa}. Plot realisasjonene i
splus.
Benytt \texttt{ sim <- read.table( ''filnavn'' ) } for \aa\ lese inn datafilen.
Ta vare p{\aa} realisasjonene, de skal brukes
videre i oppgaven.  

\textbf{(b)} ``Observer'' det Gaussiske feltet i 36  punkter 
tilfeldig fordelt utover i gridet.  
Bruk funksjonen \texttt{correlogram} til \aa\ estimere 
den empiriske korrelasjons funksjonen.
Plott estimatet og sammenlign med den sanne korrelasjons funksjonen.

\textbf{(c)} Estimer korrelasjons funksjonen ved \aa\ bruke  
9, 25, 49, 64 og 100  punkter  tilfeldig fordelt utover i gridet.
Omtrent hvor mange observasjoner trenger du for {\aa} et
tilfredsstillende resultat?  

\textbf{(d)} La observasjonspunktene v{\ae}re systematisk
plassert p{\aa} et grid.  Hvordan fungerer dette i forhold til {\aa}
la observasjonenspunktene v{\ae}re tilfeldig fordelt?  (Bruk
tilsvarende antall punkter som du gjorde i forrige delsp{\o}rsm{\aa}l.)

\textbf{(e)} Hvilke problemer kan du se i forbindelse med estimering av
kovariansfunksjon dersom  funksjonsformen  er ukjent?

\end{document}
