# [R] looping problem

Richard.Cotton at hsl.gov.uk Richard.Cotton at hsl.gov.uk
Mon Apr 14 11:27:22 CEST 2008

```> 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.

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

```