6
$\begingroup$

The function dnorm(x) in R gives you the value of the probability density function in the points x of a certain normal distribution (mean = 0 and SD = 1 by default), returning a vector of the same length than x.

However, I want to do the opposite: given a vector that approximates a probability density function (like the result of dnorm), I want to get the mean and standard deviation of the normal distribution represented by that given probability density. Code example of what I would like to do:

pdf = dnorm(seq(-3,3,0.1), mean = 0, sd=1)
## get_normal would be supposes to return an list/vector containing the mean and SD, which in this particular case should be close to 0 and 1 respectively.
var_parameters = get_normal(pdf)
$\endgroup$
13
  • $\begingroup$ Are the evaluation points known? If yes, then it's fairly trivial. If not, there is not much you can do. $\endgroup$ Jul 18, 2021 at 19:47
  • $\begingroup$ And could you elaborate on what you mean by "a vector that approximate a probability density function"? $\endgroup$ Jul 18, 2021 at 19:51
  • $\begingroup$ Check out "Riemann sums". $\endgroup$ Jul 18, 2021 at 19:54
  • $\begingroup$ It sounds like you're asking about a maximum-likelihood estimator. $\endgroup$
    – Sycorax
    Jul 18, 2021 at 19:59
  • $\begingroup$ @jarle-tufto Yes, the evaluation points are known. If you take a look at the code, pdf is a vector of evaluation points (in this specific case, the evaluation points of the density function f(x) of gaussian of mean = 0 and sd = 1, with x taken values from -3 to 3 with a step of 0.1). I wrote "a vector that approximate a probability density function" because the R vector pdf is discrete, and thus it is an approximation of the actual density function (whose parameters I don't know and want to find) If anything is still unclear, just tell me. Thanks! $\endgroup$
    – BN-stats
    Jul 18, 2021 at 20:05

2 Answers 2

3
$\begingroup$

For a normal density function $f,$ if you have a grid of points X and corresponding density values $y = f(x),$ then you can use numerical integration to find $\mu$ and $\sigma.$ [See Note (2) at the end.]

If you have many realizations $X_i$ from the distribution, you can estimate the population mean $\mu$ by the sample mean $\bar X$ and the population SD $\sigma$ by the sample SD $S.$

Another possibility is to use a kernel density estimator (KDE) of $f$ based on a sufficiently large sample. In R the procedure density gives points $(x, y)$ that can be used to plot a density estimator.

set.seed(718)
x = rnorm(100, 50, 7)
mean(x);  sd(x)
[1] 50.62287
[1] 6.443036

hist(x, prob=T, col="skyblue2")
 rug(x);  lines(density(x), col="red")

enter image description here

In R, the KDE consists of 512 points with values summarized as below:

density(x)

Call:
        density.default(x = x)

Data: x (100 obs.);     Bandwidth 'bw' = 2.309

       x               y            
 Min.   :31.36   Min.   :1.974e-05  
 1st Qu.:41.69   1st Qu.:3.239e-03  
 Median :52.03   Median :2.371e-02  
 Mean   :52.03   Mean   :2.417e-02  
 3rd Qu.:62.36   3rd Qu.:4.378e-02  
 Max.   :72.70   Max.   :5.566e-02  

You can estimate $\mu$ and $\sigma$ corresponding to the KDE as follows:

xx = density(x)$x
yy = density(x)$y              # (xx, yy) is KDE plot point
sum(xx*yy)/sum(yy)
[1] 50.62329                   # aprx pop mean = 50
sum((xx-50.62)^2 * yy)/sum(yy)
[1] 46.42294                   # aprx pop variance = 49
sqrt(sum((xx-50.62)^2 * yy)/sum(yy))
[1] 6.813438                   # aprx pop SD = 7

Because $\bar X$ and $S$ are sufficient statistics for $\mu$ and $\sigma,$ it is hard to imagine that $\hat \mu$ and $\hat \sigma$ re-claimed from a KDE (based on data) would be systematically better than the sample mean $\bar X = 50.62$ and SD $S = 6.44.$ I mention the KDE method because it seems possibly related to your question.

Notes: (1) Of course there are also methods for estimating $\bar X$ and $S$ from a histogram, but they can be very inaccurate for small samples.

(2) Here is a numerical evaluation of $\mu \approx \int_0^{100} x\varphi(x,50,7)\, dx,$ using the sum of areas of 1000 rectangles.

m = 1000
w = (100-0)/m
x = seq(0+w/2, 100-w/2, len=m) 
f = x*dnorm(x, 50, 7)
sum(w*f)
[1] 50   # mu
f2 = (x-50)^2*dnorm(x,50,7)
sum(w*f2)
[1] 49   # sigma^2
$\endgroup$
2
  • 1
    $\begingroup$ Hi! I found your answer the most complete by far. I finally opted by just finding the mean and sd using integration in R or numerical integration (your last piece of code). Both work well, but the function integrate() of R seems to have problem integrating between -Inf and Inf when the mean is big and the sd small, you have to manually tune the integration interval (I set mean-5*sd and mean+5*sd). KDE seems more suitable when handled a sample instead of a grid of evaluation points of the function. Still, much appreciated! $\endgroup$
    – BN-stats
    Jul 18, 2021 at 21:56
  • $\begingroup$ Numerical integration for $\sigma^2$ added to Note (2). $\endgroup$
    – BruceET
    Jul 18, 2021 at 22:16
4
$\begingroup$

A very simple, general-purpose solution: First, write a function that takes parameters as an input, and returns the different between the predicted PDF for those parameters and the actual PDF (I've used the sum of squared differences here). Then, use optim() to find the parameters than minimise this function.

x = seq(-3,3,0.1)
pdf = dnorm(x, mean = -.5, sd = .2)
f = function(pars){
  pred_pdf = dnorm(x, mean = pars[1], sd = pars[2])
  err = sum((pdf - pred_pdf)^2)
}
result = optim(c(0, 1), f) # c(0, 1) are initial values
round(result$par, 3)
# [1] -0.5  0.2
$\endgroup$
1
  • 1
    $\begingroup$ This answer is particularly interesting to me, since I'm doing this in relation with an optimization and Machine Learning project! Specifically, I do this to find the parameters of a geometrical mean of two normal density function. I think I will start using the solution proposed by @BruceET, but your seems perfectly valid for my case as well $\endgroup$
    – BN-stats
    Jul 18, 2021 at 22:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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