Hi R-users,
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
## To extract all the positive values
x1 <- day_data1[,1]
x2 <- day_data1[,2]
x3 <- day_data1[,3]
x4 <- day_data1[,4]
x5 <- day_data1[,5]
tol <- 1E-6
a1 <- x1[x1>tol]
a2 <- x2[x2>tol]
a3 <- x3[x3>tol]
a4 <- x4[x4>tol]
a5 <- x5[x5>tol]
library(MASS)
## Example January Charleville 1943-2007
fitdistr(a1,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
fitdistr(a2,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
fitdistr(a3,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
fitdistr(a4,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
fitdistr(a5,dgamma, list(shape = 1, rate = 0.1), lower = 0.01)
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.
Thank you for your attention.
__________________________________
t spam protection around
http://mail.yahoo.com
> 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}}
Dear r-expert,
I would like to generate 30 sets of random numbers from uniform distribution
(0,1). Each set of random numbers should have 90 data.? Here is my code:
rand.no <- function(n,itr)
{ for (i in 1:itr)
? {rand.1 <- runif(n,0,1)
? if (i ==1) rand.2 <- rand.1
? else rand.2 <- cbind(rand.2,rand.1)
? }
? rand.2
}
Question: The code is okay is just that it give me 31 sets instead of 30 sets.
Can anybody any adjustment to my code.
Your help is really appreciated.
Thank you.
____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now.
http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ