Seite 1 von 2
Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Do Okt 25, 2018 9:10 pm
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))
}
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Do Okt 25, 2018 9:33 pm
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
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Fr Okt 26, 2018 6:55 am
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
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Fr Okt 26, 2018 8:37 am
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.
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Fr Okt 26, 2018 11:45 am
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
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Fr Okt 26, 2018 1:25 pm
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
kopieren, in die man die gewünschten Änderungen einbringt.
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
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Mo Okt 29, 2018 9:27 am
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
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Mo Okt 29, 2018 9:39 am
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
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Mo Okt 29, 2018 10:30 am
von maf7wz
Es funktioniert! Recht herzlichen Dank dafür, auch für die Zeit! Das hat mir sehr weitergeholfen.
Viele Grüße
Re: Area under Curve (auc) Befehl - Fehler bzgl. number of subdivisions
Verfasst: Mo Okt 29, 2018 11:56 am
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