The R Project SVN R

Rev

Rev 42302 | Rev 64611 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 42302 Rev 42303
Line -... Line 1...
-
 
1
/*
-
 
2
 *  Mathlib : A C Library of Special Functions
-
 
3
 *  Copyright (C) 1998 Ross Ihaka
-
 
4
 *  Copyright (C) 2000	    The R Development Core Team
-
 
5
 *  Copyright (C) 2003	    The R Foundation
-
 
6
 *
-
 
7
 *  This program is free software; you can redistribute it and/or modify
-
 
8
 *  it under the terms of the GNU General Public License as published by
-
 
9
 *  the Free Software Foundation; either version 2 of the License, or
-
 
10
 *  (at your option) any later version.
-
 
11
 *
-
 
12
 *  This program is distributed in the hope that it will be useful,
-
 
13
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
-
 
14
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-
 
15
 *  GNU General Public License for more details.
-
 
16
 *
-
 
17
 *  You should have received a copy of the GNU General Public License
-
 
18
 *  along with this program; if not, a copy is available at
-
 
19
 *  http://www.r-project.org/Licenses/
-
 
20
 *
-
 
21
 *  SYNOPSIS
-
 
22
 *
-
 
23
 *	double dnorm4(double x, double mu, double sigma, int give_log)
-
 
24
 *	      {dnorm (..) is synonymous and preferred inside R}
-
 
25
 *
-
 
26
 *  DESCRIPTION
-
 
27
 *
-
 
28
 *	Compute the density of the normal distribution.
-
 
29
 */
-
 
30
 
-
 
31
#include "nmath.h"
-
 
32
#include "dpq.h"
-
 
33
 
-
 
34
double dnorm4(double x, double mu, double sigma, int give_log)
-
 
35
{
-
 
36
#ifdef IEEE_754
-
 
37
    if (ISNAN(x) || ISNAN(mu) || ISNAN(sigma))
-
 
38
	return x + mu + sigma;
-
 
39
#endif
-
 
40
    if(!R_FINITE(sigma)) return R_D__0;
-
 
41
    if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */
-
 
42
    if (sigma <= 0) {
-
 
43
	if (sigma < 0) ML_ERR_return_NAN;
-
 
44
	/* sigma == 0 */
-
 
45
	return (x == mu) ? ML_POSINF : R_D__0;
-
 
46
    }
-
 
47
    x = (x - mu) / sigma;
-
 
48
 
-
 
49
    if(!R_FINITE(x)) return R_D__0;
-
 
50
    return (give_log ?
-
 
51
	    -(M_LN_SQRT_2PI  +	0.5 * x * x + log(sigma)) :
-
 
52
	    M_1_SQRT_2PI * exp(-0.5 * x * x)  /	  sigma);
-
 
53
    /* M_1_SQRT_2PI = 1 / sqrt(2 * pi) */
-
 
54
}