[R] Issue with asin()

2012-03-19 Thread Letnichev
Hello everyone,

I am working for a few days already on a basic algorithm, very common in
applied agronomy, that aims to determine the degree-days necessary for a
given individual to reach a given growth stade. The algorithm (and context)
is explained here:  http://www.oardc.ohio-state.edu/gdd/glossary.htm , and
so I implemented my function in R as follows:

DD - function(Tmin, Tmax, Tseuil, meanT, method = DDsin)
### function that calculates the degree-days based on
### minimum and maximum recorded temperatures and the
### minimal threshold temperature (lower growth temperature)
{
### method arcsin
if(method == DDsin){
cond1 - (Tmax = Tseuil)
cond2 - (Tmin = Tseuil)
amp - ((Tmax - Tmin) / 2)
print((Tseuil-meanT)/amp)
alpha - asin((Tseuil - meanT) / amp)
DD_ifelse3 - ((1 / pi) * ((meanT - Tseuil) * ((pi/2) - alpha)) 
+
amp*cos(alpha))

DD - ifelse(cond1, 0, ifelse(cond2, (meanT - Tseuil), 
DD_ifelse3))
}

### method (Tmin + Tmax) / 2
else if(method == DDt2){
cond1 - (meanT  Tseuil)
DD - ifelse(cond1,(meanT - Tseuil),0)
}

else{
stop(\nMethod name is invalid.\nMethods available = DDsin 
(sinus) or DDt2
(mean)\n)
}
return(DD)
}

BUT! When I try to process random data:

library(reshape2)
library(plyr)

station - rep(c(station1,station2,station3), 20)
values_min - sample(-5:20, size = 60, replace = T)
values_max - sample(20:40, size = 60, replace = T)
meanT - ((values_min+values_max)/2)
d - data.frame(station,values_min,values_max,meanT)
names(d) - c(station, values_min,values_max,meanT)

x-ddply(d, .(station), transform, t1 =
cumsum(DD(values_min,values_max,0,meanT)))

I get a warning on my alpha calculation (NaN produced); indeed, the values I
give as argument to asin() are out of the range [-1:1], as the print()
reveals. I can't figure out how to solve this issue, because the same
algorithm works in Excel (visual basic).
It is very annoying, especially because it seems that no occurence of such
error using that algorithm can be found on Internet.
Any help is welcome :) Thanks for your time

P.

--
View this message in context: 
http://r.789695.n4.nabble.com/Issue-with-asin-tp4484462p4484462.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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Issue with asin()

2012-03-19 Thread Duncan Murdoch

On 12-03-19 7:42 AM, Letnichev wrote:

Hello everyone,

I am working for a few days already on a basic algorithm, very common in
applied agronomy, that aims to determine the degree-days necessary for a
given individual to reach a given growth stade. The algorithm (and context)
is explained here:  http://www.oardc.ohio-state.edu/gdd/glossary.htm , and
so I implemented my function in R as follows:

DD- function(Tmin, Tmax, Tseuil, meanT, method = DDsin)
### function that calculates the degree-days based on
### minimum and maximum recorded temperatures and the
### minimal threshold temperature (lower growth temperature)
{
### method arcsin
if(method == DDsin){
cond1- (Tmax= Tseuil)
cond2- (Tmin= Tseuil)


These look like useful diagnostics of out-of-range values, but you don't 
use them before the arcsin transformation.



amp- ((Tmax - Tmin) / 2)
print((Tseuil-meanT)/amp)
alpha- asin((Tseuil - meanT) / amp)




DD_ifelse3- ((1 / pi) * ((meanT - Tseuil) * ((pi/2) - alpha)) +
amp*cos(alpha))

DD- ifelse(cond1, 0, ifelse(cond2, (meanT - Tseuil), 
DD_ifelse3))
}

### method (Tmin + Tmax) / 2
else if(method == DDt2){
cond1- (meanT  Tseuil)
DD- ifelse(cond1,(meanT - Tseuil),0)
}

else{
stop(\nMethod name is invalid.\nMethods available = DDsin 
(sinus) or DDt2
(mean)\n)
}
return(DD)
}

BUT! When I try to process random data:


It's a good idea to use set.seed when trying to debug problems like 
this.  Then you can construct a reproducible example.  I'd also suggest 
getting rid of ddply at least for debugging; it makes it harder to see 
what's going on.


Duncan Murdoch




library(reshape2)
library(plyr)

station- rep(c(station1,station2,station3), 20)
values_min- sample(-5:20, size = 60, replace = T)
values_max- sample(20:40, size = 60, replace = T)
meanT- ((values_min+values_max)/2)
d- data.frame(station,values_min,values_max,meanT)
names(d)- c(station, values_min,values_max,meanT)

x-ddply(d, .(station), transform, t1 =
cumsum(DD(values_min,values_max,0,meanT)))

I get a warning on my alpha calculation (NaN produced); indeed, the values I
give as argument to asin() are out of the range [-1:1], as the print()
reveals. I can't figure out how to solve this issue, because the same
algorithm works in Excel (visual basic).
It is very annoying, especially because it seems that no occurence of such
error using that algorithm can be found on Internet.
Any help is welcome :) Thanks for your time

P.

--
View this message in context: 
http://r.789695.n4.nabble.com/Issue-with-asin-tp4484462p4484462.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
and provide commented, minimal, self-contained, reproducible code.


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


Re: [R] Issue with asin()

2012-03-19 Thread Berend Hasselman

On 19-03-2012, at 12:42, Letnichev wrote:

 Hello everyone,
 
 I am working for a few days already on a basic algorithm, very common in
 applied agronomy, that aims to determine the degree-days necessary for a
 given individual to reach a given growth stade. The algorithm (and context)
 is explained here:  http://www.oardc.ohio-state.edu/gdd/glossary.htm , and
 so I implemented my function in R as follows:
 
 DD - function(Tmin, Tmax, Tseuil, meanT, method = DDsin)
 ### function that calculates the degree-days based on
 ### minimum and maximum recorded temperatures and the
 ### minimal threshold temperature (lower growth temperature)
   {
   ### method arcsin
   if(method == DDsin){
   cond1 - (Tmax = Tseuil)
   cond2 - (Tmin = Tseuil)
   amp - ((Tmax - Tmin) / 2)
   print((Tseuil-meanT)/amp)
   alpha - asin((Tseuil - meanT) / amp)
   DD_ifelse3 - ((1 / pi) * ((meanT - Tseuil) * ((pi/2) - alpha)) 
 +
 amp*cos(alpha))
   
   DD - ifelse(cond1, 0, ifelse(cond2, (meanT - Tseuil), 
 DD_ifelse3))
   }
 
   ### method (Tmin + Tmax) / 2
   else if(method == DDt2){
   cond1 - (meanT  Tseuil)
   DD - ifelse(cond1,(meanT - Tseuil),0)
   }
 
   else{
   stop(\nMethod name is invalid.\nMethods available = DDsin 
 (sinus) or DDt2
 (mean)\n)
   }
   return(DD)
 }
 
 BUT! When I try to process random data:
 
 library(reshape2)
 library(plyr)
 
 station - rep(c(station1,station2,station3), 20)
 values_min - sample(-5:20, size = 60, replace = T)
 values_max - sample(20:40, size = 60, replace = T)
 meanT - ((values_min+values_max)/2)
 d - data.frame(station,values_min,values_max,meanT)
 names(d) - c(station, values_min,values_max,meanT)
 
 x-ddply(d, .(station), transform, t1 =
 cumsum(DD(values_min,values_max,0,meanT)))
 
 I get a warning on my alpha calculation (NaN produced); indeed, the values I
 give as argument to asin() are out of the range [-1:1], as the print()
 reveals. I can't figure out how to solve this issue, because the same
 algorithm works in Excel (visual basic).

That doesn't mean that Excel and/or Visual Basic gives correct answers.

With the same input?
Then what does Excel say that asin(-7.4) evaluates to?
I tried asin(-1.2) and asin(-7.4)  in LibreOffice Calc (3.5.0) and got #VALUE! 
(Error: wrong data type) twice.

You'll have to present correct input to asin()  if you want to avoid the NaN's.

Berend

 It is very annoying, especially because it seems that no occurence of such
 error using that algorithm can be found on Internet.
 Any help is welcome :) Thanks for your time
 
 P.
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Issue-with-asin-tp4484462p4484462.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
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] Issue with asin()

2012-03-19 Thread Sarah Goslee
Hi,

You're not following the algorithm as given. The asin step shouldn't
be done for all values, but only for the ones that don't meet the
previous conditions. You're trying to calculate that step for ALL
values, then only use certain ones. You must instead subset the
values, THEN calculate that step. I would guess that your working
Excel version does follow the correct algorithm, but it's hard to know
for certain.

Here's a version that more closely follows the given reference:

MaxDailyTemp - values_max
MinDailyTemp - values_min
k - 0

GDD - rep(0, length(Tmin))
AvgDailyTemp - (MaxDailyTemp + MinDailyTemp)/2

# if MaxDailyTemp  k
# GDD = GDD + 0
# - add 0

# if MaxDailyTemp  k  MinDailyTemp  k
# GDD = GDD + AvgDailyTemp - k
GDD[MaxDailyTemp  k  MinDailyTemp  k] - AvgDailyTemp[MaxDailyTemp
 k  MinDailyTemp  k] - k

# if MaxDailyTemp  k  MinDailyTemp  k
# GDD = GDD + (1/pi) * [ (AvgDailyTemp – k) * ( ( pi/2 ) – arcsine(
theta ) ) + ( a * cos( arcsine( theta ) ) ) ]
a - (MaxDailyTemp - MinDailyTemp)/2
theta - ((k - AvgDailyTemp)/a)
GDD[MaxDailyTemp  k  MinDailyTemp  k] - (1/pi) * (
(AvgDailyTemp[MaxDailyTemp  k  MinDailyTemp  k] - k) * ( ( pi/2 ) -
asin( theta[MaxDailyTemp  k  MinDailyTemp  k] ) ) + (
a[MaxDailyTemp  k  MinDailyTemp  k] * cos( asin( theta[MaxDailyTemp
 k  MinDailyTemp  k] ) ) ) )

sum(GDD)

Sarah

On Mon, Mar 19, 2012 at 7:42 AM, Letnichev chatelain.p...@gmail.com wrote:
 Hello everyone,

 I am working for a few days already on a basic algorithm, very common in
 applied agronomy, that aims to determine the degree-days necessary for a
 given individual to reach a given growth stade. The algorithm (and context)
 is explained here:  http://www.oardc.ohio-state.edu/gdd/glossary.htm , and
 so I implemented my function in R as follows:

 DD - function(Tmin, Tmax, Tseuil, meanT, method = DDsin)
 ### function that calculates the degree-days based on
 ### minimum and maximum recorded temperatures and the
 ### minimal threshold temperature (lower growth temperature)
        {
        ### method arcsin
        if(method == DDsin){
                cond1 - (Tmax = Tseuil)
                cond2 - (Tmin = Tseuil)
                amp - ((Tmax - Tmin) / 2)
                print((Tseuil-meanT)/amp)
                alpha - asin((Tseuil - meanT) / amp)
                DD_ifelse3 - ((1 / pi) * ((meanT - Tseuil) * ((pi/2) - 
 alpha)) +
 amp*cos(alpha))

                DD - ifelse(cond1, 0, ifelse(cond2, (meanT - Tseuil), 
 DD_ifelse3))
        }

        ### method (Tmin + Tmax) / 2
        else if(method == DDt2){
                cond1 - (meanT  Tseuil)
                DD - ifelse(cond1,(meanT - Tseuil),0)
                }

        else{
                stop(\nMethod name is invalid.\nMethods available = DDsin 
 (sinus) or DDt2
 (mean)\n)
                }
        return(DD)
 }

 BUT! When I try to process random data:

 library(reshape2)
 library(plyr)

 station - rep(c(station1,station2,station3), 20)
 values_min - sample(-5:20, size = 60, replace = T)
 values_max - sample(20:40, size = 60, replace = T)
 meanT - ((values_min+values_max)/2)
 d - data.frame(station,values_min,values_max,meanT)
 names(d) - c(station, values_min,values_max,meanT)

 x-ddply(d, .(station), transform, t1 =
 cumsum(DD(values_min,values_max,0,meanT)))

 I get a warning on my alpha calculation (NaN produced); indeed, the values I
 give as argument to asin() are out of the range [-1:1], as the print()
 reveals. I can't figure out how to solve this issue, because the same
 algorithm works in Excel (visual basic).
 It is very annoying, especially because it seems that no occurence of such
 error using that algorithm can be found on Internet.
 Any help is welcome :) Thanks for your time

 P.


-- 
Sarah Goslee
http://www.functionaldiversity.org

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


Re: [R] Issue with asin()

2012-03-19 Thread Letnichev
Hello,

you're totally right, I tried first to control the flow with if
(MaxDailyTemp  k  MinDailyTemp  k){statement} but it was a bit messy.
Then ifelse() was supposed to help me out, but it didn't.

Thank you for your time, your code works exactly as I want :)

P.

--
View this message in context: 
http://r.789695.n4.nabble.com/Issue-with-asin-tp4484462p4485206.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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Issue with asin()

2012-03-19 Thread Letnichev
Yes, with the same input I had two different outputs with Excel and R. When
printing a debug report of Excel, it showed no anomalies and I am certain it
didn't calculate odd values (such as NaNs).
The way I coded was wrong, as Sarah said, I didn't follow completely the
algorithm. The solution she suggested works perfectly, so I am out of
trouble (for now :p ).

Thanks for your time,


P.

--
View this message in context: 
http://r.789695.n4.nabble.com/Issue-with-asin-tp4484462p4485185.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
and provide commented, minimal, self-contained, reproducible code.