0

Suppose somebody draw me an histogram and I want to smooth it, and get the smoothed function. Is their a way to do so in R? (The histogram is not coming from data, so kernel density estimators don't seem adapted. Please tell me if you think I am wrong on this.)

So far, I choose to fit a parametric distribution to my histogram. To do so I minimize the integrated square error between my histogram and a beta distribution. Here is my code, where h is a piece-wise constant function with support [0;1].

h<-function(x) (x>0 & x<1)*1

    fit.beta<-function(h){
      dist<-function(alpha,beta){
            diff2<-function(x)(h(x)-dbeta(x,alpha,beta))^2          
            return(integrate(diff2,0,1))
      }
      res<-constrOptim(theta = c(1,1), f = dist,grad=NULL, ui = matrix(c(1,1),1,2), ci = c(0,0)) 
      return<-res
    }

And R says:

 Error in dbeta(x, alpha, beta) : 
  argument "beta" is missing, with no default

I don't understand why R doesn't understand dbeta(x, alpha, beta). I also tryed with dbeta(x, shape1=alpha,shape2=beta) it doesn't work. Could you help me?

Tetsujin no Oni
  • 7,300
  • 2
  • 29
  • 46
Mothas
  • 1
  • 1
  • I bet your code doesn't work because of your constraint matrix `ui` (your argument `ui` should be a 2x2 matrix, see `?constrOptim`). I can not help you further without a reproducible example... Otherwise, I would advise you to use the `fitdistr` package rather than trying to do fitting from scratch, see http://stat.ethz.ch/R-manual/R-patched/library/MASS/html/fitdistr.html :) – Jealie Apr 22 '14 at 14:51
  • The fitdistr function doesn't seem adapted, it needs data to compute the MLE estimate. – Mothas Apr 23 '14 at 12:15
  • that's right, I mis-understood your problem. A "quick and dirty" approach would be to generate data from your histogram (say, `10000*h[i]`) for each bin `i`, and then to call functions from the `fitdistr` package. – Jealie Apr 23 '14 at 15:15

1 Answers1

0

I found the solution to the syntax problem. The constrOptim function only optimize the first argument so it works if the optimized function as only one argument.

   fit.norm<-function(h){
  dist<-function(ab){
    diff2<-function(x)(h(x)-dnorm(x,ab[1], ab[2]))^2
    return(integrate(diff2,0,1)$value)
  }
  res<-constrOptim(theta = c(0.5,1), f = dist,grad=NULL, ui = rbind(c(1,0),c(-1,0),c(0,1)), ci = c(0,-1,0)) 
  return<-list(res)
}
Mothas
  • 1
  • 1