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