pro mc

  ; input variables
  a_mean = 10.0  &  a_err = 1.0
  ;a_mean = 1.0  &  a_err = 0.5
  b_mean = 1.0  &  b_err = 0.1
  
  ; derived variable
  c_mean = a_mean^2 / b_mean^2
  
  ; standard (Gaussian) error propagation
  c_err = c_mean * sqrt( (2d0*a_err/a_mean)^2 + (2d0*b_err/b_mean)^2 )

  ; print result of c with standard error propagation
  print, 'c = ', c_mean, ' +/- ', c_err 

  ; number of random numbers
  n = long(1d6)

  ; generate Gaussian random numbers based on a seed
  seed = 1
  x1 = randomn(seed+111, n, /double)
  x2 = randomn(seed+222, n, /double)

  ; define variables with errors
  a = a_mean + x1*a_err
  b = b_mean + x2*b_err

  ; plot PDF of a
  plot_pdf, a, XTITLE='a', YTITLE='PDF(a)'

  ; plot PDF of b
  plot_pdf, b, XTITLE='b', YTITLE='PDF(b)'

  ; now define c
  c = a^2/b^2

  ; plot PDF of c
  plot_pdf, c, XTITLE='c', YTITLE='PDF(c)', PDF=pdf, CDF=cdf, LOC=loc, BINSIZE=binsize

  ; get some statistics of c
  
  ; get the mode (most probable value)
  index = where(pdf eq max(pdf))
  mode_c = loc[index]
  print, 'mode_c (from the PDF) = ', mode_c

  ; get the mean
  mean_c = mean(c,/double)
  print, 'mean_c (from the data directly) = ', mean_c
  mean_c = total(pdf*loc)*binsize
  print, 'mean_c (from the PDF of c) = ', mean_c

  ; compute the standard deviation from the PDF
  stddev_c = sqrt ( total(pdf * (loc-mean_c)^2 )*binsize )
  print, 'stddev_c (from the PDF of c) = ', stddev_c

  ; alternative way of computing the standard deviation
  mean_squared = total(pdf*loc^2)*binsize
  print, 'stddev_c (from the PDF of c; alternative) = ', sqrt ( mean_squared - mean_c^2 )

  ; now compute the 16th and 84th percentile
  index = where(cdf ge 0.16)
  sixteenth = loc[index[0]]
  index = where(cdf ge 0.84)
  eightyfourth = loc[index[0]]
  
  print, '16th and 84th percentiles of c = ', sixteenth, eightyfourth

  ; compute the median
  index = where(cdf ge 0.5)
  median_c = loc[index[0]]
  print, 'median_c = (from the CDF of c) = ', median_c

  ; now pick your preferred best value; here let's pick the median
  err_low = median_c - sixteenth
  err_high = eightyfourth - median_c

  ; print the final c
  print, 'c = ' + string(median_c,format='(F5.1)') + ' (' + string(stddev_c,format='(F4.1)') + ') -' $
                + string(err_low,format='(F4.1)') + ' +' + string(err_high,format='(F4.1)')


  stop

end

; procedure to plot the PDF
pro plot_pdf, variable, _EXTRA=extra, PDF=pdf, CDF=cdf, LOC=loc, BINSIZE=binsize

  ; call out previous PDF procedure
  get_pdf, variable, pdf=pdf, loc=loc, cdf=cdf, binsize=binsize, nbins=1000

  ; plot the PDF
  window, /free
  plot, loc, pdf, chars=2, psym=10, _EXTRA=extra

  return

end