Airy-Funktionen kommen in der Physik an verschiedenen Stellen vor. Die Airy-Funktion wurde von George Biddell Airy als er Regenbögen untersuchte und dabei feststellte, dass beim genaueren hinschauen sich die Farben wiederholen. In R ist die Airy-Funktion über das Paket Bessel verfügbar. Wie eine Messung eine Airy-Funktion aussehen kann, ist im folgenden zu sehen:

library(Bessel)
library(ggplot2)

x <- -(1:10000)*1e-2
t <- 0.1 #length
y0 <- 0.5 #y offset
α <- 0.3# proportianality factor
y <- y0 + α * ( AiryA(x*t) + runif(length(x))/100 )

df <- data.frame(x=x,y=y)
plt <- ggplot(df,aes(x=x,y=y)) + geom_line()
plt

Wenn man jetzt irgendwo auf eine Airy-Funktion stößt, oder man vermutet dass man auf eine Airy-Funktion gestoßen ist, dann möchte man diese auch anfitten. Eine Funktion die hierfür in Frage kommt, ist die Funktion nls in R. Dabei benötigt man aber Startwerte, die idealerweise schon ziemlich gut mit dem gesuchten Wert übereinstimmen müssen. Im Beispiel ist t=0.1, wenn der Startwert nur um 0.016 abweicht, dann funktioniert der Fit mit nls nicht mehr.

startValues <- list(y0 = 1, α = 0.5, t = 0.116)

fit <- nls(y ~ y0 + α * ( AiryA(x*t) ),
           data = df,
           start = startValues)

pars <- fit$m$getPars()
pars
##        y0         α         t 
## 0.5178934 0.1121109 0.1221287
df$yFit <- pars[1] + pars[2] * ( AiryA(x*pars[3]) )
plt <- ggplot(df,aes(x=x,y=y)) + geom_line()
plt + geom_line(aes(x=x,y=yFit),colour="red")

Um an sinnvolle Startwerte zu kommen, braucht man eine Nutzbare Heuristik. Für y0 lässt sich einfach der Mittelwert aller Funktionswerte nutzen, den Proportionalitätsfaktor α kann man über den maximalen und minimalen Funktionswert abschätzen.

y0Guess <- mean(df$y)
αGuess <- 3/5 * (max(df$y) + min(df$y))/2

Um ein passendes t zu finden, gibt es jetzt die Möglichkeit, verschiedene t ‘auszuprobieren’. Um zu beurteilen, welches t am besten ist, lässt sich der Kreuzkorrelation zwischen den berechneten Werten und den zu fittenden Werten nutzen:

tSteps <- ((1:200)*5e-3)
#tSteps
airy <- function(t) {
  return( cor( 
    y0Guess + αGuess * ( AiryA(df$x*(t)) ),
    df$y)
    )
}


library(parallel) 
foo <- unlist(mclapply(tSteps,airy)) #use two cores

#max(foo)
#which.max(foo)
tSteps[which.max(foo)]
## [1] 0.1

Damit hat man nun geeignete Startwerte gefunden.

startValues <- list(y0 = y0Guess, 
                     α = αGuess, 
                     t = tSteps[which.max(foo)])

startValues
## $y0
## [1] 0.524474
## 
## $α
## [1] 0.3114035
## 
## $t
## [1] 0.1

Mit den Startwerten lässt sich nls füttern, und liefert anschließend gute Ergebnisse:

fit <- nls(y ~ y0 + α * ( AiryA(x*t) ),
           data = df,
           start = startValues)

pars <- fit$m$getPars()
pars
##        y0         α         t 
## 0.5015069 0.3000118 0.1000002
df$yFit <- pars[1] + pars[2] * ( AiryA(x*pars[3]) )
plt <- ggplot(df,aes(x=x,y=y)) + geom_line()
plt + geom_line(aes(x=x,y=yFit),colour="red")

Damit klappt der Fit dann auch jedes mal, wenn man am Anfang die Parameter der Airy-Funktion zufällig generiert.

t <- runif(1) #length
y0 <- runif(1) #offset
α <- runif(1)# proportianality factor