Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Wie erweitere ich R um eigene Funktionen oder Pakete? Welches Paket ist passend für meine Fragestellung?

Moderatoren: EDi, jogo

maf7wz
Beiträge: 7
Registriert: Do Okt 25, 2018 8:11 pm

Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von maf7wz »

Guten Abend,

ich habe eine Beispieldatei, bei der die etwas weiter unten stehenden 3 Spalten von Bedeutung sind: 1. Zeit, 3. eine vorgegebene Sinuskurve (Sollkurve), und 2. eine Kurve, die die Sinuskurve nachzuahmen versucht (Istkurve). Ich möchte ab einem Startpunkt, den ich festgelegt habe, die Fläche unter der Kurve aus der 2. Spalte berechnen, konkret für eine Wellenlänge. Dazu nutze ich den auc-Befehl. Ich habe das auch für die "1. Wellenlänge" der Kurve aus Spalte 3 gemacht, da klappt alles. Bei Spalte 2 jedoch erscheint der Fehler "Error in integrate(myfunction, lower = from, upper = to) : maximum number of subdivisions reached".

Ich habe bereits gelesen, dass man die number of subdivisions erhöhen kann, allerdings bekomme ich es nicht hin, das auf dieses Beispiel zu übertragen. Über eine Hilfe wäre ich sehr dankbar. Den Beispieldatensatz_1 habe ich hoch geladen. (die Routine ist sicherlich kürzer bzw. einfacher möglich, bitte nicht wundern).
Vielen Dank im Voraus!


Dies ist der aktuelle Code:

Code: Alles auswählen

library(MESS)
require(MESS)


a<-list.files("C:/Users/schlg/Desktop/AUC_Test", 
              full.names=TRUE)
Name<-list.files("C:/Users/schlg/Desktop/AUC_Test",    
                 full.names=FALSE)


for (i in 1:length(a)){
  options(warn=1)
  wo=scan(a[i], sep="\n", what=" ")
  ab=grep("X_Value", wo)
  länge=length(wo)
  segment=länge-ab
  b<-read.table(a[i], skip=ab, nrows=segment, dec=",", fill=TRUE)
  
  
  Name1=strsplit(Name[i], "_")[[1]][2]
  
  xWerte=b[[1]]
  yIst=b[[2]]
  ySoll=b[[3]]

  
  trigger1soll=(which(ySoll>12.5)[1])-1 #determie start of sine
  
  length=length(xWerte)-trigger1soll
  lengthh=as.numeric(length)
  
  Sollkurve=ySoll[trigger1soll:(trigger1soll+lengthh)]
  Istkurve=yIst[trigger1soll:(trigger1soll+lengthh)]
  xWerte_neu=xWerte[trigger1soll:(trigger1soll+lengthh)]
  
  
  
  #horizontal graph to detect intersections with sine
  x=rep(12.5, each=(length(Sollkurve)))

  
  #determining intersection Points
  above=Sollkurve>x
  intersect.points<-which(diff(above)!=0)
  
  Sollkurve.slopes<-Sollkurve[intersect.points+1]-
    Sollkurve[intersect.points]
  x.slopes<-x[intersect.points+1]-x[intersect.points]
  
  x.points<-intersect.points + ((x[intersect.points] - 
                                   Sollkurve[intersect.points]) / (Sollkurve.slopes-x.slopes))
  y.points<-Sollkurve[intersect.points] + 
    (Sollkurve.slopes*(x.points-intersect.points))
  
  
  #Integral Istkurve 1. Wellenlänge
  
  
  xx=x.points[c(TRUE, FALSE)]#Remove every other element to 
  #obtain full wave length (red circles)
  yy=y.points[c(TRUE, FALSE)]#only necessary for plot
  
  xx1=xx[1]
  xx2=xx[2]
  
  xWerte_1=xWerte_neu[xx1:xx2]
  yIst_1=Istkurve[xx1:xx2]
  
  #head(xWerte_1)#zeigt erste Daten an
  #tail(xWerte_1)#zeigt letzte Daten an
  
  Int_Ist_Total1=auc(xWerte_1,
                      yIst_1, 
                      type = 'spline')

  
  
  #Plot Gesamtdaten  
  plot(yIst[trigger1soll:(trigger1soll+lengthh)], 
       type="l", col="red", main=
         "Kurvenverläufe", ylab='', 
       xlab="Time", xlim=c(0, 9000))
  lines(Sollkurve, type="l", col="blue")
  lines(x, type="l")
  points(x.points, y.points, col="green")
  points(xx, yy, col="red")
  legend(6000, 2, c("Istkurve", "Sollkurve"),
         fill=c("red", "blue"), cex=0.85)
  
  abline(v=xx)
  
  axis(side=1, at=seq(0, 8000, by = 1000))
  
}
Dateianhänge
Beispieldatensatz_1.txt
(248.69 KiB) 67-mal heruntergeladen
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von jogo »

Hallo maf7wz,

willkommen im Forum!
Schau mal in die Hilfe von integrate(). Dort steht subdivisions = 100L. Die Funktion auc() kenne ich nicht; die stammt wohl aus dem Paket MESS. Wenn die Programmierer der Funktion auc() clever waren, dann hat auch die Funktion auc() einen Parameter "..." und kann vielleicht den Parameter subdivisions= an integrate() weiterreichen.
Bitte probier doch, den zusätzlichen Parameter beim Aufruf von auc() zu setzen, z.B. subdivisions = 1000L

Gruß, Jörg
maf7wz
Beiträge: 7
Registriert: Do Okt 25, 2018 8:11 pm

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von maf7wz »

Hallo Jörg,

danke und auch für die Antwort! Bisher hat dies leider noch nicht geklappt, da integrate(...) ein sog. "f"-Element verlangt, welches den Typ "function" erfüllt. Die Kurve bei mir ist ja sozusagen stochastisch und zufällig, daher habe ich bisher Probleme, diese Bedingung ("function") zu erfüllen. Ich probiere es noch weiter, vielen Dank jedenfalls schon einmal!

Viele Grüße
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von jogo »

Hallo maf,

Du hattest geschrieben:
Bei Spalte 2 jedoch erscheint der Fehler "Error in integrate(myfunction, lower = from, upper = to) : maximum number of subdivisions reached".
Da steht ganz klar, dass die Funktion integrate() aufgerufen wird. Wenn Du selber den Aufruf nicht in Deinem Code stehen hast, dann wird integrate() von einer der Funktionen aufgerufen, die Du in Deinem Code verwendest.
Hast Du eine Vermutung, welche Funktion das sein könnte?

Gruß, Jörg
p.s.:
Ich scheue mich, mich direkt mit Deinem Code auseinanderzusetzen, weil ich ihn für schlecht lesbar halte. Er macht den Eindruck, als ob er nicht dafür geschrieben wurde, dass andere ihn begutachten sollten.
maf7wz
Beiträge: 7
Registriert: Do Okt 25, 2018 8:11 pm

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von maf7wz »

Hallo Jörg,

leider nicht. Ich werde entweder eine anderen Befehl für die Flächenberechnung probieren, oder probieren die Spalte 1 und 2 als Funktion f darzustellen. Dann könnte ich über "integrate(f, subdivisions = ...) evtl. zur Lösung kommen.

In der Tat, die Struktur des Codes ist nicht optimal und schwer nachzuvollziehen. Ich dachte, ich spiele von vornherein mit offenen Karten :-)

Danke
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von jogo »

Hallo maf,
maf7wz hat geschrieben: Fr Okt 26, 2018 11:45 am leider nicht.
aber Du weißt, dass der Wert eines bestimmten Integrals nichts anderes ist als die Fläche unter einer Kurve.
Mir scheint (laut Hilfetext der Funktion auc()), dass die Funktion auc() die beiden Vektoren nimmt, approx() aufruft, um eine Funktion zu konstruieren, und dann integrate() verwendet um das bestimmte Integral auszurechnen.
Genaue Information bekommt man, wenn man sich den Quellcode von auc() anschaut: Notfalls kann man sich diesen Code auch für eine selbstgestrickte Funktion

Code: Alles auswählen

myauc <- function(...) ...
kopieren, in die man die gewünschten Änderungen einbringt. :idea:
Im Quelltext von auc() sieht man auch, dass die Zusatzargumente nur an approx() übergeben werden (der Autor der Funktion war doch nicht so clever :? ).
Ich werde entweder eine anderen Befehl für die Flächenberechnung probieren, oder probieren die Spalte 1 und 2 als Funktion f darzustellen. Dann könnte ich über "integrate(f, subdivisions = ...) evtl. zur Lösung kommen.
ich würde den obigen Weg gehen, da auc() komplett als R-Code verfügbar ist.
In der Tat, die Struktur des Codes ist nicht optimal und schwer nachzuvollziehen. Ich dachte, ich spiele von vornherein mit offenen Karten :-)
Ist diesbezüglich eine Überarbeitung geplant?

Gruß, Jörg
maf7wz
Beiträge: 7
Registriert: Do Okt 25, 2018 8:11 pm

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von maf7wz »

Hallo Jörg,

sorry für die späte Antwort und Danke für deine Zeit und Hilfe. Hier ist zunächst einmal eine deutlich reduzierter Code, abgestimmt auf die erste Wellenlänge der Beispieldatei, siehe Excel-Mappe ("AUC"). Hier der Code passend zu AUC.xlsx:

Code: Alles auswählen

library(MESS)
require(MESS)
library(readxl)

AUC <- read_excel("C:/Users/schlg/Desktop/AUC.xlsx", 
                  col_names = TRUE)

Integral_Istwerte=auc(AUC$"Zeit", AUC$"Ist", type = 'spline')
Integral_Sollwerte=auc(AUC$"Zeit", AUC$"Soll", type = 'spline')
Bezüglich der Struktur von "auc" inkl. der erstellten Näherungskurve muss ich mich noch intensiver beschäftigen. Leider sind meine R-Kenntnisse auch beschränkt.

Danke und viele Grüße
Dateianhänge
AUC.xlsx
(42.96 KiB) 89-mal heruntergeladen
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von jogo »

Code: Alles auswählen

library(MESS)
library(readxl)

setwd("~/Desktop/R.Zeug/maf")
AUC <- read_excel("AUC.xlsx", col_names = TRUE)

myauc <- function (x, y, from = min(x), to = max(x), subdivisions=100L,
                   type = c("linear", "spline"), absolutearea = FALSE, ...) 
{
  type <- match.arg(type)
  if (length(x) != length(y)) 
    stop("x and y must have the same length")
  if (length(unique(x)) < 2) 
    return(NA)
  if (type == "linear") {
    if (absolutearea) 
      y <- y - min(y)
    values <- approx(x, y, xout = sort(unique(c(from, to, 
                                                x[x > from & x < to]))), ...)
    res <- 0.5 * sum(diff(values$x) * (values$y[-1] + values$y[-length(values$y)]))
    if (absolutearea) 
      res <- res - min(y) * (max(x) - min(x))
  }
  else {
    if (absolutearea) 
      myfunction <- function(z) {
        abs(splinefun(x, y, method = "natural")(z))
      }
    else myfunction <- splinefun(x, y, method = "natural")
    res <- integrate(myfunction, lower = from, upper = to, subdivisions=subdivisions)$value
  }
  res
}

Integral_Istwerte=myauc(AUC$"Zeit", AUC$"Ist", type = 'spline') ## Standardwert subdivisions=100L, Fehler
Integral_Istwerte=myauc(AUC$"Zeit", AUC$"Ist", type = 'spline', subdivisions=100L)   ## identisch
Integral_Istwerte=myauc(AUC$"Zeit", AUC$"Ist", type = 'spline', subdivisions=1000L)  ## anderer Wert, jetzt geht es
Integral_Sollwerte=myauc(AUC$"Zeit", AUC$"Soll", type = 'spline')
Das eigentliche Arbeitstier für Deinen Zweck sind diese beiden Zeilen:

Code: Alles auswählen

myfunction <- splinefun(x, y, method = "natural")
res <- integrate(myfunction, lower = from, upper = to, subdivisions=subdivisions)$value
(also splinefun() konstuiert die Funktion)

Gruß, Jörg
maf7wz
Beiträge: 7
Registriert: Do Okt 25, 2018 8:11 pm

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von maf7wz »

Es funktioniert! Recht herzlichen Dank dafür, auch für die Zeit! Das hat mir sehr weitergeholfen.
Viele Grüße
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions

Beitrag von jogo »

Der Zugriff auf die Spalten kann noch aufgehübscht werden:

Code: Alles auswählen

Integral_Istwerte  <- myauc(AUC$Zeit, AUC$Ist, type = 'spline', subdivisions=1000L)  ## anderer Wert, jetzt geht es
Integral_Sollwerte <- myauc(AUC$Zeit, AUC$Soll, type = 'spline')
Gruß, Jörg
Antworten