1
$\begingroup$

I am trying to fit a curve to a histogram that looks roughly like exponential decay. histogram

Since this is roughly similar to exponential decay, I figured a good model was

counts ~ counts[1] * exp(-a * mids)

Which is calculated as

h1 <- hist(data[1,])
#find position that is closest to 10, our cutoff, and save only mids before 
#that (these are now our x values)
mids <- h1$mids[1:which.min(abs(h1$mids-10))]
#find y values
counts <- h1$counts[1:which.min(abs(h1$mids-10))]
mdl <- nls(counts ~ counts[1] * exp(- a * mids))

But when I plot this, I do not get a very good fit:

lines(mids, predict(mdl, list(x = mids)))

after fit

Is the problem related to my model? Or am I just fitting to this histogram incorrectly?

This what the data are like:

> mids
  [1] 1.025 1.075 1.125 1.175 1.225 1.275 1.325 1.375 1.425 1.475 1.525 1.575 1.625 1.675
 [15] 1.725 1.775 1.825 1.875 1.925 1.975 2.025 2.075 2.125 2.175 2.225 2.275 2.325 2.375
 [29] 2.425 2.475 2.525 2.575 2.625 2.675 2.725 2.775 2.825 2.875 2.925 2.975 3.025 3.075
 [43] 3.125 3.175 3.225 3.275 3.325 3.375 3.425 3.475 3.525 3.575 3.625 3.675 3.725 3.775
 [57] 3.825 3.875 3.925 3.975 4.025 4.075 4.125 4.175 4.225 4.275 4.325 4.375 4.425 4.475
 [71] 4.525 4.575 4.625 4.675 4.725 4.775 4.825 4.875 4.925 4.975 5.025 5.075 5.125 5.175
 [85] 5.225 5.275 5.325 5.375 5.425 5.475 5.525 5.575 5.625 5.675 5.725 5.775 5.825 5.875
 [99] 5.925 5.975 6.025 6.075 6.125 6.175 6.225 6.275 6.325 6.375 6.425 6.475 6.525 6.575
[113] 6.625 6.675 6.725 6.775 6.825 6.875 6.925 6.975 7.025 7.075 7.125 7.175 7.225 7.275
[127] 7.325 7.375 7.425 7.475 7.525 7.575 7.625 7.675 7.725 7.775 7.825 7.875 7.925 7.975
[141] 8.025 8.075 8.125 8.175 8.225 8.275 8.325 8.375 8.425 8.475 8.525 8.575 8.625 8.675
[155] 8.725 8.775 8.825 8.875 8.925 8.975 9.025 9.075 9.125 9.175 9.225 9.275 9.325 9.375
[169] 9.425 9.475 9.525 9.575 9.625 9.675 9.725 9.775 9.825 9.875 9.925 9.975
> counts
  [1] 83429     0     0     0     0     0     0     0     0     0     0     0     0     0
 [15]     0     0     0     0     0 43939     0     0     0     0     0     0     0     0
 [29]     0     0     0     0     0     0     0     0     0     0     0 25437     0     0
 [43]     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [57]     0     0     0 14877     0     0     0     0     0     0     0     0     0     0
 [71]     0     0     0     0     0     0     0     0     0  8876     0     0     0     0
 [85]     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [99]     0  5410     0     0     0     0     0     0     0     0     0     0     0     0
[113]     0     0     0     0     0     0     0  3450     0     0     0     0     0     0
[127]     0     0     0     0     0     0     0     0     0     0     0     0     0  2081
[141]     0     0     0     0     0     0     0     0     0     0     0     0     0     0
[155]     0     0     0     0     0  1275     0     0     0     0     0     0     0     0
[169]     0     0     0     0     0     0     0     0     0     0     0   788
$\endgroup$
2
$\begingroup$

Many points to make.

1) you're fitting the zeroes between your actual data (it looks like the data can take only integer values), fitting zeroes where you can't have data is pulling your curve down.

2) you're dealing with a discrete distribution (in fact, insertion length sounds like a count of how many whatevers were inserted in whatever), so the relevant distribution with exponential-like decay is called the geometric.

3) Maximum likelihood estimation is easy in this case. $\hat{p}=1/\bar{x}$. (what's mean(data[1,])?). From that estimate you can then obtain the fitted exponential curve.

4) Even if you do want to directly fit an exponential curve, the problem with nonlinear least squares (even if you do it right by only fitting where the data can occur) is that the counts in each of the bins are not equally variable (you're putting relatively too much weight on observations that are less certain).

One could use a GLM with log-link instead, and take proper account of the multinomial nature of the data. But it won't do any better than just fitting the geometric via maximum likelihood.

$\endgroup$
5
  • $\begingroup$ Thanks. I am confused on two things - 1) via that MLE estimate, p=0.376. But with fitdistr in MASS, I got p=0.273. Second, once I have my fit, how do I scale it so that the values match the data? I am trying plot(dgeom(seq(1,10), p)) which gives a correct looking curve which is obviously on a very different x scale than the data are. $\endgroup$ – Ian Fiddes Sep 4 '14 at 1:41
  • 1
    $\begingroup$ What did you do to use fitsdistr in MASS? Note that your smallest value is 1, while R's geometric distribution starts at 0. (If you look at the Wikipedia link to the geometric distribution it lists two possible parameterizations of the distribution, one for a smallest value of 0, and one for a smallest value of 1. If you confuse them, you're going to have a bad time.) 2) The fitted distribution gives the probability of an observation at each value. Did you try multiplying by the total number of observations? ... (ctd) $\endgroup$ – Glen_b Sep 4 '14 at 1:44
  • 1
    $\begingroup$ (ctd)... Note also that if you make your histogram bins of width $h$ that will scale the height by $1/h$, so you also have to adjust for that effect (but seriously, don't use histograms for discrete data). $\endgroup$ – Glen_b Sep 4 '14 at 1:48
  • 1
    $\begingroup$ Rather than a histogram, try hx=hist(x,breaks=(min(x):(max(x)+1))-0.5,plot=FALSE); with(hx,plot(mids,counts,type="h",lwd=2)) for example. $\endgroup$ – Glen_b Sep 4 '14 at 1:57
  • $\begingroup$ Thanks. I got a good looking fit with creating my own geometric function using the other parameterization function(p, x) {p*(1-p)^(x-1)} then plotting it as lines(g(p, seq(1:10))*length(data[1,])) over my (admittedly bad) use of histograms for discrete data $\endgroup$ – Ian Fiddes Sep 4 '14 at 4:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.