2016-05-09 1 views
2

Je voudrais calculer l'intégrale double suivante, avec la borne inférieure = -Inf et la borne supérieure = Inf pour les deux intégrales. Comment puis-je définir la valeur de la fonction à zéro pour M = 0 et/ou omettre l'intégration numérique sur M = 0? Il s'agit d'une fonction de densité et je veux définir la zone sur M = 0 comme zéro lors du calcul des probabilités.Comment définir une valeur de fonction à un point pour une intégrale dans R?

La fonction peut être vu ici:

enter image description here

Le code R est ici, où l'on peut voir le "saut" à cause de mon problème:

mum<-5.16 
mub<-1.5 
mur<-2.764 
sm<-3.37 
sb<-0.2 
sr<-0.056 

dk<-function(k){ 
integrate(function(r) { 
sapply(r, function(r) { 
1/(2*pi*sr^2)^0.5*exp(-(r-mur)^2/(2*sr^2))*integrate(function(m) {1/(abs(m)*2*pi*sm*sb)*exp(-(m-mum)^2/(2*sm^2)-((k-r)/m-mub)^2/(2*sb^2))},lower=Inf, upper=Inf)$value 
}) 
}, lower=Inf, upper=Inf)$value 

} 

curve(sapply(x,dk), from=-2, to=18, lwd=2,col="red") 
+0

votre code ne fonctionne pas comme cela est. Il y a un '}' inattendu. En outre, aucune donnée fournie pour tester la fonction pour votre comportement souhaité. Pourriez-vous modifier la question pour résoudre ces problèmes. Je m'attends à remplacer votre fonction intégrée par la fonction (m) { si (m == 0) 0 sinon 1/abs (m) * exp (- (m-m3)^2 - ((kr)/m- m2)^2) } va résoudre votre problème, mais besoin de données pour être sûr – dww

+0

S'il vous plaît voir http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example comment pour créer un MWE. – dww

+0

Merci pour votre contribution, j'ai corrigé les crochets. Cependant, je n'ai pas de données à adapter, ce sont juste des paramètres que je dois définir. La fonction que j'utilise en réel et le code est beaucoup plus long. Pour présenter le problème de la division par le jeu, je l'ai simplifié. C'est une fonction de distribution continue et je peux définir la zone sur le point comme zéro, que je dois implémenter dans R. – Seb

Répondre

2

Une approche serait simplement ne pas évaluer votre fonction dk à des valeurs où sa valeur est invalide comme dans par exemple

mum<-5.16 
mub<-1.5 
mur<-2.764 
sm<-3.37 
sb<-0.2 
sr<-0.056 

    dk<-function(k){ 
    integrate(function(r) { 
     sapply(r, function(r) { 
     1/(2*pi*sr^2)^0.5*exp(-(r-mur)^2/(2*sr^2))*integrate(function(m) {1/(abs(m)*2*pi*sm*sb)*exp(-(m-mum)^2/(2*sm^2)-((k-r)/m-mub)^2/(2*sb^2))}, lower=Inf, upper=Inf)$value 
     }) 
    }, lower=Inf, upper=Inf)$value 
    } 

x <- c(seq(-2,-0.1, by=0.2), seq(0.1,18,by=0.2)) 
y <- sapply(x,dk) 
plot(y~x, type="l", col="red") 

enter image description here

L'autre façon est d'avoir votre fonction NA dk retour à des valeurs pour lesquelles la fonction échoue comme par exemple

dk<-function(k){ 
    if ((k-mur)<0.4 & (k-mur)>-0.4) NA 
    else { 
    integrate(function(r) { 
     sapply(r, function(r) { 
     1/(2*pi*sr^2)^0.5*exp(-(r-mur)^2/(2*sr^2))*integrate(function(m) {1/(abs(m)*2*pi*sm*sb)*exp(-(m-mum)^2/(2*sm^2)-((k-r)/m-mub)^2/(2*sb^2))}, lower=Inf, upper=Inf)$value 
     }) 
    }, lower=Inf, upper=Inf)$value 
    } 
    } 

curve(sapply(x,dk), from=-2, to=18, lwd=2,col="red") 

enter image description here