Lösen Gleichungssystem mit Integralen

[anna]

Hallo zusammen,

ich möchte ein Gleichungssystem mit 3 Gleichungen nach 3 Variablen (m1, m2, alpha) auflösen. Hat jemand bitte eine Idee? Ich verzweifele

Ich habe die Gleichungen definiert (funktioniert das so?):

Code: Alles auswählen

#define constants
gamma = 0.9
rho = -0.7
d = 1.5
sigma1 = 2.0
sigma2= 3.0
mu1 = 0.5
mu2 = 0.5
d1=2
d2=2

#define functions
my_func1 = function(x,z) 
{1/(2*pi*sigma1(sigma1+sigma+2*rho*sigma1*sigma2*sqrt(1-rho^2)))
  *exp(-(1/(2*(1-rho^2)))*((x-mu1)/sigma1^2 + (x+y-mu1-mu2)^2/(sigma1+sigma2+2*rho*sigma1*sigma2)^2)
       -(2*rho*(x-mu1)*(x+y-mu1-mu2))/sigma1*(sigma1+sigma2+2*rho*sigma1*sigma2))}

my_func2 = function(y,z)
{1/(2*pi*sigma2(sigma1+sigma+2*rho*sigma1*sigma2*sqrt(1-rho^2)))
  *exp(-(1/(2*(1-rho^2)))*((y-mu2)/sigma2^2 + (x+y-mu1-mu2)^2/(sigma1+sigma2+2*rho*sigma1*sigma2)^2)
       -(2*rho*(y-mu2)*(x+y-mu1-mu2))/sigma2*(sigma1+sigma2+2*rho*sigma1*sigma2))}

#define equations
eq1 = function(m1,m2, alpha)
{integrate(my_func1, lower = -Inf, upper = d2-m2+alpha) - integrate(my_func2, lower = -Inf, upper = d1-m1+alpha)}

eq2 = function(m1,m2, alpha)
{(integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower = -Inf, upper= (d1-mu1-m1)/sigma1) 
  + integrate(integrate(my_func1, lower= -Inf, upper= d) lower=-Inf, upper=d1-m1))
  -(integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower = -Inf, upper= (d2-mu2-m2)/sigma2) 
    + integrate(integrate(my_func1, lower= -Inf, upper= d) lower=-Inf, upper=d2-m2)}

eq3 =  function(m1,m2, alpha)
{-simga1/(sqrt(2*pi))*exp(-(d1*mu1-m1)^2/(2*sigma1^2))+(m1-mu1-d1)*integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower=-Inf, upper=(d1-mu1-m1)/sigma1)
  +(m1-d1+alpha)*integrate(integrate(my_func1), lower=-Inf, upper=d)lower=-Inf, upper=d1-m1
  +integrate(integrate(x*my_func1, lower=-Inf, upper=d),lower=d1-m1-alpha, upper=d1-m1)
  +integrate(integrate(my_func2, lower=-Inf, upper=d),lower=-Inf, upper=d2-m2)
  -(m2-alpha-d2)*integrate(integrate(my_func2, lower=-Inf, upper=d),lower=-Inf, upper=d2-m2-alpha)
  -integrate(integrate(y*my_func2, lower=-Inf, upper=d), lower=d2-m2+alpha, upper=d2-m2) - gamma}

Vielen Dank und liebe Grüße
Anna

[jogo]

Hallo Anna,

bisher hast Du die drei Gleichungen nur als Funktionen definiert.
Aber welche Werte sollen die drei Funktionen annehmen?
... nur so wird auch wirklich ein Gleichungssystem daraus.

Gruß, Jörg
p.s.:
Bei Deiner vorherigen Frage viewtopic.php?f=7&t=2038 hast Du die letzten Antworten nicht kommentiert, sondern bist einfach abgetaucht. Das ist kein nettes Verhalten gegenüber den Hilfswilligen.

[anna]

Hallo Jörg,

die Funktionen sind so definiert, dass sie 0 geben sollen. Hast du eine Idee wie ich das machen kann?

Vielen Dank
Anna

[jogo]

ja

[bigben]

Code: Alles auswählen

#define constants
gamma = 0.9
rho = -0.7
d = 1.5
sigma1 = 2.0
sigma2= 3.0
mu1 = 0.5
mu2 = 0.5
d1=2
d2=2

#define functions
my_func1 = function(x,z) 
{1/(2*pi*sigma1(sigma1+sigma+2*rho*sigma1*sigma2*sqrt(1-rho^2)))
  *exp(-(1/(2*(1-rho^2)))*((x-mu1)/sigma1^2 + (x+y-mu1-mu2)^2/(sigma1+sigma2+2*rho*sigma1*sigma2)^2)
       -(2*rho*(x-mu1)*(x+y-mu1-mu2))/sigma1*(sigma1+sigma2+2*rho*sigma1*sigma2))}

my_func2 = function(y,z)
{1/(2*pi*sigma2(sigma1+sigma+2*rho*sigma1*sigma2*sqrt(1-rho^2)))
  *exp(-(1/(2*(1-rho^2)))*((y-mu2)/sigma2^2 + (x+y-mu1-mu2)^2/(sigma1+sigma2+2*rho*sigma1*sigma2)^2)
       -(2*rho*(y-mu2)*(x+y-mu1-mu2))/sigma2*(sigma1+sigma2+2*rho*sigma1*sigma2))}

#define equations
eq1 = function(m1,m2, alpha)
{integrate(my_func1, lower = -Inf, upper = d2-m2+alpha) - integrate(my_func2, lower = -Inf, upper = d1-m1+alpha)}

eq2 = function(m1,m2, alpha)
{(integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower = -Inf, upper= (d1-mu1-m1)/sigma1) 
  + integrate(integrate(my_func1, lower= -Inf, upper= d) lower=-Inf, upper=d1-m1))
  -(integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower = -Inf, upper= (d2-mu2-m2)/sigma2) 
    + integrate(integrate(my_func1, lower= -Inf, upper= d) lower=-Inf, upper=d2-m2)}

eq3 =  function(m1,m2, alpha)
{-simga1/(sqrt(2*pi))*exp(-(d1*mu1-m1)^2/(2*sigma1^2))+(m1-mu1-d1)*integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower=-Inf, upper=(d1-mu1-m1)/sigma1)
  +(m1-d1+alpha)*integrate(integrate(my_func1), lower=-Inf, upper=d)lower=-Inf, upper=d1-m1
  +integrate(integrate(x*my_func1, lower=-Inf, upper=d),lower=d1-m1-alpha, upper=d1-m1)
  +integrate(integrate(my_func2, lower=-Inf, upper=d),lower=-Inf, upper=d2-m2)
  -(m2-alpha-d2)*integrate(integrate(my_func2, lower=-Inf, upper=d),lower=-Inf, upper=d2-m2-alpha)
  -integrate(integrate(y*my_func2, lower=-Inf, upper=d), lower=d2-m2+alpha, upper=d2-m2) - gamma}

Hallo Anna,

also ich bekomme da einen Haufen Syntaxfehler, wenn ich das in R copypaste. Schlage vor, Du postest erstmal lauffähigen Code.

LG,
Bernhard

[anna]

Hallo Bernhard,

entschuldige! Läuft es jetzt?

Code: Alles auswählen

#define constants
gamma = 0.9
rho = -0.7
d = 1.5
sigma1 = 2.0
sigma2= 3.0
mu1 = 0.5
mu2 = 0.5
d1=2
d2=2

#define functions
my_func1 = function(x,z) 
{(1/(2*pi*sigma1*(sigma1+sigma+2*rho*sigma1*sigma2*sqrt(1-rho^2))))*exp(-(1/(2*(1-rho^2)))*((x-mu1)/sigma1^2 + (x+y-mu1-mu2)^2/(sigma1+sigma2+2*rho*sigma1*sigma2)^2)
       -(2*rho*(x-mu1)*(x+y-mu1-mu2))/sigma1*(sigma1+sigma2+2*rho*sigma1*sigma2))}

my_func2 = function(y,z)
{1/(2*pi*sigma2*(sigma1+sigma+2*rho*sigma1*sigma2*sqrt(1-rho^2)))*exp(-(1/(2*(1-rho^2)))*((y-mu2)/sigma2^2 + (x+y-mu1-mu2)^2/(sigma1+sigma2+2*rho*sigma1*sigma2)^2)
       -(2*rho*(y-mu2)*(x+y-mu1-mu2))/sigma2*(sigma1+sigma2+2*rho*sigma1*sigma2))}

#define equations
eq1 = function(m1,m2, alpha)
{integrate(my_func1, lower = -Inf, upper = d2-m2+alpha) - integrate(my_func2, lower = -Inf, upper = d1-m1+alpha)}

eq2 = function(m1,m2, alpha)
{integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower = -Inf, upper= ((d1-mu1-m1)/sigma1)) + integrate(integrate(my_func1, lower= -Inf, upper= d), lower=-Inf, upper=(d1-m1))-integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower = -Inf, upper= ((d2-mu2-m2)/sigma2)) + integrate(integrate(my_func1, lower= -Inf, upper= d), lower=-Inf, upper=(d2-m2))}


eq3 =  function(m1,m2, alpha)
{-simga1/(sqrt(2*pi))*exp(-(d1*mu1-m1)^2/(2*sigma1^2))+(m1-mu1-d1)*integrate(1/(sqrt(2*pi))*exp(-x^2/2), lower=-Inf, upper=((d1-mu1-m1)/sigma1))+(m1-d1+alpha)*integrate(integrate(my_func1, lower=-Inf, upper=d), lower=-Inf, upper=(d1-m1))
  +integrate(integrate(x*my_func1, lower=-Inf, upper=d),lower=(d1-m1-alpha), upper=(d1-m1))+integrate(integrate(my_func2, lower=-Inf, upper=d),lower=-Inf, upper=(d2-m2))-(m2-alpha-d2)*integrate(integrate(my_func2, lower=-Inf, upper=d),lower=-Inf, upper=(d2-m2-alpha))-integrate(integrate(y*my_func2, lower=-Inf, upper=d), lower=(d2-m2+alpha), upper=(d2-m2)) - gamma}

Vielen Dank schonmal
Anna

[anna]

ja

Wärst du so lieb mir zu helfen?

[bigben]

ja

Wärst du so lieb mir zu helfen?

Wie ich Jörg kenne, ist er grundsätzlich bestimmt bereit. Warum Du es ihm aber nicht leicht machst, indem Du wenigstens lauffähigen Code postest, steht in den Sternen. Hast Du Dir die Mühe gemacht, den Code in R mal aufzurufen? Bestimmt nicht. Sonst wärest Du darüber gestolpert, dass es in my_func folgenden Abschnitt gibt: sigma1+sigma+2*rho*s...

Dieses "sigma" ist aber vorher nirgends definiert worden und in eq3 wäre Dir ein Tippfehler wie simga anstelle von sigma bestimmt aufgefallen, hättest Du es mal laufen lassen.

Wie soll man Dir helfen, Funktionen zu untersuchen, wenn Du die Funktion nicht lesbar aufschreibst?

[jogo]

Hallo Anna,

es gibt grundsätzlich zwei Varianten für eine Lösung:
1. als nichtlineares Gleichungssystem
2. als kleinste-Quadrate-Problem (Minimierung einer Summe von Quadraten)

In beiden Fällen hängt das Gelingen von dem Vorhandensein von guten Startwerten ab.
Kannst Du uns passende Startwerte anbieten?

Dein Code ist aber immer noch nicht fehlerfrei:

Code: Alles auswählen

> gesamt <- function(m1,m2, alpha) {
+   eq1(m1,m2, alpha)^2 + eq2(m1,m2, alpha)^2 + eq3(m1,m2, alpha)^2
+ }
> 
> gesamt(0, 0, .05)

Fehler in sigma1 + sigma : nicht-numerisches Argument für binären Operator

Ich habe noch versucht, diesen Fehler zu beseitigen, aber dann tauchen andere Fehler auf. Damit entstehen auch Zweifel, dass selbst wenn Du syntaktisch fehlerfreien Code lieferst, ob der dann auch das macht, was Du willst oder ob wegen fehlerhafter Umsetzung der Code etwas ganz anderes rechnet, als Du willst.
Wahrscheinlich ist es besser, Du beschreibst, was Du mit den Berechnungen erreichen möchtest.
Gibt es irgendein Aufgabenblatt als PDF-Datei oder andere Quellen der Information?

ja

Wärst du so lieb mir zu helfen?

Bist Du bitte so freundlich und hältst Dich in zukünftigen Nachrichten an:
viewtopic.php?f=20&t=29 (Bisher habe ich die Formatierungen eingefügt.)
Gibst Du noch einen Kommentar zu den Antworten bei Deiner vorherigen Frage ab
viewtopic.php?f=7&t=2038

Warum verwendest Du immer noch handgestricktes integrate() über der Dichtefunktion der Normalverteilung, obwohl ich Dir schon bei Deiner vorherigen Frage gezeigt habe, wie Du die in R vorhandenen Funktionen einsetzen kannst, um eine übersichtliche Schreibweise zu erreichen
viewtopic.php?p=9826#p9826
Dein Code sieht so aus, als ob Du es stellenweise mit einer 2-dimensionalen Normalverteilung zu tun hast. Eventuell kannst Du das Paket mvtnorm einsetzen.
https://cran.r-project.org/web/packages ... index.html
https://de.wikipedia.org/wiki/Mehrdimen ... verteilung

So, meine Lust an der Bearbeitung dieser Frage ist erloschen.
Gruß, Jörg