Re: [R] looping problem

From: <Richard.Cotton_at_hsl.gov.uk>
Date: Mon, 14 Apr 2008 10:27:22 +0100


> I would like to do looping for this process below to estimate alpha
> beta from gamma distribution:
> Here are my data:
> day_data1 <-
> 1 2 3 4 5 6 7 8 9 10 11 12
> 13 14 15 16 17 18 19 20 21 22 23
> 1943 48.3 18.5 0.0 0.0 18.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0
> 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.8 2.8
> 1944 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
> 0.0 0.0 0.0 0.0 6.6 0.0 0.0 0.0 0.0 0.0 0.0
> 1945 5.3 0.0 0.0 0.0 0.0 2.5 0.0 0.0 0.0 0.0 0.0 0.0
> 0.0 0.0 0.0 0.0 5.3 0.0 0.0 0.0 0.0 0.0 0.0
> 1946 0.0 0.0 0.0 0.0 0.0 2.3 0.0 0.0 0.0 0.0 4.8 0.3
> 1.5 0.0 0.8 0.0 0.0 5.8 70.6 12.4 0.5 23.6 0.0
> 1947 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
> 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.8 0.0 0.0 0.0
> 1948 0.3 20.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.5 0.5
> 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 

> library(MASS)
> ## Example January Charleville 1943-2007
>
> Here is my code:
>
> alpha.beta <- function(data,n)
> { tol <- 1E-6
> { for (i in 1:n)
> xi <- data[,i]
> ai <- xi [xi > tol]
> fit <- fitdistr(ai,dgamma, list(shape = 1, rate = 0.1), lower =
0.01)
> }
> fit
> }
>
> I?m not sure what went wrong since it gives only one output by right
> 31 outputs.

A few points about your code
1. The opening brace for the for loop is in the wrong place. It should be for (i in 1:n)
{

2. In each iteration of this for loop, you overwrite the value of fit. To assign different values to fit in the loop, you could make it a list, e.g. alpha.beta <- function(data,n)
{
  tol <- 1E-6
  fit <- list()
  for (i in 1:n)
  {
    xi <- data[,i]
    ai <- xi [xi > tol]
    fit[[i]] <- fitdistr(ai,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
  }
  fit
}

3. You could tidy up the code a lot by getting rid of the explicit for loop, and using apply instead.
tol <- 1e-6
myfitdistr <- function(x)
{
  x <- x[x > tol]
  if(length(x)==0) fit <- NULL else fit <- fitdistr(x,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
}
apply(day_data, 2, myfitdistr)

(You'll need to play about with the parameters given to fitdistr, since this code curerntly fails on column 16, because fitdistr can't optimise.)

Regards,
Richie.

Mathematical Sciences Unit
HSL



ATTENTION: This message contains privileged and confidential inform...{{dropped:21}}

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 14 Apr 2008 - 11:04:40 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 14 Apr 2008 - 11:30:29 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.

list of date sections of archive