0

I am trying to fit my data to a gaussian function using the following codes in R. However, it generates a function with maximum value reaching 1.1 while I want the cap it at 1 because the y-value is representing the percentage of correct responses and it should not be higher than 1. How can I set a maximum value of this graph to be 1?

#Generatng Gaussian Function
  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)
 }

#Modeling the Neutral responses
NeutralModel1 <- fitG(FullBehData[which(FullBehData$Respond=="Neutral"),1],FullBehData[which(FullBehData$Respond=="Neutral"),2],50,10,.2)

#Plotting the graph
p3 <- ggplot(data = FullBehData[which(FullBehData$Respond=="Neutral"),],aes(x=AngryLevel,y=Control,colour="Control"))+geom_point()


xseq <- seq(0,100, len = 1000)

yseq1 <- NeutralModel1$par[3]*dnorm(xseq,NeutralModel1$par[1],NeutralModel1$par[2])
curve1 <- data.frame(xseq,yseq1)
curve1$Block <- as.character(ifelse(1,BlockLevels[1],0))
p3 <- p3 + geom_line(data = curve1, aes(x = xseq, y = yseq1,color = Block))
p3

Result

Ar Sou
  • 1
  • 2
    The answer is simply dnorm(x, m, s)/dnorm(m, m, s), but why do you want it to be 1 ..? – Tim Oct 26 '16 at 13:32
  • 2
    If you would tell us more about the data we would be able to give you appropriate answers. As an example, if your data are proportions you could use a GLM as described at http://stats.stackexchange.com/questions/70870. For general-purpose curve fitting, use the solution at http://stats.stackexchange.com/a/70184/919 by fixing $b=0$ and $a=\sqrt{2\pi}$ and varying only $m$ and $s$. – whuber Oct 26 '16 at 14:30
  • Hi, thank you for the replies. I want to cap the maximum value at 1 because the y value is actually representing percentage of responses and it should not be higher than 1. – Ar Sou Nov 06 '16 at 03:49

1 Answers1

1

The issue seems to be that the scale parameter is unrestricted. That allows the fitted function to have an arbitrary high maximum. The maximum of the Gaussian is:
$\frac{scale}{\sigma\sqrt{2\cdot \pi}}$
So if we want the maximum to be at most 1, we want the scale to be in the range:
$0\leq scale \leq \sigma\sqrt{2\cdot \pi}$
This is of course assuming the Gaussian is positive. An easy way to force the scale to be in this range is by using a sigmoid function. That is a nice smooth function which forces the scale to be between 0 and 1. So, one way to rewrite your objective function would be:

fitG =  
function(x,y,mu,sig,scale){

sigmoid = function(x) { 1 / (1 + exp(-x)) }
f = function(p){
d = sqrt(2*pi) * p[2] * sigmoid(p[3])*dnorm(x,mean=p[1],sd=p[2])
sum((d-y)^2)
}
optim(c(mu,sig,scale),f)  
}
dimpol
  • 1,032
  • 5
  • 13
  • Using the sigmoid changes the problem and therefore leads to the wrong objective function. A much simpler--and correct--way to resolve the issue is to parameterize the Gaussian so it always has a unit vertical scale, as in $$f(x;\mu,\sigma)=\exp(-((x-\mu)/\sigma)^2).$$ – whuber Oct 27 '16 at 14:51
  • Well, it depends on what OP means with a Gaussian 'capped at 1'. Does she/he want a Gaussian with a maximum of exactly 1, or a maximum of at most 1? That determines which is the 'wrong objective function' – dimpol Oct 28 '16 at 07:36
  • Sorry; I relied on your prose rather than your code. Your prose suggests you are trying to fit a "sigmoid" rather than fitting a Gaussian. The code reveals you are instead constraining the amplitude of the Gaussian by reparameterizing it. You might want to clarify that in your text. – whuber Oct 28 '16 at 13:33
  • Hi, thank you for all of your responses. Whuber, can you be more specific on how to change the code? – Ar Sou Nov 06 '16 at 05:32