Hello, the more general thing I'd like to learn here is how to compute Function of Data on the basis of grouping determiend by n variables. In terms of the reason why I am interested in this, I need to compute the average of my data based on the value of the month and day across years. I have come up withy the code below which, as far as I can see, does what I need but getting either a more elegant or a more versatile way to do this would be nice. Thanks Paolo. Days = format(as.Date(Data[["GasDays"]], format = "%d/%m/%Y"), "%d") Months = format(as.Date(Data[["GasDays"]], format = "%d/%m/%Y"), "%m") MonthDayCombs = paste(Months, Days) AvgDemand = data.matrix(by(Data$RescaledDemand, DayMonthCombs, mean)) On 4 July 2011 10:34, EdBo <n.bowora@gmail.com> wrote:> Hi > > May you help me correct my loop function. > > I want optim to estimates al_j; au_j; sigma_j; b_j by looking at 0 to 20, > 21 to 40, 41 to 60 data points. > > The final result should have 4 columns of each of the estimates AND 4 rows > of each of 0 to 20, 21 to 40, 41 to 60. > > ###MY code is > > n=20 > runs=4 > out=matrix(0,nrow=runs) > > llik = function(x) > { > al_j=x[1]; au_j=x[2]; sigma_j=x[3]; b_j=x[4] > sum(na.rm=T, > ifelse(a$R_j< 0, -log(1/(2*pi*(sigma_j^2)))- > (1/(2*(sigma_j^2))*(a$R_j+al_j-b_j*a$R_m))^2, > ifelse(a$R_j>0 , -log(1/(2*pi*(sigma_j^2)))- > (1/(2*(sigma_j^2))*(a$R_j+au_j-b_j*a$R_m))^2, > > -log(pnorm(au_j,mean=b_j*a$R_m,sd=sqrt(sigma_j^2))- > pnorm(au_j,mean=b_j*a$R_m,sd=sqrt(sigma_j^2))))) > > ) > > } > > start.par = c(0, 0, 0.01, 1) > out1 = optim(llik, par=start.par, method="Nelder-Mead") > > > for (i in 1: runs) > { > index_start=20*(i-1)+1 > index_end= 20*i > out[i]=out1[index_start:index_end] > } > out > > > Thank you in advance > > Edward > UCT > ####My data > > R_j R_m > -0.0625 0.002320654 > 0 -0.004642807 > 0.033333333 0.005936332 > 0.032258065 0.001060848 > 0 0.007114057 > 0.015625 0.005581558 > 0 0.002974794 > 0.015384615 0.004215271 > 0.060606061 0.005073116 > 0.028571429 -0.006001279 > 0 -0.002789594 > 0.013888889 0.00770633 > 0 0.000371663 > 0.02739726 -0.004224228 > -0.04 0.008362539 > 0 -0.010951605 > 0 0.004682924 > 0.013888889 0.011839993 > -0.01369863 0.004210383 > -0.027777778 -0.04658949 > 0 0.00987272 > -0.057142857 -0.062203157 > -0.03030303 -0.119177639 > 0.09375 0.077054642 > 0 -0.022763619 > -0.057142857 0.050408775 > 0 0.024706076 > -0.03030303 0.004043701 > 0.0625 0.004951088 > 0 -0.005968731 > 0 -0.038292548 > 0 0.013381097 > 0.014705882 0.006424728 > -0.014492754 -0.020115626 > 0 -0.004837891 > -0.029411765 -0.022054654 > 0.03030303 0.008936428 > 0.044117647 8.16925E-05 > 0 -0.004827246 > -0.042253521 0.004653096 > -0.014705882 -0.004222151 > 0.029850746 0.000107267 > -0.028985507 -0.001783206 > 0.029850746 -0.006372981 > 0.014492754 0.005492374 > -0.028571429 -0.009005846 > 0 0.001031683 > 0.044117647 0.002800551 > > > > > > > > > > > > > -- > View this message in context: > http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3643230.html > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > 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<http://www.r-project.org/posting-guide.html> > and provide commented, minimal, self-contained, reproducible code. >[[alternative HTML version deleted]]
Actually to get exactly what I want I need to add
no.dimnames(AvgDemand )
where
no.dimnames <- function(a) {
## Remove all dimension names from an array for compact printing.
d <- list()
l <- 0
for(i in dim(a)) {
d[[l <- l + 1]] <- rep("", i)
}
dimnames(a) <- d
a
}
Thanks
Paolo
On 6 July 2011 11:06, Paolo Rossi <statmailinglists@googlemail.com> wrote:
> Hello,
>
> the more general thing I'd like to learn here is how to compute
Function of
> Data on the basis of grouping determiend by n variables.
>
> In terms of the reason why I am interested in this, I need to compute the
> average of my data based on the value of the month and day across years. I
> have come up withy the code below which, as far as I can see, does what I
> need but getting either a more elegant or a more versatile way to do this
> would be nice.
>
> Thanks
>
> Paolo.
>
> Days = format(as.Date(Data[["GasDays"]], format =
"%d/%m/%Y"), "%d")
> Months = format(as.Date(Data[["GasDays"]], format =
"%d/%m/%Y"), "%m")
> MonthDayCombs = paste(Months, Days)
> AvgDemand = data.matrix(by(Data$RescaledDemand, DayMonthCombs, mean))
>
>
>
> On 4 July 2011 10:34, EdBo <n.bowora@gmail.com> wrote:
>
>> Hi
>>
>> May you help me correct my loop function.
>>
>> I want optim to estimates al_j; au_j; sigma_j; b_j by looking at 0 to
20,
>> 21 to 40, 41 to 60 data points.
>>
>> The final result should have 4 columns of each of the estimates AND 4
rows
>> of each of 0 to 20, 21 to 40, 41 to 60.
>>
>> ###MY code is
>>
>> n=20
>> runs=4
>> out=matrix(0,nrow=runs)
>>
>> llik = function(x)
>> {
>> al_j=x[1]; au_j=x[2]; sigma_j=x[3]; b_j=x[4]
>> sum(na.rm=T,
>> ifelse(a$R_j< 0, -log(1/(2*pi*(sigma_j^2)))-
>> (1/(2*(sigma_j^2))*(a$R_j+al_j-b_j*a$R_m))^2,
>> ifelse(a$R_j>0 , -log(1/(2*pi*(sigma_j^2)))-
>> (1/(2*(sigma_j^2))*(a$R_j+au_j-b_j*a$R_m))^2,
>>
>> -log(pnorm(au_j,mean=b_j*a$R_m,sd=sqrt(sigma_j^2))-
>>
pnorm(au_j,mean=b_j*a$R_m,sd=sqrt(sigma_j^2)))))
>>
>> )
>>
>> }
>>
>> start.par = c(0, 0, 0.01, 1)
>> out1 = optim(llik, par=start.par, method="Nelder-Mead")
>>
>>
>> for (i in 1: runs)
>> {
>> index_start=20*(i-1)+1
>> index_end= 20*i
>> out[i]=out1[index_start:index_end]
>> }
>> out
>>
>>
>> Thank you in advance
>>
>> Edward
>> UCT
>> ####My data
>>
>> R_j R_m
>> -0.0625 0.002320654
>> 0 -0.004642807
>> 0.033333333 0.005936332
>> 0.032258065 0.001060848
>> 0 0.007114057
>> 0.015625 0.005581558
>> 0 0.002974794
>> 0.015384615 0.004215271
>> 0.060606061 0.005073116
>> 0.028571429 -0.006001279
>> 0 -0.002789594
>> 0.013888889 0.00770633
>> 0 0.000371663
>> 0.02739726 -0.004224228
>> -0.04 0.008362539
>> 0 -0.010951605
>> 0 0.004682924
>> 0.013888889 0.011839993
>> -0.01369863 0.004210383
>> -0.027777778 -0.04658949
>> 0 0.00987272
>> -0.057142857 -0.062203157
>> -0.03030303 -0.119177639
>> 0.09375 0.077054642
>> 0 -0.022763619
>> -0.057142857 0.050408775
>> 0 0.024706076
>> -0.03030303 0.004043701
>> 0.0625 0.004951088
>> 0 -0.005968731
>> 0 -0.038292548
>> 0 0.013381097
>> 0.014705882 0.006424728
>> -0.014492754 -0.020115626
>> 0 -0.004837891
>> -0.029411765 -0.022054654
>> 0.03030303 0.008936428
>> 0.044117647 8.16925E-05
>> 0 -0.004827246
>> -0.042253521 0.004653096
>> -0.014705882 -0.004222151
>> 0.029850746 0.000107267
>> -0.028985507 -0.001783206
>> 0.029850746 -0.006372981
>> 0.014492754 0.005492374
>> -0.028571429 -0.009005846
>> 0 0.001031683
>> 0.044117647 0.002800551
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> --
>> View this message in context:
>> http://r.789695.n4.nabble.com/loop-in-optim-tp3643230p3643230.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> ______________________________________________
>> 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<http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
[[alternative HTML version deleted]]