Re: [R] Truncated normal distribution

From: Ben Bolker <bolker_at_ufl.edu>
Date: Sun, 16 Dec 2007 17:14:22 -0800 (PST)

ravirrrr wrote:
>
> I have the following code, where we need to solve for mu and sigma, when
> we have mut and sdt. Can you suggest how to use a solve function in R to
> do that? I am new to R and am not sure how to go from defining the
> functions, to solving for them.
>
> Thanks
>
>
> truncated <- function(x)
> {
>
> mu=x[1];
> sigma=x[2];
>
>
> f <- function(x) (1/(sigma*sqrt(2*pi)))*exp(-(x-mu)^2/(2*sigma^2));
>
> pdf.fun <- function(x) x*f(x);
>
> sd.fun <- function(x) (x)^2*f(x);
>
> st=integrate(sd.fun,lower=-Inf,upper=1)$value;
>
> a=integrate(pdf.fun,lower=-Inf,upper=1)$value;
>
> a1=integrate(f,lower=-Inf,upper=1)$value;
>
> mut <- a/a1;
> sdt <- sqrt((st/a1)-(a/a1)^2);
>
> }
>

truncated <- function(x)
{
  mu <- x[1]
  sigma <- x[2]
  pdf.fun <- function(x) x*dnorm(x,mu,sigma)

  sd.fun <- function(x) x^2*dnorm(x,mu,sigma)
  st <- integrate(sd.fun,lower=-Inf,upper=1)$value;
  a <- integrate(pdf.fun,lower=-Inf,upper=1)$value;
  a1 <- pnorm(1,mu,sigma)
  mut <- a/a1;
  sdt <- sqrt((st/a1)-(a/a1)^2)
  c(mut,sdt)
}

truncated(c(0,1)) ## sensible: mean <0, sd<1 truncated(c(0,0.1)) ## sensible: approx 0,1 ## trouble for small values
truncated(c(0,0.001))
truncated(c(0,0.0001))

optfun <- function(p,target=c(0,1)) {
  sum((truncated(p)-target)^2)
}

target <- c(mu=-0.5,sd=2)
fit1 <- optim(fn=optfun,

      par=c(-0.5,2),
      target=target)

fit1

truncated(fit1$par) ## didn't succeed

## let's do something easier -- can we
## work backward to a known value?

t1 <- truncated(c(0,1))

optim(fn=optfun,

      par=c(0,1),
      target=t1)
          

optim(fn=optfun,
      par=c(0.5,2),
      target=t1)

-- 
View this message in context: http://www.nabble.com/Truncated-normal-distribution-tp14348951p14368641.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Mon 17 Dec 2007 - 01:17:56 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 17 Dec 2007 - 01:30:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.