Nichtlineares Optimierungsproblem

einfachster ansatz, ist für mich der einsatz von straffunktionen

das heißt du schreibst dir ne funktion, in der du die parameter übergibst, sie werden anhand der nebenbedingugen auf gültigkeit geprüft und bei verletzung gibst du "Inf" zurück und wenn alles gut ausgeht, den wert der quadratischen gleichung

das kannst schließlich durch "optim" in R jagen

€: wollte es gerade nachbauen für optim ... wo zur hölle ist M??? und geh mal etwas näher auf die Notation ein. also was zeile spalte und tralala bedeuten soll

also es gilt die Beziehung D(tilde) = D+M ....das M siehst du in der Summenformel nicht, weil es im D(tilte) schon steckt. Da D festgelegt wird ist M sozusagen die Variable.

Hier nochmal als Bild :

also das n steht für die n-te Potenz, die andern beiden Symbole für Zeile und Spalte.
Somit wird über alle Elemente der Matrizen in den jeweiligen Potenzen summiert. (genau genommen wird die quadratische Abweichung immer summiert).

Ich frag mich wie du alle Nebenbedingungen erfassen willst.

auf die gefahr hin, dass ich zu blöd bin, das zu verstehen

hier bedeutet das n die potenz des elementes in der Zeile µ und spalte \nu (wieso wird hier keine latexunterstützung angeboten ... echt ma ps.), wenn ich richtig verstehe

aber was soll es dann hier bedeuten?

und die nebenbedingungen würde ich bei der holzhammermethode über straffunktionen einführen. d.h. ist die nebenbedingung verletzt gibst du in R einfach den wert unendlich zurück. sind alle NB erfüllt gibst du den funktionswert zurück

€: M scheint ja ebenfalls eine Matrix zu sein ... diese in der Funktion optim als Parameter zu nutzen sollte aber nicht das Problem sein (habe derzeit keine muße mich in nloptr genauer rein zu wurschteln)

also es bedeutet, dass n die n-te Potenz der Matrix ist, und dann davon die Zeile müh und Spalte nüh, also erst Matrix potenzieren und dann die Elemente betrachten.

Zu deiner zweiten frage: das n bei dem D welches unten steht, bedeutet genau das selbe, also n-te Potenz der Matrix. Warum die bei dem D(tilde) und beim D das n an unterschiedlichen Stellen indexiert haben ist mir auch rätselhaft aber es soll wirklich das selbe bedeuten.

Mir ist nicht klar, wie du die Strafffunktionen mit Optim kombinieren willst.

alles klar ... ich glaube ich habs nun gerafft
und ich glaube nen kleinen fehler in deiner erklärung gefunden zu haben:

\tilde D_{1,mü,nü} ^n ist das element deiner matrix hoch n und
D_{n,mü,nü} is die matrix hoch n und dann das element müh und nüh

somit ergeben sich unterschiedliche werte, sonst wäre das optimale M glaube immer 0

 
# cmpfun ist nicht notwendig, aber bytecode bietet sich bei ständigen aufrufen einfach mal an  potmat <- cmpfun( function(A,n){
  if(NROW(A)!=NCOL(A)) stop("Matrix not quadratic")
  B <- diag(1,NROW(A))
  for(i in 1:n) B <- B%*%A  
  B
})

minfun <- function(para,data,N=2){
    # data >= 0
  if(any(data<0)) return(Inf)
  
  # zeilensummen = 1
  rows <- sapply(1:NROW(data),function(i)sum(data[i,]))  
  if(any(!(rows==1))) return(Inf)
  # \tilde D erzeugen
  D2 <- data + matrix(para,ncol = NCOL(data))
  
  # Bedingung an letzte Spalte in D2
  if(any(!sapply(1:(NROW(D2)-1),function(j) D2[j,NCOL(D2)]<=D2[j+1,NCOL(D2)]))) 
    return(Inf)
  # Ordnungs in D2 von der Diagonale aus
  crit <- sapply(1:NROW(D2),function(i){
    erg <- c()
    if(i > 1) erg <- c(erg,sapply(1:(i-1),function(j) D2[i,j]<=D2[i,j+1]))
    if(i < NROW(D2))erg <- c(erg,sapply(i:(NCOL(D2)-1),function(j) D2[i,j]>=D2[i,j+1]))
    any(!erg)
  })
  if(any(crit)) return(Inf)
  
  # Funktionswert selbst
  return(
    sum(sapply(1:NROW(D2),function(i)sapply(1:NCOL(D2),function(j) sum(
      sapply(2:N,function(n)
        # Angabe der Summenfunktion
        ( (D2[i,j]^n - potmat(data,n)[i,j])/potmat(data,n)[i,j] )^2
      )))))  
  )
  
}

ich habe jetzt alle NB als strafffunktion eingeführt, d.h. sind diese nicht erfüllt, gebe ich als wert unendlich zurück

habe mir dementsprechend ein bsp. gebaut
um optim sinnvoll auszuführen, sollte die zielfunktion für den startparameter lösbar sein, sonst mosert er rum

Test <- matrix(c(0.7,0.1,0.3,0.9),2,2)
optim(par = rep(0,length(Test)),fn = minfun,data=Test,N=4)
 

der R-Code is sicherlich nicht der beste, shcönste und performanteste, aber er geht
sollte ich das problem getroffen haben, dann bitte sehr
wenn nicht, dann sorry für die wall of text

sollte noch was unklar sein, und ich helfen könnte, sag bescheid

und wenn alles scheiße is, dann is das halt so

€: im paket "expm" ist der operator %^% für potenzieren von Matrizen vorhanden, so dass du diesen nutzen kannst und meine Hilfsfunktion weglassen kannst

Hallo erstmal vielen Dank für die Mühe, die du dir gibst!

\tilde D_{1,mü,nü} ^n ist nicht das Element der Matrix hoch n, sondern erst die Matrix hoch n (siehe Bild)

also tilde D_n ist gleich (tilde D_1)^n

das M wird nicht 0, schließlich haben wir zig Nebenbedingungen die gelten müssen. Wenn D also diesen Nebenbedingungen nicht gerecht wird, dann ist M nicht 0.

 
psoptprob <- function(data,N=2){
  # Kompilieren von R-Funktionen
  require(compiler)
  # Matrixpotenz
  require(expm)
  
  # Eingangsdaten prüfen
  # data >= 0
  if(any(data<0))stop("Eingangsdaten nicht korrekt")
  
  # zeilensummen = 1
  rows <- sapply(1:NROW(data),function(i)sum(data[i,]))  
  if(any(rows!=1))stop("Eingangsdaten nicht korrekt") 
  
  # Strafsystem
  rek <- 1
  punish <- 10
  
  # Konsturktion eines InitialParameters der einfach nur die NB erfüllt
  initpara <- matrix(0,NROW(data),NCOL(data))
  
  # Erfüllung der Diagonalbedingung
  for(i in 1:NROW(initpara)){
    if(i>1){
      for(j in ((i-1):1)){
        initpara[i,j] <- min(data[i,j+1]+initpara[i,j+1]-data[i,j],0)
      }
    }
    if(i < NROW(initpara)){      
      for(j in ((i+1):NROW(initpara))){
        initpara[i,j] <- min(data[i,j-1]+initpara[i,j-1]-data[i,j],0)
      }      
    }
  }
  
  # Spaltenbedingung erfüllen
  s <- NCOL(data)
  for(i in 1:(NROW(initpara)-1)){
    diff <- data[i+1,s] + initpara[i+1,s] - data[i,s] - initpara[i,s]
    initpara[i,s] <- ifelse(diff>=0,initpara[i,s],initpara[i,s]+diff)
  };rm(s,i,j)
  
  # Anzeige des Initialparameters
  print(initpara)
  # Anzeige von von Dtilde
  print(data+initpara)
  
  
  # Zielfunktion als Bytecode
  quaddist <- cmpfun(function(Dtilde){
    sum(sapply(1:NROW(Dtilde),function(i)sapply(1:NCOL(Dtilde),function(j) sum(
      sapply(2:N,function(n)
        # Angabe der Summenfunktion
        (( (Dtilde%^%n)[i,j] - (data%^%n)[i,j])/( (data%^%n)[i,j] ))^2
      )))))
  })
  
  # Funktion für optim als Bytecode
  minfun <- cmpfun(function(para){
        
    # Starfe verschärfen
    rek <<-rek+1
    
    # D2 = Dtilde
    D2 <- data + matrix(c(para),ncol = NCOL(data))
    
    # NB das Parameters
    # Bedingung an letzte Spalte in D2
    if(any(sapply(1:(NROW(D2)-1),function(j) D2[j,NCOL(D2)]>D2[j+1,NCOL(D2)]))){
      # print("Spaltenbeding. verletzt")
      return(quaddist(D2)+punish^rek)
    }
    # Ordnungs in D2 von der Diagonale aus
    crit <- sapply(1:NROW(D2),function(i){
      erg <- c()
      if(i > 1) erg <- c(erg,sapply(1:(i-1),function(j) D2[i,j]<=D2[i,j+1]))
      if(i < NROW(D2))erg <- c(erg,sapply(i:(NCOL(D2)-1),function(j) D2[i,j]>=D2[i,j+1]))
      any(!erg)
    })
    if(any(crit)){
      # print("Diagbed. verletzt")
      return(quaddist(D2)+punish^rek)
    }
    
    # Funktionswert selbst
    # print("Kein Fehler")
    return(quaddist(D2))
    
  })
  
  
  
  print(erg<-optim(par = c(initpara),fn = minfun))
  erg <- matrix(erg$par,NROW(data))
  
  print(data + erg)
  
  erg

  
}

Test <- matrix(c(0.3,0.1,0.7,0.9),2,2)
psoptprob(data = Test,N = 2)

Test <- matrix(runif(9),3)
Test <- lapply(1:NROW(Test),function(i)Test[i,]/sum(Test[i,]))
Test <- do.call(rbind,Test)
psoptprob(data = Test,N = 4)

aaaaaaaaaaaaaaaalso

ich habe nun versucht alles zu berücksichtigen, was geht
meine straffunktion werfen nun nicht sofort unendlich zurück, sondern werden mit jedem optimierungsschritt "schlimmer"

desweiteren erstelle ich mir vor der optimierungsausführung einen Initialparameter, so dass alle Nebenbedingungen erfüllt werden

die funktion ist dadurch ziemlich voll geworden
Im Prinzip folgendes:
1. Eingabe der Matrix D
2. Prüfen ob die Matrix D die Nebenbedingungen erfüllt, da D immer gleich bleibt. sollte das nicht der fall sein, wird der Vorgang direkt abgebrochen
3. Konstruktion eines gültigen Initialparameters anhand der NB für Dtilde
4. Definition der Minimierungsfunktion, die den Einhalt der Nebenbedingung prüft
5. Ausführung von Optim
6. Rückgabe des Zielparameters

Die Funktion spuckt bei der Ausführung dieverse Matrizen aus
zuerst den initialparameter
dann Dtilde (zum guggn ob die NB auch wirklich erfüllt werden)
ergebnis von optim
dann das ergebnis von D tilde (ob NB immer noch erfüllt sind)
und gibt schließlich den zielparameter zurück

bei unklarheiten und weiteren fehlern einfach melden ... bzw, es sieht so aus, als ob die Anwendung von optim hier doch sinnloser ist, als ich dachte

Wieso wird der Vorgang abgebrochn wenn D die NB nicht erfüllt? Die Nebenbedingungen müssen nur für D_tilde gelten. Der Sinn dieser Optimierung ist quasi dass die Migrationsmatrix D, die aus statistischen Daten erhoben wird durch die optimierung intuitiv geglättet wird. Wenn D die Nebenbedingunen erfüllt, dann ist M doch immer 0.

Für D muss doch Zeilensumme = 1 gelten und alle Daten größer Null
diese Anforderungen werden nie an Dtilde gestellt

da nur Dtilde (D+M) mit Hilfe von M sich ändert, müssen die NB für D nur einmal geprüft werden. sollten diese nicht erfüllt sein, wird das mit bel. M auch nie der fall sein. Oder ich hab wieder irgendwas übersehen ?

€: hab noch ne kleine Korrektur am Code vorgenommen. Da ich 2 verschiedene Initialparameter aufstellen kann und von da aus die Optimale Lösung suche.
Ich habe festgestellt, dass ich somit beide Initialisierungen durchlaufen lasse und dann den Zielparameter wähle, welcher den niedrigeren Funktionswert auswirft

psoptprob <- function(data,N=2){
  # Kompilieren von R-Funktionen
  require(compiler)
  # Matrixpotenz
  require(expm)
  
  # Eingangsdaten prüfen
  # data >= 0
  if(any(data<0))stop("negative Daten")
  
  # zeilensummen = 1
  rows <- sapply(1:NROW(data),function(i)sum(data[i,])!=1)  
  if(any(rows))stop("Zeilensumme falsch") 
  
  # Strafsystem
    rek <- 1
    punish <- 10
  
  # Konsturktion eines InitialParameters der einfach nur die NB erfüllt
  initpara <- matrix(0,NROW(data),NCOL(data))
  
  # Initalisierung mit negaitven werten in M
  # Erfüllung der Diagonalbedingung
  for(i in 1:NROW(initpara)){
    if(i>1){
      for(j in ((i-1):1)){
        initpara[i,j] <- min(data[i,j+1]+initpara[i,j+1]-data[i,j],0)
      }
    }
    if(i < NROW(initpara)){      
      for(j in ((i+1):NROW(initpara))){
        initpara[i,j] <- min(data[i,j-1]+initpara[i,j-1]-data[i,j],0)
      }      
    }
  }
  # Spaltenbedingung erfüllen
  s <- NCOL(data)
  for(i in 1:(NROW(initpara)-1)){
    diff <- data[i+1,s] + initpara[i+1,s] - data[i,s] - initpara[i,s]
    initpara[i,s] <- ifelse(diff>=0,initpara[i,s],initpara[i,s]+diff)
  };rm(s,i,j)
  
  IP <- list(initpara)
  
  # Konsturktion eines InitialParameters der einfach nur die NB erfüllt
  initpara <- matrix(0,NROW(data),NCOL(data))  
  # Initalisierung mit positiven werten in M
  # Spaltenbedingung erfüllen
  s <- NCOL(data)
  for(i in (2:NROW(initpara))){
    diff <- data[i,s] + initpara[i,s] - data[i-1,s] - initpara[i-1,s]
    initpara[i,s] <- ifelse(diff>=0,initpara[i,s],initpara[i,s]-diff)
  }
  # Erfüllung der Diagonalbedingung
  for(i in 1:NROW(initpara)){    
    if(i < NROW(initpara)){      
      for(j in ((NROW(initpara)-1):i)){
        diff <- data[i,j]+initpara[i,j]-(data[i,j+1]+initpara[i,j+1])
        initpara[i,j] <- initpara[i,j] + ifelse(diff>=0,0,-diff)
      }      
    }
    if(i>1){
      for(j in (2:i)){
        diff <- data[i,j]+initpara[i,j]-(data[i,j-1]+initpara[i,j-1])
        initpara[i,j] <- initpara[i,j] + ifelse(diff>=0,0,-diff)
      }
    }    
  };rm(s,i,j)
  
  IP <- c(IP,list(initpara));rm(initpara)
  
  # lapply(IP,function(initpara){
    # # Anzeige des Initialparameters
    # print(initpara)
    # # Anzeige von von Dtilde
    # print(data+initpara)
  # })
  
  # Zielfunktion als Bytecode
  quaddist <- cmpfun(function(Dtilde){
    sum(sapply(1:NROW(Dtilde),function(i)sapply(1:NCOL(Dtilde),function(j) sum(
      sapply(2:N,function(n)
        # Angabe der Summenfunktion
        (( (Dtilde%^%n)[i,j] - (data%^%n)[i,j])/( (data%^%n)[i,j] ))^2
      )))))
  })
  
  # Funktion für optim als Bytecode
  minfun <- cmpfun(function(para){
        
    # Starfe verschärfen
    rek <<-rek+1
    
    # D2 = Dtilde
    D2 <- data + matrix(c(para),ncol = NCOL(data))
    
    # NB das Parameters
    # Bedingung an letzte Spalte in D2
    if(any(sapply(1:(NROW(D2)-1),function(j) D2[j,NCOL(D2)]>D2[j+1,NCOL(D2)]))){
      # print("Spaltenbeding. verletzt")
      return(quaddist(D2)+punish^rek)
    }
    # Ordnungs in D2 von der Diagonale aus
    crit <- sapply(1:NROW(D2),function(i){
      erg <- c()
      if(i > 1) erg <- c(erg,sapply(1:(i-1),function(j) D2[i,j]<=D2[i,j+1]))
      if(i < NROW(D2))erg <- c(erg,sapply(i:(NCOL(D2)-1),function(j) D2[i,j]>=D2[i,j+1]))
      any(!erg)
    })
    if(any(crit)){
      # print("Diagbed. verletzt")
      return(quaddist(D2)+punish^rek)
    }
    
    # Funktionswert selbst
    # print("Kein Fehler")
    return(quaddist(D2))
    
  })
  
  
  erg <- lapply(IP,function(initpara){
    erg<-optim(par = c(initpara),fn = minfun)
    rek <<- 1
    erg
  })
  # print(erg)
  index<-which.min(sapply(erg,function(e)e$value))
  erg <- matrix(erg[[index]]$par,NROW(data))
  
  # print(data + erg)
  
  erg

  
}

Test <- matrix(c(0.3,0.1,0.7,0.9),2,2)
print(psoptprob(data = Test,N = 10))

Test <- matrix(runif(16),4)
Test <- matrix(rexp(25),5)
# um die Zeilensumme richtig abzusichern ... blödes R
for(i in 1:NROW(Test)){
  Test <- lapply(1:NROW(Test),function(i)Test[i,]/sum(Test[i,]))
  Test <- do.call(rbind,Test)
}
print(psoptprob(data = Test,N = 2))
 

also die beiden Nebenbedingungen, dass Zeilensumme=1 gelten muss sowie alle Einträge größer 0 gelten für D und Dtilde. Nun kommen die restlichen bereits erwähnten Ungleichheitsnebenbedingungen ins Spiel. D stellt eine Transition Matrix dar, also es werden quasi statistische Daten erhoben und zwar die Verschiebung von Ratingstufen, zum Beispiel sowas wie von AAA nach AA. Da statistische Daten aber durch eine geringe Samplesize nicht den stochastischen Gesetzen entsprechen müssen, wird quasi über diese neuen Nebenbedingungen die entsprechende stochastik intuitiv angepasst, die Ergebnisse unterliegen dann den Gesetzen.