[R] looping problem

[Richard.Cotton at hsl.gov.uk]

> 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}}