Discussion :

[CupOfAppleTea]

Bonjour,

Je dispose de données auxquelles je souhaiterais appliquer une Gaussienne. Cependant, mes données ne décrivent pas la distribution dans son ensemble. Je m'explique:
Je peux disposer de données qui ne vont décrire que la moitié de la courbe de Gauss (cf.graphe en pièce jointe). Je voudrais donc, pardon si je ne m'exprime pas correctement, fitter une gaussienne (ce qui va je pense nécessiter d'imposer quelques contraintes à la fonction de fit puisque dans mon cas la Gaussienne est en quelques sortes "implicite/cachée" pour moitié) pour ensuite récupérer la valeur "x" pour laquelle "y" prend sa valeur maximale (si je ne fais pas d'erreur, il s'agit de mu, l'un des paramètres de la loi normale). Je précise que l'échelle des "y" est en pourcentage (de 0 à 1).

Auriez-vous une idée de la méthode que je pourrais employer ici?
Merci pour votre aide.

[CupOfAppleTea]

Je précise que la relation que j'étudie est très couramment caractérisée par une gaussienne (phénomène biologique). C'est la raison pour laquelle je me permets d'utiliser cette méthode, bien que n'ayant que la moitié de cette courbe théorique au travers de mes points expérimentaux.

[CupOfAppleTea]

Voilà ce que j'ai fait pour l'instant! Avec des estimateurs fourni par l'utilisateur, ce script permet via la fonction optim de déterminer un couple (mu; sigma).

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
a=read.table("marsillargues.txt", header=TRUE)
 
fitG =
  function(x,y,mu,sig,scale){
 
    f = function(p){
      d = p[3]*dnorm(x,mean=p[1],sd=p[2])
      sum((d-y)^2)
    }
 
    optim(c(mu,sig,scale),f)
  }
 
fit2P = fitG(a$Tleaf,a$CondProp,25.02,14.5,0.2)

J'obtiens les résultats suivants:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
> fit2P
$par
[1] 24.57984 10.09999 24.98141
 
$value
[1] 0.0540938
 
$counts
function gradient 
     114       NA 
 
$convergence
[1] 0
 
$message
NULL

24.98 semble tout à fait cohérent avec mes données. Maintenant, ça coince lorsque j'essaie de ploter tout ça (points expé + fit):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
plot(a$Tleaf,a$CondProp)
p = fit2P$par
lines(a$Tleaf,p[3]*dnorm(a$Tleaf,p[1],p[2]))

Voici le graphe que j'obtiens en sortie:
https://drive.google.com/file/d/0Bwe...it?usp=sharing
Je souhaiterais corriger le fait qu'il y ait 4 graphes superposés.
Egalement, je souhaiterais affiner l'échelle de mon fit.

Dites moi ce que vous pensez de tout ça!
Merci.