Re: [Rd] Error message for infinite probability parameters in rbinom() and rmultinom()

[Christophe Dutang]

Thanks Duncan and Martin for your answer.

> Le 9 avr. 2023 à 00:33, Duncan Murdoch  a écrit :
> 
> On 08/04/2023 5:53 p.m., Martin Maechler wrote:
>>>>>>> Christophe Dutang
>>>>>>> on Sat, 8 Apr 2023 14:21:53 +0200 writes:
>> > Dear all,
>> > Using rmultinom() in a stochastic model, I found this function returns 
>> an error message 'NA in probability' for an infinite probability.
>> > Maybe, a more precise message will be helpful when debugging.
>> >> rmultinom(1, 3:5, c(1/2, 1/3, Inf))
>> > Error in rmultinom(1, 3:5, c(1/2, 1/3, Inf)) : NA in probability vector
>> >> rmultinom(1, 3:5, c(1/2, 1/3, NA))
>> > Error in rmultinom(1, 3:5, c(1/2, 1/3, NA)) : NA in probability vector
>> Thank you.
>> I agree the first ('Inf') should not do what it currently does,
>> and probably the 2nd one should neither give an error.
>> Note that in rmultinom,  the 'prob' is allowed to be *NOT*
>> scaled to sum(.) = 1.
>> Therefore 'Inf' makes sense as the limit (of a sequence) of (a)
>> very large number(s).
>> I claim that
>>   rmultinom(1, 3, c(1/2, 1/3, Inf))
>> should give the same as
>>   rmultinom(1, 3, c(1/2, 1/3, 1e300))
>> even without a warning,
> 
> That case makes sense, but is it worth the effort?  Certainly
> 
>rmultinom(1, 3, c(1/2, Inf, Inf))
> 
> can't give a useful answer because we don't know the relative size of the two 
> infinities.  


I think any infinite values in the probability vector or in the size vector 
should return NaN.


> I imagine the first NA comes from computing prob/sum(prob), which is c(0, 0, 
> NaN).

Line 73 of rmultinom.c check for probability in [0,1] and throw an error by 
calling C macro ML_WARN_ret_NAN() which return NA or -1.

I also notice that rbinom() throws a warning and not an error (for prob not in 
[0,1] or missing), e.g. 

> rbinom(4, 2:5, c(1/2, -1, Inf, NA))
[1]  2 NA NA NA

Maybe, only the man page should be updated in the Value field and state that NA 
is returned when parameter are not correct. By the way, in ?dbinom, it is said 
'If size is not an integer, NaN is returned.' which is not true for rbinom() 
and qbinom(), but only for dbinom() and pbinom().

Thanks in advance for your work

Kind regards, Christophe 

> 
> Duncan Murdoch
> 
> > and OTOH,  an NA in prob may return NA (and signal a warning)
> > instead of an error.
>> > For rgeom() or rbinom(), we got a warning for infinite probability :
>> Yes, but there, prob must be in [0,1] ... so that's somewhat differnt.
>> >> rbinom(1, 3, Inf)
>> > [1] NA
>> > Warning message:
>> > In rbinom(1, 3, Inf) : NAs produced
>> >> rbinom(1, 3, NA)
>> > [1] NA
>> > Warning message:
>> > In rbinom(1, 3, NA) : NAs produced
>> >> rgeom(1, Inf)
>> > [1] NA
>> > Warning message:
>> > In rgeom(1, Inf) : NAs produced
>> >> rgeom(1, NA)
>> > [1] NA
>> > Warning message:
>> > In rgeom(1, NA) : NAs produced
>> > Maybe, it could be better to harmonize the behavior for infinite 
>> probability.
>> > Kind regards, Christophe
>> >> sessionInfo()
>> > R version 4.2.3 (2023-03-15)
>> > Platform: aarch64-apple-darwin20 (64-bit)
>> > Running under: macOS Ventura 13.2.1
>> > Matrix products: default
>> > BLAS:   
>> /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
>> > LAPACK: 
>> /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
>> > locale:
>> > [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>> > attached base packages:
>> > [1] stats graphics  grDevices utils datasets  methods   base
>> > loaded via a namespace (and not attached):
>> > [1] compiler_4.2.3 tools_4.2.3
>> > -
>> > Christophe DUTANG
>> > LJK, Ensimag, Grenoble INP, UGA, France
>> > Web: http://dutangc.free.fr
>> > __
>> > R-devel@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-devel
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
> 

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Error message for infinite probability parameters in rbinom() and rmultinom()

[Christophe Dutang]

Dear all,

Using rmultinom() in a stochastic model, I found this function returns an error 
message 'NA in probability' for an infinite probability. 

Maybe, a more precise message will be helpful when debugging. 

>  rmultinom(1, 3:5, c(1/2, 1/3, Inf))
Error in rmultinom(1, 3:5, c(1/2, 1/3, Inf)) : NA in probability vector
>  rmultinom(1, 3:5, c(1/2, 1/3, NA))
Error in rmultinom(1, 3:5, c(1/2, 1/3, NA)) : NA in probability vector


For rgeom() or rbinom(), we got a warning for infinite probability :

> rbinom(1, 3, Inf)
[1] NA
Warning message:
In rbinom(1, 3, Inf) : NAs produced
> rbinom(1, 3, NA)
[1] NA
Warning message:
In rbinom(1, 3, NA) : NAs produced
> rgeom(1, Inf)
[1] NA
Warning message:
In rgeom(1, Inf) : NAs produced
> rgeom(1, NA)
[1] NA
Warning message:
In rgeom(1, NA) : NAs produced


Maybe, it could be better to harmonize the behavior for infinite probability.

Kind regards, Christophe


> sessionInfo()
R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1

Matrix products: default
BLAS:   
/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: 
/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base 

loaded via a namespace (and not attached):
[1] compiler_4.2.3 tools_4.2.3   

---------
Christophe DUTANG
LJK, Ensimag, Grenoble INP, UGA, France
Web: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z

[Christophe DUTANG]

You should take a look at https://CRAN.R-project.org/package=Rmpfr

Regards, Christophe 


Christophe Dutang
CEREMADE, Univ. Paris Dauphine, PSL Univ., France
Web : http://dutangc.free.fr


> Le 21 juin 2019 à 17:11, jing hua zhao  a écrit :
> 
> Dear Rui,
> 
> Thanks for your quick reply -- this allows me to see the bottom of this. I 
> was hoping we could have a handle of those p in genmoics such as 1e-300 or 
> smaller.
> 
> Best wishes,
> 
> 
> Jing Hua
> 
> 
> From: Rui Barradas 
> Sent: 21 June 2019 15:03
> To: jing hua zhao; r-devel@r-project.org
> Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
> 
> Hello,
> 
> Well, try it:
> 
> p <- .Machine$double.eps^seq(0.5, 1, by = 0.05)
> z <- qnorm(p/2)
> 
> pnorm(z)
> # [1] 7.450581e-09 1.22e-09 2.026908e-10 3.343152e-11 5.514145e-12
> # [6] 9.094947e-13 1.500107e-13 2.474254e-14 4.080996e-15 6.731134e-16
> #[11] 1.110223e-16
> p/2
> # [1] 7.450581e-09 1.22e-09 2.026908e-10 3.343152e-11 5.514145e-12
> # [6] 9.094947e-13 1.500107e-13 2.474254e-14 4.080996e-15 6.731134e-16
> #[11] 1.110223e-16
> 
> exp(z*z/2)
> # [1] 9.184907e+06 5.301421e+07 3.073154e+08 1.787931e+09 1.043417e+10
> # [6] 6.105491e+10 3.580873e+11 2.104460e+12 1.239008e+13 7.306423e+13
> #[11] 4.314798e+14
> 
> 
> p is the smallest possible such that 1 + p != 1 and I couldn't find
> anything to worry about.
> 
> 
> R version 3.6.0 (2019-04-26)
> Platform: x86_64-pc-linux-gnu (64-bit)
> Running under: Ubuntu 19.04
> 
> Matrix products: default
> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
> 
> locale:
>  [1] LC_CTYPE=pt_PT.UTF-8   LC_NUMERIC=C
>  [3] LC_TIME=pt_PT.UTF-8LC_COLLATE=pt_PT.UTF-8
>  [5] LC_MONETARY=pt_PT.UTF-8LC_MESSAGES=pt_PT.UTF-8
>  [7] LC_PAPER=pt_PT.UTF-8   LC_NAME=C
>  [9] LC_ADDRESS=C   LC_TELEPHONE=C
> [11] LC_MEASUREMENT=pt_PT.UTF-8 LC_IDENTIFICATION=C
> 
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods
> [7] base
> 
> other attached packages:
> 
> [many packages loaded]
> 
> 
> Hope this helps,
> 
> Rui Barradas
> 
> �s 15:24 de 21/06/19, jing hua zhao escreveu:
>> Dear R-developers,
>> 
>> I am keen to calculate exp(z*z/2) with z=qnorm(p/2) and p is very small. I 
>> wonder if anyone has experience with this?
>> 
>> Thanks very much in advance,
>> 
>> 
>> Jing Hua
>> 
>>   [[alternative HTML version deleted]]
>> 
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] 64-bit integer type warning on windows

[Christophe DUTANG]

Dear list,

During the last two months, I spent a lot of time trying to remove the 
following warnings of my package randtoolbox:

congruRand.c:180:3: warning: ISO C does not support the 'I64' ms_scanf length 
modifier [-Wformat=]
   sscanf(params[0], "%" PRIu64 "\n", &inp_mod);

Please see 
https://www.r-project.org/nosvn/R.check/r-devel-windows-ix86+x86_64/randtoolbox-00install.html
 

The warning is raised by the option -pedantic added when checking the package 
by R-devel. This warning is due to the use of the macro PRIu64. Note that this 
warning does not show up when I use PRIu64 with printf via Rprintf in other 
parts of the file congruRand.c. 

My first reaction to look at 
https://msdn.microsoft.com/fr-fr/library/f9ezhzhs(v=vs.110).aspx is useless, 
because R on windows uses Mingw-W64 integer types and not the MS version. 

Then, I tried to find information in such as Write Portable Code by Brian Hook 
and C11 (http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf).

They do recommend to use the macro SCNu64 (reserved for scanf family) rather 
than PRIu64 reserved for printf family functions. 

Unfortunately when I check the package with R-devel, I still have the warning.

On GitHub, there are a version of mingw64 sources : 
https://github.com/Alexpux/mingw-w64/ . It is explicitly written that SCNu64 is 
I64u, see line 210 of 
https://github.com/Alexpux/mingw-w64/blob/master/mingw-w64-headers/crt/inttypes.h
 

So, I was wondering if this warning was removable at all? Does anyone encounter 
this issue?

Any help is appreciated.

Thanks in advance, Christophe

PS: latest version of randtoolbox source code is here 
https://r-forge.r-project.org/scm/viewvc.php/pkg/randtoolbox/?root=rmetrics
--------
Christophe Dutang
CEREMADE, Univ. Paris Dauphine, France
Web : http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] registering Fortran routines in R packages

[Christophe Dutang]

Thanks Jari and Berend.

I was not aware of that function. I do use it : indeed it may extern 
declarations for Fortran routines. see 
https://r-forge.r-project.org/scm/viewvc.php/pkg/randtoolbox/src/init.c?view=markup&root=rmetrics
 
<https://r-forge.r-project.org/scm/viewvc.php/pkg/randtoolbox/src/init.c?view=markup&root=rmetrics>

However I get a new error when building 

Error: package or namespace load failed for ‘randtoolbox’ in dyn.load(file, 
DLLpath = DLLpath, ...):
 impossible de charger l'objet partagé 
'/Library/Frameworks/R.framework/Versions/3.4/Resources/library/randtoolbox/libs/randtoolbox.so':
  
dlopen(/Library/Frameworks/R.framework/Versions/3.4/Resources/library/randtoolbox/libs/randtoolbox.so,
 6): Symbol not found: _halton_f_
  Referenced from: 
/Library/Frameworks/R.framework/Versions/3.4/Resources/library/randtoolbox/libs/randtoolbox.so
  Expected in: flat namespace
 in 
/Library/Frameworks/R.framework/Versions/3.4/Resources/library/randtoolbox/libs/randtoolbox.so
Erreur : le chargement a échoué
Exécution arrêtée
ERROR: loading failed

Note that I change the name of Fortran routines to halton_f and sobol_f to 
avoid conflict with C version also renamed halton_c and sobol_c.

In the NAMESPACE, I put at the top useDynLib(randtoolbox, .registration = 
TRUE). 

Berend, I do look at your package… I still don’t figure out why it works for 
you!

Regards, Christophe

-----------
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr <http://dutangc.free.fr/>
> Le 10 mai 2017 à 08:56, Jari Oksanen  a écrit :
> 
> Have you tried using tools:::package_native_routine_registration_skeleton()? 
> If you don't like its output, you can easily edit its results and still avoid 
> most pitfalls.
> 
> Cheers, Jari Oksanen
> 
> From: R-devel  on behalf of Berend Hasselman 
> 
> Sent: 10 May 2017 09:48
> To: Christophe Dutang
> Cc: r-devel@r-project.org
> Subject: Re: [Rd] registering Fortran routines in R packages
> 
> Christophe,
> 
>> On 10 May 2017, at 08:08, Christophe Dutang  wrote:
>> 
>> Thanks for your email.
>> 
>> I try to change the name in lowercase but it conflicts with a C 
>> implementation also named halton. So I rename the C function halton2() and 
>> sobol2() while the Fortran function are HALTON() and SOBOL() (I also try 
>> lower case in the Fortran code). Unfortunately, it does not help since I get
>> 
>> init.c:97:25: error: use of undeclared identifier 'halton_'; did you mean 
>> 'halton2'?
>>  {"halton", (DL_FUNC) &F77_SUB(halton),  7},
>> 
>> My current solution is to comment FortEntries array and use 
>> R_useDynamicSymbols(dll, TRUE) for a dynamic search of Fortran routines.
> 
> Have a look at my package geigen and its init.c.
> Could it be that you are missing extern declarations for the Fortran routines?
> 
> 
> Berend
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel


[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] registering Fortran routines in R packages

[Christophe Dutang]

Thanks for your email.

I try to change the name in lowercase but it conflicts with a C implementation 
also named halton. So I rename the C function halton2() and sobol2() while the 
Fortran function are HALTON() and SOBOL() (I also try lower case in the Fortran 
code). Unfortunately, it does not help since I get

init.c:97:25: error: use of undeclared identifier 'halton_'; did you mean 
'halton2'?
  {"halton", (DL_FUNC) &F77_SUB(halton),  7},

My current solution is to comment FortEntries array and use 
R_useDynamicSymbols(dll, TRUE) for a dynamic search of Fortran routines.

Regards, Christophe
-------
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr <http://dutangc.free.fr/>
> Le 9 mai 2017 à 14:32, Berend Hasselman  a écrit :
> 
> 
>> On 9 May 2017, at 13:44, Christophe Dutang  wrote:
>> 
>> Dear list,
>> 
>> I’m trying to register Fortran routines in randtoolbox (in srt/init.c file), 
>> see 
>> https://r-forge.r-project.org/scm/viewvc.php/pkg/randtoolbox/src/init.c?view=markup&root=rmetrics.
>>  
>> 
>> Reading 
>> https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Registering-native-routines
>>  and looking at what is done in stats package, I first thought that the 
>> following code will do the job:
>> 
>> static const R_FortranMethodDef FortEntries[] = {
>> {"halton", (DL_FUNC) &F77_NAME(HALTON),  7},
>> {"sobol", (DL_FUNC) &F77_NAME(SOBOL),  11},
>> {NULL, NULL, 0}
>> };
>> 
>> But I got error messages when building : use of undeclared identifier 
>> ‘SOBOL_’. I also tried in lower case sobol and halton.
>> 
>> Looking at expm package 
>> https://r-forge.r-project.org/scm/viewvc.php/pkg/src/init.c?view=markup&revision=94&root=expm,
>>  I try  
>> 
>> static const R_FortranMethodDef FortEntries[] = {
>> {"halton", (DL_FUNC) &F77_SUB(HALTON),  7},
>> {"sobol", (DL_FUNC) &F77_SUB(SOBOL),  11},
>> {NULL, NULL, 0}
>> };
>> 
>> But the problem remains the same.
>> 
>> Is there a way to have header file for Fortran codes? how to declare 
>> routines defined in my Fortran file src/LowDiscrepancy.f?
>> 
> 
> lowercase routine names? manual does mention that.
> 
> Berend Hasselman
> 
> 
>> Any help appreciated
>> 
>> Regards, Christophe
>> ---
>> Christophe Dutang
>> LMM, UdM, Le Mans, France
>> web: http://dutangc.free.fr
>> 
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
> 


[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] registering Fortran routines in R packages

[Christophe Dutang]

Dear list,

I’m trying to register Fortran routines in randtoolbox (in srt/init.c file), 
see 
https://r-forge.r-project.org/scm/viewvc.php/pkg/randtoolbox/src/init.c?view=markup&root=rmetrics.
 

Reading 
https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Registering-native-routines
 and looking at what is done in stats package, I first thought that the 
following code will do the job:

static const R_FortranMethodDef FortEntries[] = {
 {"halton", (DL_FUNC) &F77_NAME(HALTON),  7},
 {"sobol", (DL_FUNC) &F77_NAME(SOBOL),  11},
 {NULL, NULL, 0}
};

But I got error messages when building : use of undeclared identifier ‘SOBOL_’. 
I also tried in lower case sobol and halton.

Looking at expm package 
https://r-forge.r-project.org/scm/viewvc.php/pkg/src/init.c?view=markup&revision=94&root=expm,
 I try  

static const R_FortranMethodDef FortEntries[] = {
 {"halton", (DL_FUNC) &F77_SUB(HALTON),  7},
 {"sobol", (DL_FUNC) &F77_SUB(SOBOL),  11},
 {NULL, NULL, 0}
};

But the problem remains the same.

Is there a way to have header file for Fortran codes? how to declare routines 
defined in my Fortran file src/LowDiscrepancy.f?

Any help appreciated

Regards, Christophe
---
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] names treatment in optim()

[Christophe Dutang]

Dear list,

A possible patch to correct the treatment of names consists in adding the 
following lines

if(!is.null(names(par))) names(res$par) <- names(par)  
if(!is.null(names(fn1(par names(res$value) <- names(fn1(par)) 

 just before returning the variable res in optim. That is 

optim <-
function(par, fn, gr = NULL, ...,
 method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", 
"Brent"),
 lower = -Inf, upper = Inf,
 control = list(), hessian = FALSE)
{
…

if (hessian)
res$hessian <- .External2(C_optimhess, res$par, fn1, gr1, con)
if(!is.null(names(par))) names(res$par) <- names(par)  
if(!is.null(names(fn1(par names(res$value) <- names(fn1(par)) 
res
}

Regards, Christophe

---
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr

> Le 17 sept. 2015 à 13:04, Christophe Dutang  a écrit :
> 
> Dear both,
> 
> I have found that names are not treated in the same way in optim() depending 
> on the optimization method (argument method). 
> 
> The example below shows the difference between the Brent method and the 
> L-BFGS-B method.
> 
> f <- function(x){ y <- x^2;names(y) <-"f(x)";y}
> optim(10, f, method="Brent", lower=-1, upper=10)$value
> optim(10, f, method="L-BFGS-B", lower=-1, upper=10)$value
> 
> z <- 10
> names(z) <- "x"
> z
> optim(z, f, method="Brent", lower=-1, upper=10)$par
> optim(z, f, method="L-BFGS-B", lower=-1, upper=10)$par
> 
> 
> Do you obtain the same behavior? would you consider interesting to have the 
> names passed to the components value and par?
> 
> I’m using R version 3.2.2 (2015-08-14)
> Platform: x86_64-apple-darwin13.4.0 (64-bit)
> Running under: OS X 10.10.5 (Yosemite)
> 
> locale:
> [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
> 
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base   
> 
> 
> 
> Kind regards, Christophe
> 
> ---
> Christophe Dutang
> LMM, UdM, Le Mans, France
> web: http://dutangc.free.fr
> 

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] names treatment in optim()

[Christophe Dutang]

Dear both,

I have found that names are not treated in the same way in optim() depending on 
the optimization method (argument method). 

The example below shows the difference between the Brent method and the 
L-BFGS-B method.

f <- function(x){ y <- x^2;names(y) <-"f(x)";y}
optim(10, f, method="Brent", lower=-1, upper=10)$value
optim(10, f, method="L-BFGS-B", lower=-1, upper=10)$value

z <- 10
names(z) <- "x"
z
optim(z, f, method="Brent", lower=-1, upper=10)$par
optim(z, f, method="L-BFGS-B", lower=-1, upper=10)$par


Do you obtain the same behavior? would you consider interesting to have the 
names passed to the components value and par?

I’m using R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base   



Kind regards, Christophe

---
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr <http://dutangc.free.fr/>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] names in optim()

[Christophe Dutang]

Dear both,

I have just found that names are treated in the same way in optim() depending 
on the optimization method. The example below shows the difference between the 
Brent method and the L-BFGS-B method.

f <- function(x){ y <- x^2;names(y) <-"f(x)";y}
optim(10, f, method="Brent", lower=-1, upper=10)$value
optim(10, f, method="L-BFGS-B", lower=-1, upper=10)$value

z <- 10
names(z) <- "x"
z
optim(z, f, method="Brent", lower=-1, upper=10)$par
optim(z, f, method="L-BFGS-B", lower=-1, upper=10)$par


Do you obtain the same behavior? 

I’m using R version 3.2.1 (2015-06-18)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)

locale:
[1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base 

loaded via a namespace (and not attached):
[1] tools_3.2.1 

Kind regards, CD
---
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Programming Tools CTV

[Christophe Dutang]

Dear Willem,

Personally, I use the R-forge project for the distribution CTV : 
https://r-forge.r-project.org/projects/ctv/

It’s an alternative option to github.

Regards, Christophe
---
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr

Le 23 janv. 2015 à 12:49, Luca Braglia  a écrit :

> Hi Willem
> 
> thanks for volounteering.
> 
> To the best of my knowledge (regarding the machinery side), if you're
> planning to use github (and maybe even if you don't) you can "stole"
> ideas from
> 
> https://github.com/ropensci/webservices
> https://github.com/lbraglia/PackageDevelopmentTaskView (minor
> modifications from webservices)
> https://github.com/eddelbuettel/ctv-finance or
> https://github.com/eddelbuettel/ctv-hpc
> 
> 
> HTH, Luca
> 
> 2015-01-23 11:13 GMT+01:00 Willem Ligtenberg
> :
>> Hi all,
>> 
>> Sorry if this doesn't end up in the thread.
>> Tobias Verbeke forwarded that e-mail to me, because he thought I would be 
>> interested in maintaining the Programming Tools CTV.
>> I wasn't subscribed to R-devel yet, but I would indeed like to volunteer to 
>> maintain the Programming Tools CTV.
>> 
>> It will be my first time creating a CTV, so some guidance on getting it 
>> setup will be appreciated.
>> I myself am very interested in better/easier ways to develop faster and 
>> nicer code.
>> 
>> Kind regards,
>> 
>> Willem
>> 
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [Wishlist] a 'PackageDevelopment' Task View

[Christophe Dutang]

Hi Luca,

Let me suggest to follow the table of contents of 
http://cran.r-project.org/doc/manuals/r-release/R-exts.html

For example, you could use the following TOC:

1 - Creating R pkg
2 - Writing documentation and vignettes
3 - Profiling code
4 - Debugging, spell checking
5 - Foreign languages interfaces
6 - GUI and other frond-ends 
7 - Unit testing
8 - Benchmarking
9 - Automation


Maybe 7 and 8 could be merged?

Regards, Christophe
—
Christophe Dutang
LMM, UdM, Le Mans, France
web: http://dutangc.free.fr

Le 26 juil. 2014 à 17:04, Luca Braglia  a écrit :

> 2014-07-25 14:29 GMT+02:00 Duncan Murdoch :
>> On 25/07/2014 8:05 AM, Luca Braglia wrote:
>>> 
>>> Hello everybody,
>>> 
>>> as a young/unexperienced R package developer (only a few, mainly
>>> for personal use) i was thinking it could be very useful having a
>>> "meta" task view for all package-development related
>>> packages and/or function.
>>> 
>>> Something like ...
>>> 
>>> Creation
>>> - utils::package.skeleton, pkgKitten, Rcpp::Rcpp.package.skeleton
>>> 
>>> Foreign languages interfaces:
>>> - Rcpp
>>> 
>>> Documentation
>>> - roxygen2
>>> 
>>> Profiling:
>>> - utils::Rprof
>>> - profr
>>> - proftools
>>> 
>>> Unit test
>>> - RUnit
>>> - testthat
>>> 
>>> Spell checking
>>> - tools::aspell_package_* functions
>>> 
>>> "Misc":
>>> - devtools
>>> 
>>> and so on.
>>> 
>>> These are only the ones i (use or know) & (remember), but for
>>> sure there is already a lot of useful code in this area and
>>> having a summary (by more experienced developers) of which good
>>> tools are available would be *very* useful, IMHO.
>>> 
>>> How about it?
>> 
>> 
>> Sounds like a good idea.  You should do it.
>> 
>> Download the "ctv" package, and read the vignette for instructions.
>> 
>> Duncan Murdoch
> 
> Fine, I've set up a github repo here
> 
> https://github.com/lbraglia/PackageDevelopmentTaskView
> 
> and a first thread about task view structure here
> 
> https://github.com/lbraglia/PackageDevelopmentTaskView/issues/1
> 
> Contribution (via e-mail or github) are, of course, really welcome.
> 
> Cheers, Luca
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Patch to add Beta binomial distribution. Mentor needed!

[Christophe Dutang]

Dear Joan,

Are you aware of this page 
http://cran.r-project.org/web/views/Distributions.html ?

If you want to contribute to R, you should write a package and submit it to 
CRAN. See http://cran.r-project.org/doc/manuals/R-exts.html

Regards

Christophe

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

Le 22 mai 2012 à 12:14, Joan Maspons a écrit :

> Hello,
> 
> I implemented the Beta binomial distribution following the patterns of the
> binomial distribution code and inspired by JAGS' code [1]. I have studied
> the code carefully but it's my first run in the R internals.
> 
> Can somebody review the code and if everything it's ok commit to the
> repository?
> 
> [1]
> http://mcmc-jags.hg.sourceforge.net/hgweb/mcmc-jags/mcmc-jags/file/15af65a4be29/src/modules/bugs/distributions/DBetaBin.cc
> 
> 
> -- 
> Joan Maspons
> CREAF (Centre de Recerca Ecològica i Aplicacions Forestals)
> Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Catalonia
> Tel +34 93 581 2915j.masp...@creaf.uab.cat
> http://www.creaf.uab.cat
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Beta binomial and Beta negative binomial

[Christophe Dutang]

Hi,

Please look at the distribution task view 
(http://cran.r-project.org/web/views/Distributions.html) and the package 
gamlss.dist.

By the way, distributions in R are implemented in /src/nmath directory 
and not /src/library/stats

Regards

Christophe

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

Le 16 mars 2012 à 18:41, Joan Maspons a écrit :

> Hi,
> I need Beta binomial and Beta negative binomial functions but in R there is
> only SuppDists package which provide this distributions using a limited
> parameter space of the generalized hypergeometric distribution (dghyper & Co.)
> which provide a limited parameter space for Beta binomial and Beta negative
> binomial functions (e.g. alpha + beta <1 in the Beta negative binomial).
> 
> I've done a checkout to R svn to seek the code for the implementation of
> distribution functions (dbinom et al.) but I haven't succeed. I don't know how
> difficult could be to implement Beta binomial and Beta negative binomial
> functions having the PDF  and CDF functions but I'm interested in implementing
> it. I have programming skills but I don't know much about R internals.
> 
> Where is the code? Can I implement these new functions inside stats
> package following the
> same patterns as other probability distributions?
> 
> Yours,
> -- 
> Joan Maspons
> CREAF (Centre de Recerca Ecològica i Aplicacions Forestals)
> Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Catalonia
> Tel +34 93 581 2915j.masp...@creaf.uab.cat
> http://www.creaf.uab.cat
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] First package submission to CRAN

[Christophe Dutang]

Hi,

By default, R CMD build makes sources, you have to use --binary if you want to 
get binaries. But you have to submit sources to the CRAN ftp server (and not 
binary). So just run a R CMD build.

C

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

Le 22 juin 2011 à 22:12, steven mosher a écrit :

> I'm preparing to submit my first package to CRAN, thanks to the help of too
> many people to mention.
> 
> I've built and checked the package on Windows  ( making a zip) and my path
> points to the 64 bit version of R.
> 
> Everything builds and checks and the final warnings have been fixed. My
> package is pure R with no source from
> 
> other languages.  My questions are  as follows. I've read the docs and just
> need a bit of clarification.
> 
> 1. For submission I should just build source  R CMD build mypkg  which
> outputs a tar.gz
> 2. Do I have to/ how do I build for 32 bit?
> 
> Thanks,
> 
> Steve
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Wish R Core had a standard format (or generic function) for "newdata" objects

[Christophe Dutang]

Among many solutions, I generally use the following code, which avoids the
ideal average individual, by considering the mean across of the predicted
values:

averagingpredict <- function(model, varname, varseq, type, subset=NULL)
{
if(is.null(subset))
mydata <- model$data
else
mydata <- model$data[subset, ]

f <- function(x)
{
mydata[, varname] <- x
mean(predict(model, newdata=mydata, type=type), na.rm=TRUE)
}

sapply(varseq, f)
}

It is time consuming, but it deals with non numeric variables.


Christophe


2011/4/26 Paul Johnson 

> Is anybody working on a way to standardize the creation of "newdata"
> objects for predict methods?
>
> When using predict, I find it difficult/tedious to create newdata data
> frames when there are many variables. It is necessary to set all
> variables at the mean/mode/median, and then for some variables of
> interest, one has to insert values for which predictions are desired.
> I was at a presentation by Scott Long last week and he was discussing
> the increasing emphasis in Stata on calculations of marginal
> predictions and "Spost" an several other packages, and,
> co-incidentally, I had a student visit who is learning to use R MASS's
> polr (W.Venables and B. Ripley) and we wrestled for quite a while to
> try to make the same calculations that Stata makes automatically.  It
> spits out predicted probabilities each independent variable, keeping
> other variables at a reference level.
>
> I've found R packages that aim to do essentially the same thing.
>
> In Frank Harrell's Design/rms framework, he uses a "data.dist"
> function that generates an object that the user has to put into the R
> options.  I think many users trip over the use of "options" there.  If
> I don't use that for a month or two, I completely forget the fine
> points and have to fight with it.  But it does "work" to give plots
> and predict functions the information they require.
>
> In  Zelig ( by Kosuke Imai, Gary King, and Olivia Lau), a function
> "setx" does the work of creating "newdata" objects. That appears to be
> about right as a candidate for a generic "newdata" function. Perhaps
> it could directly generalize to all R regression functions, but right
> now it is tailored to the models in Zelig. It has separate methods for
> the different types of models, and that is a bit confusing to me,since
> the "newdata" in one model should be the same as the newdata in
> another, I'm guessing. But his code is all there, I'll keep looking.
>
> In Effects (by John Fox), there are internal functions to create
> newdata and plot the marginal effects. If you load effects and run,
> for example, "effects:::effect.lm" you see Prof Fox has his own way of
> grabbing information from model columns and calculating predictions.
>
> I think it is time the R Core Team would look at this tell "us" what
> is the right way to do this. I think the interface to setx in Zelig is
> pretty easy to understand, at least for numeric variables.
>
> In R's termplot function, such a thing could be put to use.  As far as
> I can tell now, termplot is doing most of the work of creating a
> newdata object, but not exactly.
>
> It seems like it would be a shame to proliferate more functions that
> do the same function, when it is such a common thing.
>
> --
> Paul E. Johnson
> Professor, Political Science
> 1541 Lilac Lane, Room 504
> University of Kansas
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>



-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] on the Distribution task view

[Christophe Dutang]

Hi useRs,

As the maintainer of the Distribution task view (
http://cran.r-project.org/web/views/Distributions.html) for more than two
years, the following feedback exercise should have been done earlier. But
late is better than never!

I start this discussion to get your feedbacks/suggestions on the task view,
as well as to list the new packages that could be added to it.

Kind regards

Christophe

-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] An update of the Distribustions man page

[Christophe Dutang]

Dear list,

I would like to suggest a small update the ?Distributions man page of the
stats package. The current version contains the following line.

  The CRAN package \pkg{SuppDists} for additional distributions.

I think it would be better to put in this man page a link to the CRAN task
view on Distributions
http://cran.r-project.org/web/views/Distributions.html. It is not true
that the SuppDists package alone contains distributions
which are not implemented in R base.

Could it be possible to modify this sentence?

I don't know if there are other overview man pages in R related to other
task views? But a link to the task views on those man pages could increase
the popularity of the CRAN task views.


Kind regards

Christophe

-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] on the usage of do.call

[Christophe Dutang]

Thank you for help.

So I will look for another implementation not involving do_docall.

Actually, I want to pass a function and a list of arguments as argument. The
example lapply2 in section 5.11 is what I would like to do.

Christophe


2010/11/2 Duncan Murdoch 

> On 02/11/2010 11:28 AM, Christophe Dutang wrote:
>
>> Hello all,
>>
>> I don't know if it is possible, but I would like to use do.call in C code
>> in
>> my package. The function do.call is defined as
>> >  do.call
>> function (what, args, quote = FALSE, envir = parent.frame())
>> {
>> if (!is.list(args))
>> stop("second argument must be a list")
>> if (quote) {
>> enquote<- function(x) as.call(list(as.name("quote"),
>> x))
>> args<- lapply(args, enquote)
>> }
>> .Internal(do.call(what, args, envir))
>> }
>> 
>> >
>>
>> In/main/names.c, the function do.call is linked to the C function
>> do_docall (line 499). (.Internal calls the good function line 1194.)
>>
>> And the do_docall function is defined in/main/coerce.c line 2217 and
>> declared as
>>   SEXP attribute_hidden do_docall(SEXP call, SEXP op, SEXP args, SEXP rho)
>> in Internal.h.
>>
>> The problem is that I do not guess the exact meaning of theses arguments,
>> and the header file Internal.h is not found when included in my C file.
>>
>> As this header is not listed in R-exts.pdf section 6.17, I think I cannot
>> use the do_docall function.
>>
>> Does anyone face this problem? or have an idea on how to solve it?
>>
>
> You would rarely want to use do.call from C.  You should use eval().  See
> the examples in the Writing R Extensions manual, in the section "Evaluating
> R expressions from C".
>
> Duncan Murdoch
>



-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] on the usage of do.call

[Christophe Dutang]

Hello all,

I don't know if it is possible, but I would like to use do.call in C code in
my package. The function do.call is defined as
> do.call
function (what, args, quote = FALSE, envir = parent.frame())
{
if (!is.list(args))
stop("second argument must be a list")
if (quote) {
enquote <- function(x) as.call(list(as.name("quote"),
x))
args <- lapply(args, enquote)
}
.Internal(do.call(what, args, envir))
}

>

In /main/names.c, the function do.call is linked to the C function
do_docall (line 499). (.Internal calls the good function line 1194.)

And the do_docall function is defined in /main/coerce.c line 2217 and
declared as
  SEXP attribute_hidden do_docall(SEXP call, SEXP op, SEXP args, SEXP rho)
in Internal.h.

The problem is that I do not guess the exact meaning of theses arguments,
and the header file Internal.h is not found when included in my C file.

As this header is not listed in R-exts.pdf section 6.17, I think I cannot
use the do_docall function.

Does anyone face this problem? or have an idea on how to solve it?

Thanks in advance



Christophe



-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] on the optim function

[Christophe Dutang]

Thanks Ravi for your answer. You convince me!

Regards

Christophe

Le 6 août 2010 à 16:48, Ravi Varadhan a écrit :

> Christophe,
> 
> This is a good question. It is conventional that most classical optimization
> schemes (e.g. gradient-descent, quasi-Newton, and Newton-type) methods do
> not return the "iteration" counter, but they do return the number of times
> the gradient is computed.  These two are equivalent in the sense that the
> gradient is evaluated once and only once during each iteration.  However,
> the optimization codes always return the number of times the objective
> function is evaluated.  This is widely considered to be the most appropriate
> metric for measuring the computational effort of the algorithm.  The reason
> is this.  The actual step length of an algorithm for each iteration depends
> upon a "safe-guard" rule such as a line-search or a trust-region approach.
> These safe-guard rules might evaluate the function multiple times to ensure
> that the step length is "safe", for example, that it results in a reduction
> of the function value.  Therefore, just monitoring the number of iterations
> does not give you an accurate picture of the computational effort expended
> by the algorithm.
> 
> Another point worth noting is that when an analytic gradient function is not
> specified by the user, most (if not all) gradient-based algorithms use a
> numerical approximation for the gradient.  This is typically done using a
> first-order forward difference scheme.  This involves `p' evaluations of the
> objective function, when you have a p-dimensional parameter vector to be
> optimized.  This greatly increases the computational effort, much more than
> that involved in safe-guarding.  Therefore, it is always a good idea to
> specify analytic gradient, whenever possible.
> 
> Here is an example (based on that from optim's help page) that illustrates
> the equivalency of # iterations and gradient count:
> 
>> optim(c(-1.2,1), fr, grr, method="CG", control=list(maxit=50))
> $par
> [1] -0.9083598  0.8345932
> 
> $value
> [1] 3.636803
> 
> $counts
> function gradient 
> 202   51 
> 
> $convergence
> [1] 1
> 
> $message
> NULL
> 
> 
>> optim(c(-1.2,1), fr, method="CG", control=list(maxit=75))$counts
> function gradient 
> 304   76 
>> 
> Note that when I specified `maxit = 50', I got gradient counts = 51, and
> when I specified maxit=75, I got gradient counts = 76.
> 
> But look at the difference in function counts.  There is 102 more function
> evals, even though I did only 25 more iterations.  This tells you that the
> function is evaluated more than once during an iteration.
> 
> The above is also true for BFGS, L-BFGS-B.  You can check it for yourself.
> 
> 
> Hope this helps,
> Ravi.
> 
> 
> -Original Message-
> From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-project.org]
> On Behalf Of Christophe Dutang
> Sent: Friday, August 06, 2010 5:15 AM
> To: r-devel@r-project.org
> Subject: [Rd] on the optim function
> 
> Dear useRs,
> 
> I have just discovered that the R optim function does not return the number
> of iterations. 
> 
> I still wonder why line 632-634 of optim C, the iter variable is not
> returned (for the BFGS method for example) ? 
> 
> Is there any trick to compute the iteration number with function call
> number?
> 
> Kind regards
> 
> Christophe
> 
> --
> Christophe Dutang
> Ph.D. student at ISFA, Lyon, France
> website: http://dutangc.free.fr
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
> 

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] on the optim function

[Christophe Dutang]

Dear useRs,

I have just discovered that the R optim function does not return the number of 
iterations. 

I still wonder why line 632-634 of optim C, the iter variable is not returned 
(for the BFGS method for example) ? 

Is there any trick to compute the iteration number with function call number?

Kind regards

Christophe

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] constrained optimization

[Christophe Dutang]

Dear Paul,

No I haven't tried to use L-BFGS-B algorithm with bounds. Let me check if a 
variable change can be used on my problem.

Thanks

Christophe

Le 7 juil. 2010 à 19:49, Paul Bailey a écrit :

> 
> Have you considered using the optim function with L-BFGS-B for bounded
> optimization. Obviously, you will have to do changes of variables so that
> everything is in terms of a rectangle (which is the type of bounding that it
> accepts).
> 
> I believe it is based on "A LIMITED MEMORY ALGORITHM FOR BOUND CONSTRAINED
> OPTIMIZATION" By Richard H Byrd,  Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.
> Technical Report NAM-08, May 1994.
> 
> They also have a paper on the exact implementation. I also believe that R
> uses the code from Nocedal et al.
> -- 
> View this message in context: 
> http://r.789695.n4.nabble.com/constrained-optimization-tp2280809p2281282.html
> Sent from the R devel mailing list archive at Nabble.com.
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] constrained optimization

[Christophe Dutang]

Hello Ravi,

Yes I'm willing to test it. Tell me when there is an R version of it.

Christophe

2010/7/7 Ravi Varadhan 

> Paul,
>
> SQUAREM is not an optimization package.  It is an "accelerator" of
> fixed-point iteration such as the EM, MM, and other algorithms.
>
> I am referring to a set of functions that I have written, which I call
> "ALABAMA" for Augmented Lagrangian Adaptive Barrier Minimization Algorithm.
>  I think it should be package up pretty soon as there seems to be a need for
> nonlinearly constrained optimization in R.
>
> Best,
> Ravi.
>
> -Original Message-
> From: Paul Gilbert [mailto:pgilb...@bank-banque-canada.ca]
> Sent: Wednesday, July 07, 2010 9:18 AM
> To: Ravi Varadhan; Christophe Dutang
> Subject: RE: [Rd] constrained optimization
>
> Ravi
>
> I think you are referring to the SQUAREM package.  The development
> version is available at https://r-forge.r-project.org/R/?group_id=395
> and can be installed with install.packages("SQUAREM",
> repos="http://R-Forge.R-project.org";) .
>
> Paul
>
> >-Original Message-
> >From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-
> >project.org] On Behalf Of Ravi Varadhan
> >Sent: July 7, 2010 9:10 AM
> >To: 'Christophe Dutang'; r-devel@r-project.org
> >Subject: Re: [Rd] constrained optimization
> >
> >Hi Christophe,
> >
> >I have an algorithm for solving nonlinearly constrained optimization.
> It
> >is
> >a combination of an interior point (for inequalities) algorithm with an
> >augmented Lagrangian (for equalities).  It is coded entirely in R, and
> >hence
> >is a bit slow, but it seems to do the job quite robustly in terms of
> >handling poor starting values.  I can send this to you, if you are
> >interested.
> >
> >Ravi.
> >
> >-Original Message-
> >From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-
> >project.org]
> >On Behalf Of Christophe Dutang
> >Sent: Wednesday, July 07, 2010 8:01 AM
> >To: r-devel@r-project.org
> >Subject: [Rd] constrained optimization
> >
> >Dear list,
> >
> >The task view on optimization does not reference a package for non
> >linear
> >constrained optimization problems. Stefan Theussl told me to look at
> the
> >Rsolnp package, but unfortunately it is not very clear what method is R
> >ported. (The authors ported the matlab code of Yinyu Ye
> >http://www.stanford.edu/~yyye/ <http://www.stanford.edu/%7Eyyye/> <
> http://www.stanford.edu/%7Eyyye/>)
> >
> >Currently I'm looking for an implementation of sequential quadratic
> >programming to replicate SNOPT*. A good reference I found on the web is
> >this
> >booklet
> >http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/5456/pdf/imm5456.p
> d
> >f .
> >
> >Does anyone know an implementation of such algorithms? Is there any
> >fortran
> >implementation available useful if I have to implement it?
> >
> >Thanks in advance
> >
> >Christophe
> >
> >
> >* SNOPT: An SQP Algorithm For Large-Scale Constrained Optimization
> >(1997) by
> >Philip E. Gill ,  Walter Murray ,  Michael ,  Michael A. Saunders
> >--
> >Christophe DUTANG
> >Ph. D. student at ISFA, Lyon, France
> >
> >   [[alternative HTML version deleted]]
> >
> >__
> >R-devel@r-project.org mailing list
> >https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> >__
> >R-devel@r-project.org mailing list
> >https://stat.ethz.ch/mailman/listinfo/r-devel
>
> 
>
> La version française suit le texte anglais.
>
>
> 
>
> This email may contain privileged and/or confidential information, and the
> Bank of
> Canada does not waive any related rights. Any distribution, use, or copying
> of this
> email or the information it contains by other than the intended recipient
> is
> unauthorized. If you received this email in error please delete it
> immediately from
> your system and notify the sender promptly by email that you have done so.
>
>
> --------
>
> Le présent courriel peut contenir de l'information privilégiée ou
> confidentielle.
> La Banque du Canada ne renonce pas aux droits qui s'y rapportent. Toute
> diffusion,
> utilisation ou copie de ce courriel ou des renseignements qu'il contient
> par une
> personne autre que le ou les destinataires désignés est interdite. Si vous
> recevez
> ce courriel par erreur, veuillez le supprimer immédiatement et envoyer sans
> délai à
> l'expéditeur un message électronique pour l'aviser que vous avez éliminé de
> votre
> ordinateur toute copie du courriel reçu.
>
>


-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] constrained optimization

[Christophe Dutang]

Thank you Jelmer,

I have just spent five minutes on the website, but I'm not sure to
understand which method is implemented. Is it *this method?
On the Implementation of a Primal-Dual Interior Point Filter Line Search
Algorithm for Large-Scale Nonlinear Programming*

Do you plan to put lpopt on CRAN? I'm willing to test it!

Christophe

2010/7/7 Jelmer Ypma 

> Hi Christophe,
>
> I wrote an R interface to Ipopt ( https://projects.coin-or.org/Ipopt
> ), which does non-linear constrained optimization for very large
> problems. Would this be of interest to you? It took some time to find
> out which license I could use to distribute it, but I'm planning to
> make it publicly available this weekend.
>
> Jelmer
>
> On Wed, Jul 7, 2010 at 13:01, Christophe Dutang  wrote:
> > Dear list,
> >
> > The task view on optimization does not reference a package for non linear
> > constrained optimization problems. Stefan Theussl told me to look at the
> > Rsolnp package, but unfortunately it is not very clear what method is R
> > ported. (The authors ported the matlab code of Yinyu Ye
> > http://www.stanford.edu/~yyye/ <http://www.stanford.edu/%7Eyyye/> <
> http://www.stanford.edu/%7Eyyye/>)
> >
> > Currently I'm looking for an implementation of sequential quadratic
> > programming to replicate SNOPT*. A good reference I found on the web is
> this
> > booklet
> >
> http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/5456/pdf/imm5456.pdf.
> >
> > Does anyone know an implementation of such algorithms? Is there any
> fortran
> > implementation available useful if I have to implement it?
> >
> > Thanks in advance
> >
> > Christophe
> >
> >
> > * SNOPT: An SQP Algorithm For Large-Scale Constrained Optimization (1997)
> by
> > Philip E. Gill ,  Walter Murray ,  Michael ,  Michael A. Saunders
> > --
> > Christophe DUTANG
> > Ph. D. student at ISFA, Lyon, France
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-devel@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
> >
>



-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] constrained optimization

[Christophe Dutang]

Dear list,

The task view on optimization does not reference a package for non linear
constrained optimization problems. Stefan Theussl told me to look at the
Rsolnp package, but unfortunately it is not very clear what method is R
ported. (The authors ported the matlab code of Yinyu Ye
http://www.stanford.edu/~yyye/ <http://www.stanford.edu/%7Eyyye/>)

Currently I'm looking for an implementation of sequential quadratic
programming to replicate SNOPT*. A good reference I found on the web is this
booklet
http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/5456/pdf/imm5456.pdf .

Does anyone know an implementation of such algorithms? Is there any fortran
implementation available useful if I have to implement it?

Thanks in advance

Christophe


* SNOPT: An SQP Algorithm For Large-Scale Constrained Optimization (1997) by
Philip E. Gill ,  Walter Murray ,  Michael ,  Michael A. Saunders
-- 
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] discrete-time survival model

[Christophe Dutang]

Dear List,

Last week, I asked on R-help if there was a package implementing the dynamic 
discrete time model of Fahrmeir (1994) - http://www.jstor.org/pss/2336962. 
Unfortunately, nobody answered. So I start to implement it. 

Does anyone interested in such development? Or is it useless because someone 
has already implemented this model?

In the coming days, I plan to move my code on R-forge, so people can easily 
contribute. If you are interested please let me know.

Thanks in advance

Christophe


--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr







[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] p-generalized normal distribution

[Christophe Dutang]

He already asked... And I don t know any pkg implementing this  
distribution...


iPhone.fan

Le 26 nov. 2009 à 21:22, Kjetil Halvorsen m> a écrit :


There is a CRAN Task View for probability Distributions. take a look  
there!


Kjetil

On Tue, Nov 24, 2009 at 1:53 PM, Steve Kalke > wrote:

Hello,

I would like to know if there is an R-package available for  
computing the

density, distribution function, quantiles and random numbers of the
p-generalized normal distribution or if somebody is already working  
on it.


Best regards,
Steve Kalke

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel



__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Fwd: AR(2) modelling

[Christophe Dutang]

Thanks Gabor,

My mistake was I used the 'cor' function rather than the function  
'acf' with "correlation" argument.


Christophe

Le 15 nov. 2009 à 18:13, Gabor Grothendieck a écrit :


Try this:

debug(stats:::ar.yw.default)

and then run the ar function to step through the ar code so you can
see the results line by line and understand what it is doing at a very
detailed level.

On Sun, Nov 15, 2009 at 10:06 AM, Christophe Dutang  
 wrote:
As you are sure of the accuracy of your code, why don't you tell me  
where is

my mistake?

Le 15 nov. 2009 à 12:03, Prof Brian Ripley a écrit :




--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595


--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel



--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Fwd: AR(2) modelling

[Christophe Dutang]

As you are sure of the accuracy of your code, why don't you tell me  
where is my mistake?


Le 15 nov. 2009 à 12:03, Prof Brian Ripley a écrit :




--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595


--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Fwd: AR(2) modelling

[Christophe Dutang]

My email does not seem to receive any attention on R-help, so I  
forward it on R-devel if someone has already faced the "problem".

Thanks in advance

Christophe

Début du message réexpédié :

> De : Christophe Dutang 
> Date : 13 novembre 2009 23:30:14 HNEC
> À : r-h...@r-project.org
> Objet : AR(2) modelling
>

Hi useRs,

I'm trying to fit a basic AR(2) model with the 'ar' function. And when  
I try to check the value of the coefficients, I could not find the  
same value as the 'ar' function.

Here is my example:
myserie <- c(212, 205, 210, 213, 217, 222, 216, 218, 220, 212, 215, 236)

#plot(myserie, type="l")

myserieminus0 <- tail(myserie, -2)
myserieminus1 <- tail(head(myserie, -1), -1)
myserieminus2 <- head(myserie, -2)


###Yule Walker equations

r1 <- cor(myserieminus0, myserieminus1)
r2 <- cor(myserieminus0, myserieminus2)

#method 1
phihat1 <- r1*(1-r2)/(1-r1^2)
phihat2 <- (r2-r1^2)/(1-r1^2)

#method 2
bigR <- cbind(c(1, r1), c(r1, 1))
smallr <- c(r1, r2)
ressolve <- solve(bigR, smallr)

resaryw <- ar(myserie, method="yule-walker", order.max=2, aic=FALSE)

cat("\t\tmanual YW 1\t\tmanual YW 2\t\tar YW\n")
cat("first coef", phihat1,"\t", ressolve[1],"\t\t", resaryw$ar[1], "\n")
cat("first coef", phihat2,"\t", ressolve[2],"\t", resaryw$ar[2], "\n\n")

 >  manual YW 1 manual YW 2 ar YW
 > first coef 0.2941808  0.2941808   0.1869641
 > first coef -0.1316839 -0.1316839  -0.1038001

I do not understand why the "yule-walker" does not solve exactly the  
Yule-Walker equations. A reference can be found here 
http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YWSourceFiles/YW-Eshel.pdf
 
  .

In the R source code (/src/library/stats/R), the file ar.R  
contains the ar.yw.default function implements the function. For  
univariate case (line 130), r_eureka function is used, which seems to  
be implemented in the eureka.f function.

  subroutine eureka (lr,r,g,f,var,a)
c
c  solves Toeplitz matrix equation toep(r)f=g(1+.)
c  by Levinson's algorithm
c  a is a workspace of size lr, the number
c  of equations
c

is supposed to implement the Yule-Walker equations...

Any help is welcome.

Just to be sure, I can do something I try to reconcile the ordinary  
least square methods. And it works!

Christophe


PS : OLS code

###Ordinary Least Square

reslm <- lm(myserieminus0 ~ myserieminus1 + myserieminus2)
#summary(reslm)
coef1ols <- reslm$coefficients["myserieminus1"]
coef2ols <- reslm$coefficients["myserieminus2"]

resarols <- ar(myserie, method="ols", order.max=2, aic=FALSE)
cat("\t\tmanual ols\t\tar ols\n")
cat("first coef", coef1ols,"\t", resarols$ar[1], "\n")
cat("first coef", coef2ols,"\t", resarols$ar[2], "\n\n")


--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr







[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [R-SIG-Mac] rnorm.halton

[Christophe Dutang]

Thanks Simon for your answers,

I understand only pointer size change between 32bit and 64bit in  
Fortran, however my problem is with this function:


C  
--

  REAL*8 FUNCTION HQNORM(P)

C USED TO CALCULATE HALTON NORMAL DEVIATES:
  REAL*8 P,R,T,A,B, EPS
  DATA P0,P1,P2,P3,P4, Q0,Q1,Q2,Q3,Q4
 &   /-0.322232431088E+0, -1.E+0, -0.342242088547E+0,
 &-0.204231210245E-1, -0.453642210148E-4, +0.993484626060E-1,
 &+0.588581570495E+0, +0.531103462366E+0, +0.103537752850E+0,
 &+0.385607006340E-2  /

C NOTE, IF P BECOMES 1, THE PROGRAM FAILS TO CALCULATE THE
C NORMAL RDV. IN THIS CASE WE REPLACE THE LOW DISCREPANCY
C POINT WITH A POINT FAR IN THE TAILS.
  EPS = 1.0D-6
  IF (P.GE.(1.0D0-EPS)) P = 1.0d0 - EPS
  IF (P.LE.EPS) P = EPS
  IF (P.NE.0.5D0) GOTO 150
  HQNORM = 0.0D0
  RETURN
150   R = P
  IF (P.GT.0.5D0) R = 1.0 - R
  T = DSQRT(-2.0*DLOG(R))
  A = T*P4 + P3)*T+P2)*T + P1)*T + P0)
  B = T*Q4 + Q3)*T+Q2)*T + Q1)*T + Q0)
  HQNORM = T + (A/B)
  IF (P.LT.0.5D0) HQNORM = -HQNORM

  RETURN
  END
C  
--


Despite this approximation of the normal quantile function is very  
basic, I think there is no pointer used? Can someone tell me if I'm  
wrong? So I do not understand why rnorm.halton does not work on 64bit  
machine... runif.halton works! (so the rest of the code is ok!)


I'm little bit curious why it does not work... but I think I will just  
remove the use of this subroutine and use directly the R function qnorm.


Christophe

Le 10 oct. 2009 à 20:02, Simon Urbanek a écrit :


Christophe,

I forgot to answer the second part of your e-mail -- see below.

On Oct 10, 2009, at 12:04 PM, Christophe Dutang wrote:


Hi all,

I need to transform classic 32bit Fortran code to 64bit Fortran  
code, see the discussion [R-SIG-Mac] rnorm.halton. But I'm clearly  
a beginner in Fortran...


Does someone already do this for his package?

From here, http://techpubs.sgi.com/library/tpl/cgi-bin/getdoc.cgi?coll=linux&db=bks&fname=/SGI_Developer/Porting_Guide/ch03.html 
 , I identify the following changes (Fortran types) to the Fortran  
code:


- INTEGER to INTEGER*8
- REAL*8 to REAL*16

The code I would like to change is available on R forge here : 
http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/randtoolbox/src/LowDiscrepancy.f?rev=4229&root=rmetrics&view=markup

Another question I have is how do I tell the R code to use the  
64bit version of my code on 64bit machine?


In the current implementation, the file quasiRNG.R calls directly  
the Fortran code with .Fortran.

How could I use the 64bit version directly in the R code?



You don't have a choice -- whether 32 or 64 bits are used is defined  
by R, because you cannot mix 32-bit and 64-bit code. So in 64-bit R  
your package will be compiled into 64-bit. In 32-bit R your package  
will be compiled in 32-bit. You can see which R you are using with


8 * .Machine$sizeof.pointer

(Note: The Mac binaries we provide are multi-arch so they [the  
Leopard build] do in fact include both 32-bit and 64-bit -- which  
one gets started depends on the --arch flag, see R docs for details  
on multi-arch installations).


Cheers,
Simon


I suspect I need to have a quasiRNG.c file where I will use  
preprocessor statement that will select the good version of the  
function to call. Is that correct?




No (see previous e-mail).



Thanks in advance

Christophe


Le 15 sept. 2009 à 18:25, Anirban Mukherjee a écrit :


Sorry, what I should have said was Halton numbers are quasi-random,
and not pseudo-random. Quasi-random is the technically appropriate
terminology.

Halton sequences are low discrepancy: the subsequence looks/smells
random. Hence, they are often used in quasi monte carlo simulations.
To be precise, there is only 1 Halton sequence for a particular  
prime.

Repeated calls to Halton should return the same numbers. The first
column is the Halton sequence for 2. the second for 3, the third  
for 5

and so on using the first M primes (for M columns). (You can also
scramble the sequence to avoid this.)

I am using them to integrate over a multivariate normal space. If  
you
take 1000 random draws, then sum f() over the draws is the  
expectation

of f(). If f() is very non-linear (and/or multi-variate) then even
with large N, its often hard to get a good integral. With quasi- 
random
draws, the integration is better for the same N. (One uses the  
inverse
distribution function.) For an example, you can look at Train's  
paper

(page 4 and 5 have a good explanation) at:

http://elsa.berkeley.edu/wp/train0899.pdf

In the context of simulated maximum likelihood estimation, such
integrals are very common. Of course true randomness

Re: [Rd] [R-SIG-Mac] rnorm.halton

[Christophe Dutang]

irban Mukherjee | Assistant Professor, Marketing | LKCSB, SMU
5062 School of Business, 50 Stamford Road, Singapore 178899 |  
+65-6828-1932


___
R-SIG-Mac mailing list
r-sig-...@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-mac


--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] rowSums, rowMean and rowCumSums?

[Christophe Dutang]

Dear list,

Don't you think it could be useful to have in R base a function  
rowCumSums, that compute cumulative sums for each row of a matrix?

My implementation of rowCumSums is

rowCumSums <- function(x) t(mapply(function(row)cumsum(x[row,]),  
1:NROW(x)))

I'm sure it can be improved to have other arguments like na.rm or dims.

Is there any hope to have this function in R?

Christophe

--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr







[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] user supplied random number generators

[Christophe Dutang]

Hello,

Le 30 juil. 09 à 08:21, Ross Boylan a écrit :


?Random.user says (in svn trunk)
 Optionally,
 functions \code{user_unif_nseed} and \code{user_unif_seedloc} can be
 supplied which are called with no arguments and should return  
pointers

 to the number of seeds and to an integer array of seeds.  Calls to
 \code{GetRNGstate} and \code{PutRNGstate} will then copy this array  
to

 and from \code{.Random.seed}.
And it offers as an example
 void  user_unif_init(Int32 seed_in) { seed = seed_in; }
 int * user_unif_nseed() { return &nseed; }
 int * user_unif_seedloc() { return (int *) &seed; }

First question: what is the lifetime of the buffers pointed to by the
user_unif-* functions, and who is responsible for cleaning them up?   
In

the help file they are static variables, but in general they might be
allocated on the heap or might be in structures that only persist as
long as the generator does.

Since the example uses static variables, it seems reasonable to  
conclude

the core R code is not going to try to free them.

Second, are the types really correct?  The documentation seems quite
explicit, all the more so because it uses Int32 in places.  However,  
the

code in RNG.c (RNG_Init) says

ns = *((int *) User_unif_nseed());
if (ns < 0 || ns > 625) {
warning(_("seed length must be in 0...625; ignored"));
break;
}
RNG_Table[kind].n_seed = ns;
RNG_Table[kind].i_seed = (Int32 *) User_unif_seedloc();
consistent with the earlier definition of RNG_Table entries as
typedef struct {
   RNGtype kind;
   N01type Nkind;
   char *name; /* print name */
   int n_seed; /* length of seed vector */
   Int32 *i_seed;
} RNGTAB;

This suggests that the type of user_unif_seedloc is Int32*, not int *.
It also suggests that user_unif_nseed should return the number of 32  
bit

integers.  The code for PutRNGstate(), for example, uses them in just
that way.

While the dominant model, even on 64 bit hardware, is probably to  
leave

int as 32 bit, it doesn't seem wise to assume that is always the case.


You can test the size of an int with a configure script. see for  
example the package foreign, the package randtoolbox (can be found in  
Rmetrics R forge project) I maintain with Petr Savicky.


By the way, I'm sure he has an answer about RNGkind because he made  
the runif interface in the randtoolbox package and in rngWELL19937  
package.


Christophe




I got into this because I'm trying to extend the rsprng code; sprng
returns its state as a vector of bytes.  Converting these to a  
vector of

integers depends on the integer length, hence my interest in the exact
definiton of integer.  I'm interested in lifetime because I believe
those bytes are associated with the stream and become invalid when the
stream is freed; furthermore, I probably need to copy them into a  
buffer

that is padded to full wordlength.  This means I allocate the buffer
whose address is returned to the core R RNG machinery.  Eventually
somebody needs to free the memory.

Far more of my rsprng adventures are on
http://wiki.r-project.org/rwiki/doku.php?id=packages:cran:rsprng.   
Feel

free to read, correct, or extend it.

Thanks.

Ross Boylan

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] cran task view

[Christophe Dutang]

Ok,

For the moment, there is only this in the CTV on finance
"
several packages provide functionality for Extreme Value Theory  
models: evd, evdbayes, evir, extRremes, ismev, POT.

"

Before the useR conference, I tried to think on a future section on  
EVT to start a discussion.


Thanks for your quick answer

Christophe

Le 11 juin 09 à 14:34, Dirk Eddelbuettel a écrit :



On 11 June 2009 at 14:11, Achim Zeileis wrote:
| > Yesterday, I try to look at a CRAN task view on extreme value  
theory (EVT).

| > Does anyone plan to do a CRAN task view on this topic?
|
| Not that I know of.
|
| > Packages related to this topic seem to be spread among finance,  
econometrics

| > and social econometrics.
|
| Maybe more on the finance side...

I also think that we shouldn't get too granular. 20-some task views  
is fine,

200-some less so. But that is Achim's editorial call to make.

| > Maybe this discussion could wait this july useR?
|
| Sure, we can talk, probably together with Dirk. I guess that the  
number
| of packages for EVT is also not excessive so that we might be able  
to

| include it in the Finance view.

Agreed. Patches to the task view are always welcome.  :-)

Dirk

--
Three out of two people have difficulties with fractions.


--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] cran task view

[Christophe Dutang]

Hi all,

Yesterday, I try to look at a CRAN task view on extreme value theory  
(EVT). Does anyone plan to do a CRAN task view on this topic?


Packages related to this topic seem to be spread among finance,  
econometrics and social econometrics.


Maybe this discussion could wait this july useR?

Thanks in advance for any advice

Christophe

--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] statutes of R Foundation

[Christophe Dutang]

Good evening all,

I realised yesterday that the R Foundation statutes doc was only  
available in English and in German.


I tried to translate into French: a first version is available here http://dutangc.free.fr/pub/statut%20R.pdf 
 .


Could you please tell me what do you think of my translation?

Thanks in advance



Christophe

PS : Tex sources are here http://dutangc.free.fr/pub/statut%20R.tex

--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Fwd: Vignettes with missing or empty \VignetteIndexEntry:

[Christophe Dutang]

I also have another warning while checking the package

* checking package vignettes in 'inst/doc' ... WARNING
Package vignettes without corresponding PDF:
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
chi.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
contextra.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
discrete.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
discrextra.Rnw
   .../These/code/R/probdistr/Rforge/distributions/pkg/inst/doc/ 
probdistr-exp.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
finitesupp.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
gaussian.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
general.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
gumbel.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
intro.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
logis.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
math.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
misc.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
mult.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
pareto.Rnw
   .../code/R/probdistr/Rforge/distributions/pkg/inst/doc/probdistr- 
student.Rnw
* checking PDF version of manual ... OK



Début du message réexpédié :

> De : Christophe Dutang 
> Date : 10 mai 2009 11:52:17 HAEC
> À : r-devel@r-project.org
> Objet : Vignettes with missing or empty \VignetteIndexEntry:
>
> Hi,
>
> I have a problem when checking the package 'probdistr' (on  
> probability distributions).
> I got this warning
>
> * checking index information ... WARNING
> Vignettes with missing or empty \VignetteIndexEntry:
> [1] "probdistr-chi""probdistr-contextra"  "probdistr-discrete"
> [4] "probdistr-discrextra" "probdistr-exp""probdistr- 
> finitesupp"
> [7] "probdistr-gaussian"   "probdistr-general""probdistr-gumbel"
> [10] "probdistr-intro"  "probdistr-logis"  "probdistr-math"
> [13] "probdistr-misc"   "probdistr-mult"   "probdistr-pareto"
> [16] "probdistr-student"
> See the information on INDEX files and package subdirectories in the
> chapter 'Creating R packages' of the 'Writing R Extensions' manual.
>
> The warning tells me \VignetteIndexEntry is missing however there is  
> one in probdistr-main.Rnw (Sweave file calling all other files via  
> \SweaveInput{probdistr-X}). Let us note this file is not listed  
> in the above files...
>
> But the build is ok
>
> * checking for file 'pkg/DESCRIPTION' ... OK
> * preparing 'pkg':
> * checking DESCRIPTION meta-information ... OK
> * installing the package to re-build vignettes
> * Installing *source* package ‘probdistr’ ...
> ** R
> ** inst
> ** help
> *** installing help indices
> >>> Building/Updating help pages for package 'probdistr'
> Formats: text html latex example
>  codeforplot   texthtmllatex   example
> ** building package indices ...
> * DONE (probdistr)
> * creating vignettes ... OK
> * removing junk files
> * checking for LF line-endings in source and make files
> * checking for empty or unneeded directories
> * building 'probdistr_1.00.tar.gz'
>
> Does someone have an idea about this?
>
> Thanks in advance
>
> Christophe
>
> PS : files are available on R-forge
>
> --
> Christophe Dutang
> Ph. D. student at ISFA, Lyon, France
> website: http://dutangc.free.fr
>
>
>
>

--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr





[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] Vignettes with missing or empty \VignetteIndexEntry:

[Christophe Dutang]

Hi,

I have a problem when checking the package 'probdistr' (on probability  
distributions).

I got this warning

* checking index information ... WARNING
Vignettes with missing or empty \VignetteIndexEntry:
 [1] "probdistr-chi""probdistr-contextra"  "probdistr-discrete"
 [4] "probdistr-discrextra" "probdistr-exp""probdistr- 
finitesupp"

 [7] "probdistr-gaussian"   "probdistr-general""probdistr-gumbel"
[10] "probdistr-intro"  "probdistr-logis"  "probdistr-math"
[13] "probdistr-misc"   "probdistr-mult"   "probdistr-pareto"
[16] "probdistr-student"
See the information on INDEX files and package subdirectories in the
chapter 'Creating R packages' of the 'Writing R Extensions' manual.

The warning tells me \VignetteIndexEntry is missing however there is  
one in probdistr-main.Rnw (Sweave file calling all other files via  
\SweaveInput{probdistr-X}). Let us note this file is not listed in  
the above files...


But the build is ok

* checking for file 'pkg/DESCRIPTION' ... OK
* preparing 'pkg':
* checking DESCRIPTION meta-information ... OK
* installing the package to re-build vignettes
* Installing *source* package ‘probdistr’ ...
** R
** inst
** help
*** installing help indices
 >>> Building/Updating help pages for package 'probdistr'
 Formats: text html latex example
  codeforplot   texthtmllatex   example
** building package indices ...
* DONE (probdistr)
* creating vignettes ... OK
* removing junk files
* checking for LF line-endings in source and make files
* checking for empty or unneeded directories
* building 'probdistr_1.00.tar.gz'

Does someone have an idea about this?

Thanks in advance

Christophe

PS : files are available on R-forge

--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] License status of CRAN packages

[Christophe Dutang]

Hi all,

I think for the common licences, we should also add BSD licence... for  
example my pkg randtoolbox (which is currently with incompatible  
licences) will probably be in a near future with the BSD licence.


Anyway I like the idea of two different repositories for GPL like  
licensed pkg and other packages.


Christophe

Le 24 avr. 09 à 18:20, Gabor Grothendieck a écrit :

On Fri, Apr 24, 2009 at 11:44 AM, Ben Goodrich > wrote:

Kurt Hornik wrote:

AGPL, unfortunately, allows supplements, and hence cannot fully be
standardized.  We've been thinking about extending the current  
scheme to

indicate a base license plus supplements, but this is still work in
progress.


This would be helpful. I would just reemphasize that a package that
includes some AGPL code and some GPL3 code is standard as far as  
the FSF

is concerned, e.g. from section 13 of the AGPL:

"Notwithstanding any other provision of this License, you have
permission to link or combine any covered work with a work licensed
under version 3 of the GNU General Public License into a single  
combined
work, and to convey the resulting work. The terms of this License  
will

continue to apply to the part which is the covered work, but the work
with which it is combined will remain governed by version 3 of the  
GNU

General Public License."

So, I think that CRAN should at least have a canonical spec that  
covers

*this* situation. Other situations may be more complicated to handle
elegantly.


Another possibility is to simply standardize the set of licenses  
that CRAN

supports.  GPL licenses (GPl-2, GPL-2.1, GPL-3, LGPL), MIT and
X11 already cover 98% of all packages on CRAN.   If there truly is an
advantage to the AGPL license perhaps a standard version could be  
offered
in the set.  Perhaps, for the 2% of packages that want a different  
license

a second repository could be made available.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] License question (RUnit)

[Christophe Dutang]

Hi,

Try to contact the other authors... I spent 5 min on google and found  
this http://www.uni-konstanz.de/FuF/Verwiss/koenig/


Regards

Christophe

Le 1 avr. 09 à 10:58, Pierre-Yves a écrit :


Dear list,

Sorry for the noise but I have a question regarding the license used  
in

RUnit [1], I contacted the maintainer( burgerm -at- users -dot-
sourceforge -dot- net ) on March 20th but I have received no answer.

Could anyone help to solve this question ?

Basically, my problem is that the website and the DESCRIPTION file say
that the license is GPLv2 while the header in the code says it is  
GPLv2

or any later version.

Thanks in advance for your help,

Best regards,

Pierre

[1] http://cran.r-project.org/web/packages/RUnit/index.html

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] dger_ in BLAS definition

[Christophe Dutang]

Yes x and y arguments are unchanged on exit, cf. 
http://www.mathkeisan.com/UsersGuide/man/dger.html

This is the work of the R core team to update those files, but I fear  
there are other functions which are not well declared.


Will you agree to take a look at the BLAS.h file? It will be very  
useful if you have time to do it.


(A year ago, I check the lapack.h and there was some wrong  
declarations like zgecon function.)


Christophe


Le 10 mars 09 à 22:49, Andrew Redd a écrit :


I'm developing some software and running into compiling warning:
conditionals.c:104: warning: passing argument 4 of 'dger_' discards
qualifiers from pointer target type
conditionals.c:104: warning: passing argument 6 of 'dger_' discards
qualifiers from pointer target type

the netlib documentation states that the arguments x and y should be
unchanged on exit.  Should should imply the defintion:

F77_NAME(dger)(const int * const m, const int * const n, const  
double *

const alpha,
  double const * const x, const int *const incx,
  double const * const y, const int *const incy,
  double * const a, const int * const lda);

the current definition is missing the appropriate consts:
F77_NAME(dger)(const int *m, const int *n, const double *alpha,
  double *x, const int *incx,
  double *y, const int *incy,
  double *a, const int *lda);

I don't want my code compiling with warnings that shouldn't be  
there.  Are

there suggestions of how to work around this?

Thanks

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] R tutorial

[Christophe Dutang]

Ok

I will stop sending you stuff I found on the web... but just asking  
the authors to contact us.


It is true that the page is quite old... So i will directly contact  
the author.


Thanks for your comments

Christophe

iPhone.fan

Le 22 févr. 09 à 23:13, Friedrich Leisch de> a écrit :



On Sun, 22 Feb 2009 14:44:27 +0100,
Christophe Dutang (CD) wrote:



Dear all,
I have just found a 'good' tutorial R for datamining. I think it
should be on the contributed docs.
http://cran.r-project.org/other-docs.html



Here is the link
http://www.liaad.up.pt/~ltorgo/DataMiningWithR/



What do you think?


Looks intersting (although may be outdated as Doug already
mentioned). In any case, we don't simply "harvest" docs off the web,
we only put docs on CRAN which get sent to us by the original
authors. So if you want to see it on CRAN, please contact the author
(and perhaps ask for the status while doing so).

Best,
Fritz

--
--- 


Prof. Dr. Friedrich Leisch

Institut für Statistik  Tel: (+49 89) 2180 3 
165
Ludwig-Maximilians-Universität  Fax: (+49 89) 2180 5 
308

Ludwigstraße 33
D-80539 München http://www.statistik.lmu.de/~lei 
sch
--- 


  Journal Computational Statistics --- http://www.springer.com/180
 Münchner R Kurse --- http://www.statistik.lmu.de/R
--- 





__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] R tutorial

[Christophe Dutang]

Dear all,

I have just found a 'good' tutorial R for datamining. I think it  
should be on the contributed docs.

http://cran.r-project.org/other-docs.html

Here is the link
http://www.liaad.up.pt/~ltorgo/DataMiningWithR/

What do you think?

Kind regards

Christophe


--
Christophe Dutang
Ph. D. student at ISFA, Lyon, France
website: http://dutangc.free.fr

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] percent concordant function

[christophe dutang]

Hi,

I would like to compute ''SAS'' output statistic for logistic regression,
namely the percent concordant, percent discordant, Sommer's D, Gamma,
Tau-a...Actually, I'm very interested in percent concordant.

Since the function 'cor' with method kendall compute the kendall's tau, I'm
wondering how can I use this function to get back the percent concordant?

Does anyone know if there is a package implementing those statistics? I look
around on the web, but did not find anything useful.


Thanks in advance

Christophe

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Can DESCRIPTION Maintainer: field contain general URL instead of only email address?

[christophe dutang]

There is a field URL in the DESCRIPTION file of the package. Could you use
it for route support request?


2008/10/3 William Dunlap <[EMAIL PROTECTED]>

> Our firm would like to route support requests through
> a website instead of using email.  However R will refuse
> to install a package if DESCRIPTION's Maintainer field
> does not have a valid email adress or the special value
> "ORPHANED".  (The check is done with
> tools:::.valid_maintainer_field_regexp.)
>
> I imagine a typical Maintainer line for the package "foo"
> might be something like
>   Maintainer: Amalgamated Widget's Support Team
> 
>
> How much does the package system (or CRAN) depend on
> the Maintainer field being an email address?  I can imagine
> programs would find mailing to an address easier than
> dealing with a website, but humans would have a better experience
> with a well designed website.  Perhaps we could add an optional
> %s field to the website so a program could post a small message.
>
> Bill Dunlap
> TIBCO Spotfire Inc.
> wdunlap tibco.com
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] package config file for windows

[Christophe Dutang]

Thanks for your answer.

You confirm what I fear, it is not easily possible to test for SSE2  
support on windows.


Can I assume there exists inttypes.h on windows platform?

Thanks again

Christophe

Le 14 sept. 08 à 01:11, Duncan Murdoch a écrit :


Christophe Dutang wrote:

Hi,

I'm maintaining randtoolbox package on CRAN and I wonder how to do  
a  windows config file? I need to test SSE2 instructions support as  
well  as inttypes.h library check.


Currently I use the trick of 'foreign' package, i.e. I have  
config.win  file with

cp -p src/config.h.win src/config.h
and config.h.win was written manually from config.h.in. There is  
no  test at all on windows.


In 'writing r extension' on this topic, we find
"
You should bear in mind that the configure script may well not work  
on  Windows systems
(this seems normally to be the case for those generated by  
Autoconf,  although simple shell scripts
do work). If your package is to be made publicly available, please   
give enough information for a
user on a non-Unix platform to configure it manually, or provide a   
‘configure.win’ script to be
used on that platform. (Optionally, there can be a  
‘cleanup.win’  script as well. Both should
be shell scripts to be executed by ash, which is a minimal version  
of  Bourne-style sh.)

"
Which tool do I need to write config script on windows? Do I need  
an  autoconf-cygwin solution? ( http://sources.redhat.com/autobook/autobook/autobook_242.html#SEC242 
   ) or just a port of autoconf available on sourceforge?


I think most configure scripts on Windows were written using a text  
editor, i.e. by hand. Windows systems are pretty consistent, so the  
kinds of tests and searches that you need to do on *nix aren't  
needed. You can assume that the Rtools are installed, but don't rely  
on anything else. We're unlikely to drop tools from that collection,  
so this advice should last quite a while.


Duncan Murdoch


Thanks in advance


Christophe
[[alternative HTML version deleted]]

  



__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel





__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] package config file for windows

[Christophe Dutang]

Hi,

I'm maintaining randtoolbox package on CRAN and I wonder how to do a  
windows config file? I need to test SSE2 instructions support as well  
as inttypes.h library check.

Currently I use the trick of 'foreign' package, i.e. I have config.win  
file with
cp -p src/config.h.win src/config.h
and config.h.win was written manually from config.h.in. There is no  
test at all on windows.

In 'writing r extension' on this topic, we find
"
You should bear in mind that the configure script may well not work on  
Windows systems
(this seems normally to be the case for those generated by Autoconf,  
although simple shell scripts
do work). If your package is to be made publicly available, please  
give enough information for a
user on a non-Unix platform to configure it manually, or provide a  
‘configure.win’ script to be
used on that platform. (Optionally, there can be a ‘cleanup.win’  
script as well. Both should
be shell scripts to be executed by ash, which is a minimal version of  
Bourne-style sh.)
"
Which tool do I need to write config script on windows? Do I need an  
autoconf-cygwin solution? ( 
http://sources.redhat.com/autobook/autobook/autobook_242.html#SEC242 
  ) or just a port of autoconf available on sourceforge?


Thanks in advance


Christophe
[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] non user-friendly error for chol2inv functions

[christophe dutang]

Yes, I do not cast the first argument as a matrix with as.matrix function.
Maybe we could detail the error message if the first argument is a numeric?

error(_("'a' is a numeric and must be coerced to a numeric matrix"));

Thanks for your answer


2008/8/29 Martin Maechler <[EMAIL PROTECTED]>

> >>>>> "cd" == christophe dutang <[EMAIL PROTECTED]>
> >>>>> on Fri, 29 Aug 2008 12:44:18 +0200 writes:
>
>cd> Hi,
>cd> In function chol2inv with the option LINPACK set to false (default),
> it
>cd> raises an error when the matrix is 1x1 matrix (i.e. just a real)
> saying
>
>cd> 'a' must be a numeric matrix
>
> It is very helpful, but you have to read and understand it.
> I'm pretty sure you did not provide a  1 x 1 matrix.
>
> Here's an example showing how things works :
>
> > m <- matrix(4,1,1)
> > cm <- chol(m)
> > cm
> [,1]
> [1,]2
> > chol2inv(cm)
> [,1]
> [1,] 0.25
> >
>
> Martin Maechler, ETH Zurich
>
>
>cd> This error is raised by the underlying C function (modLa_chol2inv in
>cd> function Lapack.c). Everything is normal, but I wonder if we could
> have
>cd> another behavior when we pass a 1x1 matrix. I spent time this
> morning
>cd> finding where was the error, and it was this "problem".
>
>cd> Thanks in advance
>
>cd> Christophe
>
>cd> [[alternative HTML version deleted]]
>
>cd> __
>cd> R-devel@r-project.org mailing list
>cd> https://stat.ethz.ch/mailman/listinfo/r-devel
>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] non user-friendly error for chol2inv functions

[christophe dutang]

Hi,

In function chol2inv with the option LINPACK set to false (default), it
raises an error when the matrix is 1x1 matrix (i.e. just a real) saying

'a' must be a numeric matrix

This error is raised by the underlying C function (modLa_chol2inv in
function Lapack.c). Everything is normal, but I wonder if we could have
another behavior when we pass a 1x1 matrix. I spent time this morning
finding where was the error, and it was this "problem".

Thanks in advance

Christophe

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] standard errors for glm

[christophe dutang]

Hi,

Thanks for answering so quickly.

Actually, I do not have the John Fox's book. On this webpage, there is a
handout for logistic regression (GLM I'm interested in) (cf.
http://socserv.socsci.mcmaster.ca/jfox/Courses/SPIDA/logistic-regression-handout.pdf).

The summary function return standard errors for coefficient estimates but
not for predictions. According to
http://www.nag.co.uk/numeric/cl/Manual/pdf/G02/g02gbc.pdf<http://www.nag.co.uk/numeric/cl/Manual/pdf/G02/g02gbc.pdfprovides>
the standard errors for the linear predictors X \hat{\beta} is X C X^t where
C is the variance/covariance matrix of coefficient estimate \hat{\beta} .
But I do not know standard errors for a response.

Christophe

2008/8/27 <[EMAIL PROTECTED]>

> Hi:  you can check John Fox's CAR book if you have it. I don't remember for
> sure but I may have  standard error related calculations for some of
> his graphics in the GLM section of the book ? But, Can you not get the
> sigma^2 hat  the summary of a GLM ? i would do
>
> glmsum<-summary("yourglmmodel")  and check there.
>
> I'd be really surprised if the variance of the error term wasn't there.
> Then, if you have that I think you can use that to calculate the prediction
> standard error as long as you assume that the parameters of the model are
> known with certainty. The formula  for that is in the regular regression (
> non GLM ) textbooks but i don't remember it off the top of my head. Good
> luck and hopefully someone else will reply with more exact info.
>
>
>
>
>
>
> On Wed, Aug 27, 2008 at  1:25 PM, christophe dutang wrote:
>
>  Hi,
>>
>> I'm currently using biglm package to compute GLM outputs on a very large
>> dataset. But no function computes standard erros of predictions. I look in
>> what is done in R, namely in the function predict.glm.R in stats package.
>> In this function, we call predict.lm to compute the standard errors (line
>> 51). The code of predict.lm (in lm.R) is very hard to understand.
>>
>> I wonder if there is any good reference and / documentation on this topic?
>> the manual at
>> http://www.nag.co.uk/numeric/cl/Manual/pdf/G02/g02gbc.pdfprovides
>> a good overview of the method used in R, but there is no reference
>> to standard errors...
>>
>> I suppose this topic have already raised in the past, but I found only
>> this
>> http://tolstoy.newcastle.edu.au/R/help/04/08/1762.html
>>
>> Thanks in advance
>>
>> Christophe
>>
>>[[alternative HTML version deleted]]
>>
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] standard errors for glm

[christophe dutang]

Hi,

I'm currently using biglm package to compute GLM outputs on a very large
dataset. But no function computes standard erros of predictions. I look in
what is done in R, namely in the function predict.glm.R in stats package.
In this function, we call predict.lm to compute the standard errors (line
51). The code of predict.lm (in lm.R) is very hard to understand.

I wonder if there is any good reference and / documentation on this topic?
the manual at http://www.nag.co.uk/numeric/cl/Manual/pdf/G02/g02gbc.pdfprovides
a good overview of the method used in R, but there is no reference
to standard errors...

I suppose this topic have already raised in the past, but I found only this
http://tolstoy.newcastle.edu.au/R/help/04/08/1762.html

Thanks in advance

Christophe

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [R] RNG Cycle and Duplication (PR#12540)

[Christophe Dutang]

Hi,

If you want I can put RANARRAY in the 'randtoolbox' package? unless  
you want in R and not in a package?


In randtoolbox (not the version on cran), I port SFMT and WELL  
generator respectively from Matsumoto and L'Ecuyer. It could be a good  
idea to add Knuth's code?


Christophe Dutang

Le 17 août 08 à 00:09, Shengqiao Li a écrit :



Knuth's double version RNG rng-double.c dose a great job. No ties  
were observed for 10M numbers ( totally 2^52 possible different  
values?). In rng-double, double modulo mod_sum replaced the integer  
version mod_diff in the integer version rng.c that is adopted by R.


The integer version uses modulus 2^30. Therefore there are only 2^30  
distinct numbers, which is confirmed by my previous test in R.


If someday Knuth's double version is also included in R, it will be  
great.



Shengqiao Li


On Fri, 15 Aug 2008, Duncan Murdoch wrote:


On 15/08/2008 10:28 AM, Shengqiao Li wrote:
Thank you for your reply and for your suggestion. So the note in  
man page could be more accurate since for an end user, man page  
should be more helpful and source code is mainly for developers.

I was also adviced to use Knuth's  double version RANARRAY from
http://www-cs-faculty.stanford.edu/~knuth/programs.html instead of  
the integer versions in R. I'm a R user. So why not also include  
the double verion in R implementation?


You can try it using kind="user-supplied" if you like, but I  
suspect it's the same as "Knuth-TAOCP-2002".


Duncan Murdoch


Thanks again,

Shengqiao Li
Research Associate
The Department of Statistics
PO Box 6330
West Virginia University
Morgantown, WV 26506-6330

On Fri, 15 Aug 2008, Duncan Murdoch wrote:

[EMAIL PROTECTED] wrote:
 This message is in MIME format.  The first part should be  
readable text,
 while the remaining parts are likely unreadable without MIME- 
aware tools.

---559023410-851401618-1218751024=:15885
Content-Type: TEXT/PLAIN; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: QUOTED-PRINTABLE
I didn't describe the problem clearly. It's about the number of  
distinct=20

values. So just ignore cycle issue.
My tests were:
RNGkind(kind=3D"Knuth-TAOCP");
sum(duplicated(runif(1e7))); #return 46552
RNGkind(kind=3D"Knuth-TAOCP-2002");
sum(duplicated(runif(1e7))); #return 46415
#These collision frequency suggested there were 2^30 distinct  
values by=20

birthday problem.
The birthday problem distribution applies to independent draws,  
but they are only pseudo-independent.  I think the only ways to  
know for sure if there are 2^30 values are to look at the code,  
or run through a complete cycle.  And you need to determine the  
cycle by looking at .Random.seed, not at the returned value.

RNGkind(kind=3D"Marsaglia-Multicarry");
sum(duplicated(runif(1e7))); #return 11682
RNGkind(kind=3D"Super-Duper");
sum(duplicated(runif(1e7))); #return 11542
RNGkind(kind=3D"Mersenne-Twister");
sum(duplicated(runif(1e7))); #return 11656
#These indicated there were 2^32 distinct values, which agrees  
with the=20

help info.
If there are 2^30 distinct values for the two generators above,  
that also agrees with the documentation.

RNGkind(kind=3D"Wichmann-Hill");
sum(duplicated(runif(1e7))); #return 0
#So for this method, there should be more than 2^32 distinct  
values.
You may not get the exact numbers, but they should be close. So  
how to=20

explain above problem?
You haven't demonstrated what you claim, but if you look at the  
source, you'll see that in fact the man page is wrong:  Wichmann- 
Hill is based on 3 integer values, which each take on  
approximately 15 bits of different values. So Wichmann-Hill could  
take nearly 2^45 different values (actually 30269*30307*30323).
The source is in https://svn.r-project.org/R/trunk/src/main/RNG.c  
if you want to check the others.

I need generate a large sample without any ties, it seems to me=20
"Wichmann-Hill" is only choice right now.
An alternative would be to construct a new value from two (or  
more) runif() values, but be careful that you don't mess up the  
distribution when you do that.

Duncan Murdoch
= 
3D 
= 
3D 
= 
3D 
= 
3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
Shengqiao Li
The Department of Statistics
PO Box 6330
West Virginia University
Morgantown, WV 26506-6330
= 
3D 
= 
3D 
= 
3D 
= 
3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=

=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
On Thu, 14 Aug 2008, Peter Dalgaard wrote:

Shengqiao Li wrote:

Hello all,
=20
I am generating large samples of random numbers. The RNG help  
page says:=

=20
"All the supplied uniform generators return 32-bit integer  
values that a=

re=20
converted to double

Re: [Rd] [R] aligned memory allocation in C

[christophe dutang]

Thanks for your remarks.

According to
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/howto-compile.html, I
need 16 bit aligned memory when using fill_array64. So I suppose I need 8
bit aligned memory. I will test what you advise me and will come back to
R-devel list after.

Thanks again

Christophe

2008/8/13 Prof Brian Ripley <[EMAIL PROTECTED]>

> On Tue, 12 Aug 2008, Jeffrey Horner wrote:
>
>  Christophe Dutang1 wrote:
>>
>>> Hi,
>>>
>>> I'm currently R porting SF Mersenne Twister algorithm of Matsumoto and
>>> Saito. To get the full power of their code, I want to use their fonction
>>> fill_array32 which need aligned memory. That is to say I need to use the C
>>> function memalign on windows, posix_memalign on linux and classic malloc on
>>> Mac OS. In 'writing R extenstion', they recommand to use R_alloc function to
>>> allocate memory in C.
>>>
>>> Does R_alloc return a pointer to aligned memory?
>>> if not how can I do this?
>>> probably no, because R crashes when I succesively R_alloc and
>>> fill_array32 (cf below) on my macbook with R 2.7.1.
>>>
>>
>> You can still do this. Just take the address returned from R_alloc and
>> test for alignment. If it's not, then just use an aligned address beyond the
>> one returned.
>>
>
> We haven't been told what the desired alignment is (and those functions
> need to be told).  On 32-bit Mac OS X, R_alloc is definitely aligned on
> 4-byte boundaries (on 64-bit OSes it is usually 8-byte aligned).
>
>  (But then the question is, which direction beyond the one returned? How
>> does one test for that?)
>>
>
> Addresses always go upwards.  So if you want 64-byte alignment you need to
> allocate a block at least 64 bytes longer than required, and go up to the
> nearest multiple of 64.
>
> BTW, this is clearly an R-devel question -- see the posting guide.
>
>
>
>> Jeff
>>
>>  Thanks in advance
>>>
>>> Kind regards
>>>
>>> Christophe
>>>
>>>
>>> PS :
>>> http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/howto-compile.htmlprovides
>>>  an example of memalign.
>>>
>>> PPS : mac os report
>>>
>>
> [removed]
>
> --
> Brian D. Ripley,  [EMAIL PROTECTED]
> Professor of Applied Statistics,  
> http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] [R] aligned memory allocation in C

[Christophe Dutang]

Hi,

You are completely right, the problem was not about what R_alloc  
returned but my use of the fill_array32. Now it works fine, but the  
block generation is longer due to a conversion to double after using  
fill_array32. I must underline the fact that I do not yet use sse2  
support.


Thanks again

Christophe


Le 13 août 08 à 14:30, Luke Tierney a écrit :


On Wed, 13 Aug 2008, Christophe Dutang1 wrote:


Hi,

I'm currently R porting SF Mersenne Twister algorithm of Matsumoto  
and Saito. To get the full power of their code, I want to use their  
fonction fill_array32 which need aligned memory. That is to say I  
need to use the C function memalign on windows, posix_memalign on  
linux and classic malloc on Mac OS. In 'writing R extenstion', they  
recommand to use R_alloc function to allocate memory in C.


Does R_alloc return a pointer to aligned memory?
if not how can I do this?
probably no, because R crashes when I succesively R_alloc and  
fill_array32 (cf below) on my macbook with R 2.7.1.


R_alloc's alignment will be appropriate for holding any data type. It
will be offset from a value returned by malloc by a multiple of 8
bytes.

My recollection, which may be wrong, is that on both Intel and PPC
unaligned access to all basic data types is permitted but may be
inefficient (in particular on Intel), so the reason for your crash is
probably elsewhere.

Best,

luke




Thanks in advance

Kind regards

Christophe


PS : http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/howto-compile.html 
 provides an example of memalign.


PPS : mac os report



__
[EMAIL PROTECTED] 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.


--
Luke Tierney
Chair, Statistics and Actuarial Science
Ralph E. Wareham Professor of Mathematical Sciences
University of Iowa  Phone: 319-335-3386
Department of Statistics andFax:   319-335-3017
  Actuarial Science
241 Schaeffer Hall  email:  [EMAIL PROTECTED]
Iowa City, IA 52242 WWW:  http://www.stat.uiowa.edu


__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel