[julia-users] ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)

[Maxwell]

Hello folks,

I am trying to run a little bootstrap in Julia using pmap and I am getting 
the following error from each process:

fatal error on fatal error on fatal error on fatal error on fatal error on 
fatal error on fatal error on fatal error on fatal error on fatal error on 
fatal error on fatal error on fatal error on fatal error on fatal error on 
2141612: : : : 46: : 8: 91071131315: : : : : : : 5: ERROR: `convert` has no 
method matching convert(::Type{Int64...}, ::Int64)
ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)
ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)
ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)
ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)

However, running the code for one case works just perfectly. I am been 
scratching my head about this for few days now. Any help will be greatly 
appreciated.
Find attached my code.

Thanks.




  _
   _   _ _(_)_ |  A fresh approach to technical computing
  (_) | (_) (_)|  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "help()" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.3.6
 _/ |\__'_|_|_|\__'_|  |  
|__/   |  x86_64-redhat-linux




Boot_code.jl
Description: Binary data

[julia-users] re ERRoR `convert` has no method matching convert(::Type{Int64...}, ::Int64) in pmap in Julia 0.3.6

[Maxwell]

Hello folks,

This is a re-post from my previous post earlier in the week. My code 
 was referring t a non-existent file. I am attaching the code I am using 
since I have no idea what part of the code is producing errors in Julia 
0.3.6.

 
I am trying to run a little bootstrap in Julia using pmap and I am getting 
the following error:
 

fatal error on fatal error on 5: 6: ERROR: `convert` has no method matching 
convert(::Type{Int64...}, ::Int64)
ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)
Worker 6 terminated.
Worker 5 terminated.

julia> versioninfo()
Julia Version 0.3.6
Platform Info:
  System: Linux (x86_64-redhat-linux)
  CPU: Intel(R) Xeon(R) CPU E5-2690 0 @ 2.90GHz
  WORD_SIZE: 64
  BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblasp.so.0
  LIBM: libopenlibm
  LLVM: libLLVM-3.4.2

backtrace()
3-element Array{Ptr{None},1}:
 Ptr{Void} @0x003c90c7ed25
 Ptr{Void} @0x003c90c7ed80
 Ptr{Void} @0x7fbf138f4010


However, when I run the same code on my personal computer in Julia 0.3.3 it 
works fine:

julia> include("Run_Boot_code.jl")

Warning: could not import Base.add! into NumericExtensions
Warning: could not import Base.add! into NumericExtensions
From worker 2:  Done!
From worker 3:  Done!
Done!!
2x12 DataFrame
| Row | n1   | n2   | p | k   | m| kap0 | pct | mu  | Log_bay  | 
pos |
|-|--|--|---|-|--|--|-|-|--|-|
| 1   | 50.0 | 50.0 | 200.0 | 0.0 | 24.0 | 0.1  | 1.0 | 1.0 | -58.5275 | 
0.0 |
| 2   | 50.0 | 50.0 | 200.0 | 1.0 | 24.0 | 0.1  | 1.0 | 1.0 | -58.2967 | 
0.0 |

| Row | stron | ver_strg |
|-|---|--|
| 1   | 0.0   | 0.0  |
| 2   | 0.0   | 0.0  |

julia> versioninfo()
Julia Version 0.3.3
Commit b24213b* (2014-11-23 20:19 UTC)
Platform Info:
  System: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-4570 CPU @ 3.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3



However, running the code for one case works just perfectly on both 
version. I am been scratching my head about this for few days now. Any help 
will be greatly appreciated.


Find attached my code. Sorry about the length of my code. 

Thanks.


Boot_code.jl
Description: Binary data


Run_Boot_code.jl
Description: Binary data

[julia-users] Non-central F distribution

[Maxwell]

Hi all,
I would like to simulate from a non-central F distribution. I don't think 
the package Distributions supports non-central F or the non-central Chi 
squared. Is there any other package(s) that can help me with that? or does 
anybody have an idea how that can be done?

Thanks,  

Re: [julia-users] Non-central F distribution

[Maxwell]

Thanks Andreas! That's awesome!

On Wednesday, December 17, 2014 12:23:54 AM UTC-6, Andreas Noack wrote:
>
> You might be right that it is not documented, but they are supported. Have 
> a look in src/univariate/continuous
>
> julia> using Distributions
>
> julia> d1 = NoncentralF(1.2,2.3,3.4)
> NoncentralF(ndf=1.2, ddf=2.3, ncp=3.4)
>
> julia> rand(d)
> 33.94546618024409
>
> julia> d2 = NoncentralChisq(3, 1.1)
> NoncentralChisq(df=3.0, ncp=1.1)
>
> julia> rand(d2)
> 4.822266647364242
>
> 2014-12-17 2:31 GMT+01:00 Maxwell >:
>>
>>
>> Hi all,
>> I would like to simulate from a non-central F distribution. I don't think 
>> the package Distributions supports non-central F or the non-central Chi 
>> squared. Is there any other package(s) that can help me with that? or does 
>> anybody have an idea how that can be done?
>>
>> Thanks,  
>>
>

Re: [julia-users] Non-central F distribution

[Maxwell]

I agree Tim! I will work on that.

On Wednesday, December 17, 2014 5:13:50 AM UTC-6, Tim Holy wrote:
>
> Maxwell, it would be great if you could submit a pull request to add that 
> documentation. 
>
> Best, 
> --Tim 
>
> On Wednesday, December 17, 2014 07:23:50 AM Andreas Noack wrote: 
> > You might be right that it is not documented, but they are supported. 
> Have 
> > a look in src/univariate/continuous 
> > 
> > julia> using Distributions 
> > 
> > julia> d1 = NoncentralF(1.2,2.3,3.4) 
> > NoncentralF(ndf=1.2, ddf=2.3, ncp=3.4) 
> > 
> > julia> rand(d) 
> > 33.94546618024409 
> > 
> > julia> d2 = NoncentralChisq(3, 1.1) 
> > NoncentralChisq(df=3.0, ncp=1.1) 
> > 
> > julia> rand(d2) 
> > 4.822266647364242 
> > 
> > 2014-12-17 2:31 GMT+01:00 Maxwell >: 
> > > Hi all, 
> > > I would like to simulate from a non-central F distribution. I don't 
> think 
> > > the package Distributions supports non-central F or the non-central 
> Chi 
> > > squared. Is there any other package(s) that can help me with that? or 
> does 
> > > anybody have an idea how that can be done? 
> > > 
> > > Thanks, 
>
>

[julia-users] Issues using the R package corpcor with RCall in Julia

[Maxwell]

Hi folks,

I am trying to play around with RCall by calling some R packages from 
Julia. Then I run into this issue (see below)

using RCall
@rimport highmean as highmean ## Import the package highmean from r
@rusing highmean
srand(1234)
Y1 = randn(5,10)
Y2 = rand(5,10)
julia> pval  = rcopy(apval_Bai1996(Y1,Y2))[:pval]
0.026691441016612916

# Now I try to obtain a shrinkage estimate of the common covariance matrix 
based on Y1 and Y2

@rimport corpcor as corpcor ## Import the package highmean from r
@rusing corpcor

julia> sig = rcopy(corpcor.shrink.cov(vcat(Y1,Y2)))
ERROR: UndefVarError: shrink not defined


 corpcor is installed in R on my machine.

Thanks

[julia-users] Re: emacs IJulia notebook?

[Maxwell Peterson]

I have just spent a few hours getting past this error and this thread is a 
top google result, so I am adding my fix:

I had the same issue as Andrew. First, note that the tkf repo hasn't been 
updated for years, and the current emacs-ipython-notebook repo is 
millejoh's fork: https://github.com/millejoh/emacs-ipython-notebook. His is 
the current MELPA package, so under emacs 24, use package-list-packages and 
install "ein" to get emacs-ipython-notebook.

Second, I had an old version of ipython. I uninstalled all the ipython on 
my system and reinstalled with pip: pip install ipython (giving me version 
4.0.0).

Then it worked.

On Wednesday, April 16, 2014 at 2:53:54 PM UTC-5, Andrew Dabrowski wrote:
>
> I'm trying iPython 2.0, but now I'm getting this error:
>
> https://github.com/tkf/emacs-ipython-notebook/issues/137
>
> That seems to still be an open issue with iPython 2.  I get the same error 
> whether I use stable emacs 24 or a bleeding edge version.
>
>
>

[julia-users] Re: how to print out part of a row in DataFrame?

[Taylor Maxwell]

If your dataframe has 6 or more columns then it prints this summary.  If 
you index less than 6 columns it will print ok.  You can force it to spit 
it all out with showall(df) but that may not be entirely legible.   There 
may be other ways to display it.

[julia-users] Re: how to print out part of a row in DataFrame?

[Taylor Maxwell]

scratch what I just said

[julia-users] Re: ANN: EM algorithm for regularized L0 regression. More regreg Julia code out there?

[Taylor Maxwell]

It looks like there is a group effort currently underway to make a nice 
interface for much of the machine learning and statistics problems out 
there.  This is primarily being headed up by Dahua Lin.  Your work 
definitely fits within the scope they have talked about. Here are some 
links:

ttp://juliastats.github.io

https://github.com/JuliaStats/Roadmap.jl/issues/14

https://github.com/JuliaStats/Roadmap.jl/issues/11

https://github.com/JuliaStats/RegERMs.jl/issues/3

[julia-users] Re: Ordinary least squares regression

[Taylor Maxwell]

I haven't seen code to do it yet but it is very simple to calculate with a 
LinearModel from GLM  Below is some code to calculate r-squared or adjusted 
r-squared from a linear model in GLM calculated from a dataframe or from a 
model calculated without a data frame.  At the bottom is some simple code 
to apply these functions.  It uses the Abs2Fun() from the NumericFuns 
package to calculate the sum of squares quickly.

using DataFrames
using GLM

import DataFrames.DataFrameRegressionModel
import GLM.LinearModel
import GLM.residuals
using NumericFuns


function 
rsquared(mod::DataFrameRegressionModel{LinearModel{DensePredQR{Float64}},Float64})
SStot=sum(Abs2Fun(),mod.model.rr.y)
SSres=sum(Abs2Fun(),residuals(mod))
return (1-(SSres/SStot))
end

function rsquared(mod::LinearModel{DensePredQR{Float64}})
SStot=sum(Abs2Fun(),mod.rr.y)
SSres=sum(Abs2Fun(),residuals(mod))
return (1-(SSres/SStot))
end

function 
adjrsquared(mod::DataFrameRegressionModel{LinearModel{DensePredQR{Float64}},Float64})
SStot=sum(Abs2Fun(),mod.model.rr.y)
SSres=sum(Abs2Fun(),residuals(mod))
n=size(mod.model.rr.y,1)  #number of samples
p=size(mod.mm.m,2)-1  #number of variables besides constant (assumes 
intercept)
return 1- ( (SSres/(n-p-1)) / (SStot/(n-1)) )
end

function adjrsquared(mod::LinearModel{DensePredQR{Float64}})
SStot=sum(Abs2Fun(),mod.rr.y)
SSres=sum(Abs2Fun(),residuals(mod))
n=size(mod.rr.y,1)#number of samples
p=size(mod.pp.X,2)-1  #number of variables besides constant (assumes 
intercept)
return 1- ( (SSres/(n-p-1)) / (SStot/(n-1)) )
end


dd=randn(1000,7);
df=convert(DataFrame,dd);
aa=fit(LinearModel,x1~x2+x3+x4+x5+x6,df)
bb=fit(LinearModel,aa.mm.m,aa.model.rr.y)

rsquared(aa)
adjrsquared(aa)

rsquared(bb)
adjrsquared(bb)


On Wednesday, August 6, 2014 5:33:15 PM UTC-6, Jason Solack wrote:
>
> Hello everyone.  I am looking for a way to run this regression and 
> retrieve the corresponding r squared value.  Is there a package out there 
> that does this?
>
> Thank you for the help!
>
> Jason
>

[julia-users] Re: Ordinary least squares regression

[Taylor Maxwell]

Sorry, I threw it up quickly, I should have spent more time making sure it 
was correct.

On Thursday, August 7, 2014 8:17:22 AM UTC-6, Simon Kornblith wrote:
>
>
>

Re: [julia-users] Has anyone successfully performed probit or logit regression in Julia?

[Taylor Maxwell]

Are you looking for the fitted values?  Is predict(OLS) what you are 
looking for?

*julia> **X = [1;2;3.]*

*3-element Array{Float64,1}:*

* 1.0*

* 2.0*

* 3.0*


*julia> **Y = [1;0;1.]*

*3-element Array{Float64,1}:*

* 1.0*

* 0.0*

* 1.0*


*julia> **data = DataFrame(X=X,Y=Y)*

*3x2 DataFrame*

*|---|-|-|*

*| Row # | X   | Y   |*

*| 1 | 1.0 | 1.0 |*

*| 2 | 2.0 | 0.0 |*

*| 3 | 3.0 | 1.0 |*


*julia> **OLS = glm(Y~X,data,Normal(),IdentityLink())*

*DataFrameRegressionModel{GeneralizedLinearModel,Float64}:*


*Coefficients:*

* Estimate Std.Error  z value Pr(>|z|)*

*(Intercept)  0.67   1.24722 0.534522   0.5930*

*X-4.16334e-16   0.57735 -7.2e-16   1.*



*julia> **predict(OLS)*

*3-element Array{Float64,1}:*

* 0.67*

* 0.67*

* 0.67*

[julia-users] Nested parametric types

[David Maxwell]

I ran into the following when tracking down a type instability issue. 
 Consider the following type hierarchy:

```julia

type T1{I<:Number} end

type T2{T<:T1} 
x::T
end

type T3{T<:T2} 
x::T
end

type T4{T<:T3} 
x::T
end

type T5{T<:T4} 
x::T
end
```
Now, two functions, one which returns a `T4`, and one which returns a `T5`.

```julia
function foo() 
t2 = T2( T1{Float64}() )
t3 = T3( t2 )
T4( t3 )
end

function bar() 
t2 = T2( T1{Float64}() )
t3 = T3( t2 )
t4 = T4( t3 )
T5( t4 )
end
```

Julia is able to infer the return type of `foo`, but not of `bar`.

```
julia> @code_warntype foo()
Variables:
  t2::T2{T1{Float64}}
  t3::T3{T2{T1{Float64}}}

Body:
  begin  #
  t2 = $(Expr(:new, T2{T1{Float64}}, :($(Expr(:new, T1{Float64}) # 
line 22:
  t3 = $(Expr(:new, T3{T2{T1{Float64}}}, :(t2::T2{T1{Float64}}))) # 
line 23:
  return $(Expr(:new, T4{T3{T2{T1{Float64, 
:(t3::T3{T2{T1{Float64}}})))
  end::T4{T3{T2{T1{Float64
Variables:
  t2::T2{T1{Float64}}
  t3::T3{T2{T1{Float64}}}
  t4::T4{T3{T2{T1{Float64

julia> @code_warntype bar()
Body:
  begin  #
  t2 = $(Expr(:new, T2{T1{Float64}}, :($(Expr(:new, T1{Float64}) # 
line 28:
  t3 = $(Expr(:new, T3{T2{T1{Float64}}}, :(t2::T2{T1{Float64}}))) # 
line 29:
  t4 = $(Expr(:new, T4{T3{T2{T1{Float64, 
:(t3::T3{T2{T1{Float64}}}))) # line 30:
  GenSym(0) = t4::T4{T3{T2{T1{Float64
  return 
((top(apply_type))(T5,T4{T3{T2{T1{Float64)::TYPE{_<:T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:NUMBER}})(GenSym(0))::T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:NUMBER}
  end::T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:NUMBER}
```
Am I running up against a genuine limit in Julia's type inference?

Re: [julia-users] Re: Nested parametric types

[David Maxwell]

Thanks everyone for the replies.

Jiahao:  Thanks especially for the pointer to MAX_TYPE_DEPTH. This is 
exactly what I hit, and now I'm trying to learn what its consequences are.

I appreciate the need for constraints to keep type inference efficient, but 
I'm not sure why what I am doing is running up against them. The root 
problem seems to be at the level of T5's constructor.  With the same type 
hierarchy as in the original post, constructing a T4 is not a problem.

*julia> **@code_warntype T4(T3(T2(T1{Float64}(*

Variables:

  x::T3{T2{T1{Float64}}}


Body:

  begin 

  return $(Expr(:new, T4{T3{T2{T1{Float64, 
:(x::T3{T2{T1{Float64}}})))

  end::T4{T3{T2{T1{Float64


But constructing a T5 is a problem, even if I tell the compiler exactly 
what kind of T5 I am constructing.

*julia> **t4 = T4(T3(T2(T1{Float64}(;*


*julia> **@code_warntype T5(t4)*

Variables:

  x::T4{T3{T2{T1{Float64


Body:

  begin 

  return ((top(apply_type))(T5,T)
*::Type{_<:T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:Number}}*
)(x::T4{T3{T2{T1{Float64)*::T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:Number}*

  end*::T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:Number}*

*julia> **@code_warntype T5{T4{T3{T2{T1{Float64}(t4)*

Variables:

  #s2::Type{T5{T4{T3{T2{T1{Float64}}

  x::T4{T3{T2{T1{Float64


Body:

  begin 

  return $(Expr(:new, 
:((top(apply_type))((top(getfield))(Main,:T5)::Type{T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:Number}},T)::Type{_<:T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:Number}}),
 
:(x::T4{T3{T2{T1{Float64)))

  end*::T5{T<:T4{T<:T3{T<:T2{T<:T1{I<:Number}*

I would have thought in this last case there was nothing left to infer.  Is 
this last case expected behavior? Is there a way to construct a T5 such 
that the type is known to the type inference system? (Without bumping up 
MAX_TYPE_DEPTH; I'd like to learn to live within the default limits!)

Without delving too far into inference.jl, I had worried that maybe type 
inference just punts when it encounters types exceeding MAX_TYPE_DEPTH. 
 But that seems to be not true.  The functions I've attempted seems to 
reason along happily about T5's once they are constructed.  Here's a 
trivial example.

*julia> **function id(x::T5) x end*

*id (generic function with 1 method)*


*julia> **t5 = T5(t4);*


*julia> **@code_warntype id(t5)*

Variables:

  x::T5{T4{T3{T2{T1{Float64}


Body:

  begin  # none, line 1:

  return x::T5{T4{T3{T2{T1{Float64}

  end::T5{T4{T3{T2{T1{Float64}

More interesting examples require modifying the types in the hierarchy so 
that they are more complex, but the few things I've tried (without pushing 
too hard) seem ok, once a T5 is in hand.

Anyway, in practice I find I use the idiom

abstract AbstractClass

type C{A<:AbstractClass}
  x::A
  more fields
end

so that the explicit type of field x is known if the type of C is known. 
 But if I understand the above correctly, it seems that C's constructor is 
not type stable, in the sense that there are constructors for kinds of C's 
such that the return type is not known to the type system, even when all 
the argument types to the constructor are known. Is this correct?

Many thanks,

David

[julia-users] Re: Problem iterating through GLMs in Monte Carlo experiments

[Taylor Maxwell]

You do not need DataFrames to run a GLM.  If you construct your now Float64 
design matrix (X) and have you dependent variable y you  can fit a GLM with 
such a logistic regression with the following command:

fit(GeneralizedLinearModel,X,y,Binomial())

to fit an LM model you can:

fit(LinearModel,X,y)

If Y is a matrix of dependent variables you can access it by y[:,i]

I hope this helps.

On Saturday, February 21, 2015 at 4:37:27 AM UTC-7, Steve Kay wrote:
>
> I'm trying to program a simple Monte Carlo experiment that after 
> generating sample data within each MC rep then runs a GLM within it, 
> accumulating the coefficients from each rep in a matrix. As the GLM command 
> needs a dataframe I convert all my matrix held simulated data into this 
> form and then try and run a for loop (see end of included code below). It 
> does not like the $i notation. Can anyone please come up with a workaround? 
> I had a quick look in the glm code to see if I could figure out how to 
> create a modelmatrix myself and call the glm directly on that - it defeated 
> me (as I'm not the greatest of programmers). If anyone could offer advice 
> on that too I'd be grateful - would save the glm command having to convert 
> the dataframe on each call - which on a large number of MC reps might be 
> worth while. 
>
> Any help much appreciated (and apologies for all the code comments - my 
> memories terrible so need to include).
>
> Best,
>
> Steve
>
> using DataFrames, GLM, Distributions
> reps=100 # number of MC sims
> n=50 # sample size in each sim
> X0 = rep(1,n) # the constant
> nLevelsX1 = 5 # distinct values in continuous X1 - make it a factor of n
> X1=gl(nLevelsX1,int(n/nLevelsX1)) # gl() is in DataFrames
> X2 = gl(2,1,n) - 1 # the treatmnt dummy
> X = [X0 X1 X2] # design matrix
> T=1.0 # the treatment effect
> B = [10, 1, T] # Beta vector
> mu = float(X*B) # varying mu parameter in lognormal distribution
> sigma = rep(1.0,n) # constant sigma parameter in lognormal distribution
>
> Obs = Array(Float64, reps,n) # holds data for analysis - each row new 
> dataset
> # loop below creates the Y (dependent) data - one row for each MC sim
> for i = 1:n
>   Obs[:,i]= rand(LogNormal(mu[i],sigma[i]),reps)
> end
>
>
> # create dataframe where first two cols hold the unvarying X
> # and the remaining cols hold each Y sim
> df= convert(DataFrame, [X[:,2:3] Obs'])
> rename!(df.colindex, [(symbol("x$i")=>symbol("Y$(i-2)")) for i in 3:(n+2)])
> names(df) # x1 x2 Y1 Y2 ...Yreps
>
> GLMresults =  Array(Float64, reps,3) # will hold all coeff
> #next loop fails - please help - won't accept $i
> for i = 1:reps
>   GLMresults[i,:] = coef(glm(Y$i ~ x1 + x2,df,Normal(),LogLink()))
> end
>

[julia-users] reconstructing a return type

[David Maxwell]

Suppose I have a heterogeneous-valued dictionary *foo_cache* where the keys 
are a tuple of types (e.g. *(S,T,R)*) and the values are the results of a 
call to a function *foo(x,S,T,R)* where the type of object *x* is known to 
type inference.  I would like to annotate the type of a call to 
*foo_cache[(S,T,R)]* so that it is known to type inference as well. I.e. I 
want to write

*f = foo_cache[(S,T,R)]::F*

where *F* is the return type of a call to *foo(x,S,T,R)*.  In effect, I 
want to tell type inference that *f* has the type of whatever type 
inference thinks a call to *foo(x,S,T,R)* would be. Is there a way to do 
this?

[julia-users] Re: Any ANOVA or ANCOVA examples for Julia?

[Taylor Maxwell]

You should be able to perform an anova or ancova with the GLM package.

https://github.com/JuliaStats/GLM.jl

You can use the lm() function and your variables with factors need to be 
PooledDataArrays in the dataframe.


On Wednesday, February 12, 2014 2:57:08 PM UTC-6, Edward Stembler wrote:
>
>
>
> Does anyone know of any ANOVA or ANCOVA examples for Julia out there?  I 
> haven't been able to find anything via Google.
>

[julia-users] Re: Any ANOVA or ANCOVA examples for Julia?

[Taylor Maxwell]

I just added a pull request with code to the GLM package to add an 
anovatable and an effects function. It gets the effects from a LmMod 
(created from lm()) and makes a standard anova table with term, df, SS, MS, 
F, and p-value.

On Wednesday, February 12, 2014 1:57:08 PM UTC-7, Edward Stembler wrote:
>
>
>
> Does anyone know of any ANOVA or ANCOVA examples for Julia out there?  I 
> haven't been able to find anything via Google.
>

[julia-users] Re: Julians in Colorado?

[Taylor Maxwell]

I just moved to Fort Collins in February (~hour north of Denver).  I am not 
a programming aficionado but I love Julia and it has opened my world of 
possibilities. I primarily do human statistical genetics.  My git name is 
timema.  Thanks for posting about the Spark/Julia night, I had not heard of 
it before.  If I have time, I might try to come down (I am in Virginia all 
of this week).  It would be great to meet more people from the Julia 
community, especially those in the local area.

On Monday, October 27, 2014 6:37:37 PM UTC-4, Sean Garborg wrote:
>
> I just moved from Chicago and am curious.
>
> Relatedly, a data science meetup in Broomfield, CO has a half Spark / 
> half Julia night 
>  
> coming up Wednesday Nov 5th, with Galen O'Neil running the Julia half. No 
> idea what the crowd is like, but I'm sure it wouldn't hurt to show up.
>

[julia-users] Re: Packages for handling pedigrees?

[Taylor Maxwell]

Pedigrees sounds like a good name. Will this package be able to generate 
and manipulate standard kinship matrices?  This is not my specialty but 
familial data is pretty common and many standard association methods for 
"unrelated" individuals have analogous method that need to take into about 
the relationships among individuals that come with familial data.  In 
particularly the the kinship matrix is often a large sparse matrix that is 
used as random effects in mixed linear models.  The kinship package in R is 
very useful for the construction and use of these matrices.  With some 
colleagues (who know the statistical theory better than I) we have 
developed a methodology to look at simultaneous mean and variance effects 
in unrelated individuals.  We recently adapted it to use family data.  We 
were able to use the kinship package in R to construct the kinship matrix 
(a sparse matrix) and applied our methodology with the help of their 
"gchol" function to implement a generalized Cholesky decomposition and used 
it along with BFGS to optimize the likelihood function.  I have implemented 
the "unrelated" method in Julia and would love to be able to implement this 
family based method in Julia as well.




On Tuesday, December 16, 2014 10:03:06 AM UTC-7, Douglas Bates wrote:
>
> Is anyone working on a Julia package or packages to work with pedigrees, 
> in the sense used in animal breeding?  Ana Vazquez and I are developing 
> such a package which we plan to call Pedigrees.  Its principle type will be 
> Pedigree, which is constructed from two Integer vectors representing the 
> 'sire' and 'dam' of each animal.  These names reflect the animal breeding 
> orientation of the package.
>
> Capabilities include sorting a pedigree so that parents occur before 
> offspring and deriving the additive relationship matrix and other 
> quantities related to it, such as inbreeding coefficients.
>
> Is using the name Pedigrees going to cause conflicts?
>

[julia-users] Re: Packages for handling pedigrees?

[Taylor Maxwell]

Sounds great, I remember the brief discussion we previously had about 
pedigreemm.  Thanks for working on this!


On Tuesday, December 16, 2014 12:57:14 PM UTC-7, Douglas Bates wrote:
>
> On Tuesday, December 16, 2014 1:32:19 PM UTC-6, Taylor Maxwell wrote:
>>
>> Pedigrees sounds like a good name. Will this package be able to generate 
>> and manipulate standard kinship matrices?  This is not my specialty but 
>> familial data is pretty common and many standard association methods for 
>> "unrelated" individuals have analogous method that need to take into about 
>> the relationships among individuals that come with familial data.  In 
>> particularly the the kinship matrix is often a large sparse matrix that is 
>> used as random effects in mixed linear models.  The kinship package in R is 
>> very useful for the construction and use of these matrices.  With some 
>> colleagues (who know the statistical theory better than I) we have 
>> developed a methodology to look at simultaneous mean and variance effects 
>> in unrelated individuals.  We recently adapted it to use family data.  We 
>> were able to use the kinship package in R to construct the kinship matrix 
>> (a sparse matrix) and applied our methodology with the help of their 
>> "gchol" function to implement a generalized Cholesky decomposition and used 
>> it along with BFGS to optimize the likelihood function.  I have implemented 
>> the "unrelated" method in Julia and would love to be able to implement this 
>> family based method in Julia as well.
>>
>
> Ana and I are familiar with using pedigree information in mixed-effects 
> models; we wrote the R package pedigreemm and I am one of the authors of 
> the lme4 package for R.  The purpose of creating tools for pedigrees in 
> Julia is to use Pedigrees and MixedModels to fit such models using Julia.
>
> I should look again at the kinship package for R.  (I get twitchy every 
> time I start to read it because Terry chose to reimplement sparse matrix 
> methods from scratch instead of using the recommended Matrix package for R 
> that Martin Mӓchler and I wrote).  I have the suspicion that the ideas 
> are the same, it's just the names that are different.
>
> On Tuesday, December 16, 2014 10:03:06 AM UTC-7, Douglas Bates wrote:
>>>
>>> Is anyone working on a Julia package or packages to work with pedigrees, 
>>> in the sense used in animal breeding?  Ana Vazquez and I are developing 
>>> such a package which we plan to call Pedigrees.  Its principle type will be 
>>> Pedigree, which is constructed from two Integer vectors representing the 
>>> 'sire' and 'dam' of each animal.  These names reflect the animal breeding 
>>> orientation of the package.
>>>
>>> Capabilities include sorting a pedigree so that parents occur before 
>>> offspring and deriving the additive relationship matrix and other 
>>> quantities related to it, such as inbreeding coefficients.
>>>
>>> Is using the name Pedigrees going to cause conflicts?
>>>
>>

[julia-users] Figure out how to get pointer to ge

[Taylor Maxwell]

I am moving some of my code to 0.4 and I am having trouble figuring out how 
to get a pointer for beginning of a UInt matrix.  In the past I did:

*julia> **bpt = convert(Ptr{Uint8},b)  #where b it my Uint8 matrix*

*Ptr{Uint8} @0x7fdd243fca70*

in 0.4 I get:

*julia> **bpt = convert(Ptr{UInt8},b)*

*ERROR: MethodError: `convert` has no method matching 
convert(::Type{Ptr{UInt8}}, ::Array{UInt8,1})*

*This may have arisen from a call to the constructor Ptr{UInt8}(...),*

*since type constructors fall back to convert methods.*

Closest candidates are:

  call{T}(::Type{T}, ::Any)

  convert{T<:Union{Int8,UInt8}}(::Type{Ptr{T<:Union{Int8,UInt8}}}, 
*::Cstring*)

  convert{T}(::Type{Ptr{T}}, *::UInt64*)

Unfortunately I am not certain any of these options match what I am looking 
for:

The method used in 0.3 is:
convert{T}(::Type{Ptr{T}},a::Array{T,N}) at pointer.jl:22
which is:
convert{T}(::Type{Ptr{T}}, a::Array{T}) = ccall(:jl_array_ptr, Ptr{T}, 
(Any,), a)
this obviously does not exist in 0.4.

My use case is that I read in a Uint8 matrix from a PLINK .bed file format 
which is a Uint8 matrix that is a dense way to pack in genetic locus 
genotypes every two bits and every column represents a new locus.  With the 
pointer at the beginning of the Uint8 matrix, I can easily work my way down 
a column and extract the genotype calls with bit operations across the 
matrix.  This kind of work is at the edge of my understanding (i.e. 
pointers and bit operations) so despite looking at some of the code for the 
suggested convert methods I am not sure how to get back to what I 
originally use with what is available in 0.4.

In 0.4 I can just make my own new convert call with the old code from 0.3 
by:
import Base.convert
convert{T}(::Type{Ptr{T}}, a::Array{T}) = ccall(:jl_array_ptr, Ptr{T}, 
(Any,), a)
and it get:

*julia> **bpt = convert(Ptr{UInt8},b)*

*Ptr{UInt8} @0x00011253ac00*


But I wanted to know if the was a technical reason for that convert call to 
be removed in 0.4 and if using this "reinstated" convert call will do what 
it originally did in 0.3.

Re: [julia-users] Figure out how to get pointer to ge

[Taylor Maxwell]

Thanks for your response.  It looks like pointer(::Array) gives me the same 
pointer as if I defined the base.convert call from 0.3.

I am not passing it to a ccall.  Below is an example of what i am doing. 
 The bedfreq function gets the genotype and missing counts (4 possibilities 
encoded in every 2 bits) where b is the UInt8 matrix and ns is the number 
of individuals (this helps me figure out how far to go for each locus.  The 
pj variable within the function moves the pointer to the beginning of each 
column (each locus). The ival function used to pull out the genotype or 
missing call for each subject which is then added up in the cc matrix.

I read the garbage collection safety section in documentation and and only 
partially understand it.  Once the b Uint8 Matrix is created, it will not 
be changed or deleted while the pointer is in use.  From your knowledge of 
the the safe and unsafe uses of pointer to an array, is this ok (it worked 
perfectly fine in 0.3).

ival(p::Ptr{UInt8},i::Integer) = i>0 ? 
(unsafe_load(p,(i-1)>>2+1)>>(2(i-1)&0x06))&0x03 : throw(DomainError)

function bedfreq(b::Matrix{UInt8},ns::Integer)
m,n = size(b)
div(ns+3,4) == m || error("ns = $ns should be in [$(4m-3),$(4m)] for 
size(b) = $(size(b))")
cc = zeros(Int, 4, n)
bpt = convert(Ptr{UInt8},b)
for j in 0:(n-1)
pj = bpt + j*m
## the ival function is inlined explicitly
for i in 1:ns cc[1+ival(pj,i),1+j] += 1 end
end
cc
end

Thanks!

Re: [julia-users] Figure out how to get pointer to ge

[Taylor Maxwell]

Thanks, I will try and figure out if I can work this out from your 
suggestions.


On Tuesday, October 27, 2015 at 8:29:07 AM UTC-6, Yichao Yu wrote:
>
> On Tue, Oct 27, 2015 at 10:08 AM, Taylor Maxwell 
> > wrote: 
> > Thanks for your response.  It looks like pointer(::Array) gives me the 
> same 
> > pointer as if I defined the base.convert call from 0.3. 
> > 
> > I am not passing it to a ccall.  Below is an example of what i am doing. 
> > The bedfreq function gets the genotype and missing counts (4 
> possibilities 
> > encoded in every 2 bits) where b is the UInt8 matrix and ns is the 
> number of 
> > individuals (this helps me figure out how far to go for each locus.  The 
> pj 
> > variable within the function moves the pointer to the beginning of each 
> > column (each locus). The ival function used to pull out the genotype or 
> > missing call for each subject which is then added up in the cc matrix. 
>
> I didn't look at the code close enough to rewrite the code myself but 
> it seems like you could just pass the column number to ival instead. 
> You could also construct a sub array for the column. 
>
> This also avod some undefined behavior since the GC/codegen might get 
> smarter and figure out that you are not actually using the array 
> anymore after casting it to a pointer and free it while you are still 
> using the pointer. 
>
> > 
> > I read the garbage collection safety section in documentation and and 
> only 
> > partially understand it.  Once the b Uint8 Matrix is created, it will 
> not be 
> > changed or deleted while the pointer is in use.  From your knowledge of 
> the 
> > the safe and unsafe uses of pointer to an array, is this ok (it worked 
> > perfectly fine in 0.3). 
>
> The separation of cconvert and unsafe_convert does not really matter 
> for this case. But in general, there are cases that the conversion 
> might allocate new buffers for the c function so you want to put that 
> part in cconvert and only do the unsafe take-pointer part in 
> unsafe_convert. I think we remove the convert method for converting to 
> Ptr because you almost never want that convertion to happen 
> automatically. (It's unsafe and you should be explicit if you want it, 
> by calling pointer() (which should also be unnecessary in most cases), 
> for example. 
>
> > 
> > ival(p::Ptr{UInt8},i::Integer) = i>0 ? 
> > (unsafe_load(p,(i-1)>>2+1)>>(2(i-1)&0x06))&0x03 : throw(DomainError) 
> > 
> > function bedfreq(b::Matrix{UInt8},ns::Integer) 
> > m,n = size(b) 
> > div(ns+3,4) == m || error("ns = $ns should be in [$(4m-3),$(4m)] for 
> > size(b) = $(size(b))") 
>
> hint, we have \div operator 
>
> > cc = zeros(Int, 4, n) 
> > bpt = convert(Ptr{UInt8},b) 
> > for j in 0:(n-1) 
> > pj = bpt + j*m 
> > ## the ival function is inlined explicitly 
> > for i in 1:ns cc[1+ival(pj,i),1+j] += 1 end 
> > end 
> > cc 
> > end 
> > 
> > Thanks! 
>

[julia-users] can I load specific lines of a script file into REPL without copy/paste?

[Taylor Maxwell]

I have looked around and haven't been able to find an answer (I am sure it 
is simple).

Say I have a script.jl file with a bunch of commands that I want to work 
with interactively in the REPL but only want to do a few lines at a time 
and I am unable to copy/paste. 

Is there a way to load/import/include a specific set of lines of a 
script.jl file into the REPL without copy/paste or putting those lines into 
a new file? 

I can include("script.jl") to read in the whole file but I would like to do 
it piecemeal if possible.

Thanks

Re: [julia-users] can I load specific lines of a script file into REPL without copy/paste?

[Taylor Maxwell]

Ordinarily that would work but my problem is that I need to do it from the 
command line on a remote workstation (the data files cannot be moved from 
that computer) and I can only have one terminal open in the vnc.  I also 
generally cannot install additional software such at atom (someone else has 
to do it and it would be a giant hassle unless absolutely necessary).

Thanks,
Taylor

On Tuesday, January 5, 2016 at 10:08:47 AM UTC-7, Spencer Russell wrote:
>
> The Atom interface [1] works really well for that sort of thing. You can 
> just highlight the lines you want to run and press ctrl+enter to execute 
> them.
>  
> -s
>  
> [1]: https://github.com/JunoLab/atom-julia-client
>  
>  
> On Tue, Jan 5, 2016, at 11:54 AM, Taylor Maxwell wrote:
>
> I have looked around and haven't been able to find an answer (I am sure it 
> is simple).
>  
> Say I have a script.jl file with a bunch of commands that I want to work 
> with interactively in the REPL but only want to do a few lines at a time 
> and I am unable to copy/paste. 
>  
> Is there a way to load/import/include a specific set of lines of a 
> script.jl file into the REPL without copy/paste or putting those lines into 
> a new file? 
>  
> I can include("script.jl") to read in the whole file but I would like to 
> do it piecemeal if possible.
>  
> Thanks
>
>  
>

Re: [julia-users] can I load specific lines of a script file into REPL without copy/paste?

[Taylor Maxwell]

Thanks, I will see if I can figure it out.

On Tuesday, January 5, 2016 at 11:54:06 AM UTC-7, Tim Holy wrote:
>
> You should be able to write your own method 
>
> include(filename, linerange) 
>
> and specify 3:17 to load line numbers 3-17. 
>
> Best, 
> --Tim 
>
> On Tuesday, January 05, 2016 09:25:40 AM Taylor Maxwell wrote: 
> > Ordinarily that would work but my problem is that I need to do it from 
> the 
> > command line on a remote workstation (the data files cannot be moved 
> from 
> > that computer) and I can only have one terminal open in the vnc.  I also 
> > generally cannot install additional software such at atom (someone else 
> has 
> > to do it and it would be a giant hassle unless absolutely necessary). 
> > 
> > Thanks, 
> > Taylor 
> > 
> > On Tuesday, January 5, 2016 at 10:08:47 AM UTC-7, Spencer Russell wrote: 
> > > The Atom interface [1] works really well for that sort of thing. You 
> can 
> > > just highlight the lines you want to run and press ctrl+enter to 
> execute 
> > > them. 
> > > 
> > > -s 
> > > 
> > > [1]: https://github.com/JunoLab/atom-julia-client 
> > > 
> > > 
> > > On Tue, Jan 5, 2016, at 11:54 AM, Taylor Maxwell wrote: 
> > > 
> > > I have looked around and haven't been able to find an answer (I am 
> sure it 
> > > is simple). 
> > > 
> > > Say I have a script.jl file with a bunch of commands that I want to 
> work 
> > > with interactively in the REPL but only want to do a few lines at a 
> time 
> > > and I am unable to copy/paste. 
> > > 
> > > Is there a way to load/import/include a specific set of lines of a 
> > > script.jl file into the REPL without copy/paste or putting those lines 
> > > into 
> > > a new file? 
> > > 
> > > I can include("script.jl") to read in the whole file but I would like 
> to 
> > > do it piecemeal if possible. 
> > > 
> > > Thanks 
>
>

Re: [julia-users] can I load specific lines of a script file into REPL without copy/paste?

[Taylor Maxwell]

Thank you!! I hadn't gotten very far. I was stuck on looking at the 
include_from_node1 function that include seems to refer to.

On Wednesday, January 6, 2016 at 3:03:15 PM UTC-7, Ethan Anderes wrote:
>
> Taylor: just for your reference, here is an implementation along the lines 
> of what Tim suggested:
>
> function Base.include(path::AbstractString, lines::AbstractVector{Int64})
> lns = open(readlines, path)
> include_string(join(lns[lines]))
> end
>
> Now include("file.jl", 10:20) will evaluate lines 10 through 20. You can 
> also do include("file.jl", [1:2;9:10]) to evaluate lines 1 through 2 and 
> then 9 through 10. 
>
> I’m going to put this is my .juliarc.jl :)
>
> On Tuesday, January 5, 2016 at 1:46:05 PM UTC-8, Taylor Maxwell wrote:
>
> Thanks, I will see if I can figure it out.
>>
>> On Tuesday, January 5, 2016 at 11:54:06 AM UTC-7, Tim Holy wrote:
>>>
>>> You should be able to write your own method 
>>>
>>> include(filename, linerange) 
>>>
>>> and specify 3:17 to load line numbers 3-17. 
>>>
>>> Best, 
>>> --Tim 
>>>
>>> On Tuesday, January 05, 2016 09:25:40 AM Taylor Maxwell wrote: 
>>> > Ordinarily that would work but my problem is that I need to do it from 
>>> the 
>>> > command line on a remote workstation (the data files cannot be moved 
>>> from 
>>> > that computer) and I can only have one terminal open in the vnc.  I 
>>> also 
>>> > generally cannot install additional software such at atom (someone 
>>> else has 
>>> > to do it and it would be a giant hassle unless absolutely necessary). 
>>> > 
>>> > Thanks, 
>>> > Taylor 
>>> > 
>>> > On Tuesday, January 5, 2016 at 10:08:47 AM UTC-7, Spencer Russell 
>>> wrote: 
>>> > > The Atom interface [1] works really well for that sort of thing. You 
>>> can 
>>> > > just highlight the lines you want to run and press ctrl+enter to 
>>> execute 
>>> > > them. 
>>> > > 
>>> > > -s 
>>> > > 
>>> > > [1]: https://github.com/JunoLab/atom-julia-client 
>>> > > 
>>> > > 
>>> > > On Tue, Jan 5, 2016, at 11:54 AM, Taylor Maxwell wrote: 
>>> > > 
>>> > > I have looked around and haven't been able to find an answer (I am 
>>> sure it 
>>> > > is simple). 
>>> > > 
>>> > > Say I have a script.jl file with a bunch of commands that I want to 
>>> work 
>>> > > with interactively in the REPL but only want to do a few lines at a 
>>> time 
>>> > > and I am unable to copy/paste. 
>>> > > 
>>> > > Is there a way to load/import/include a specific set of lines of a 
>>> > > script.jl file into the REPL without copy/paste or putting those 
>>> lines 
>>> > > into 
>>> > > a new file? 
>>> > > 
>>> > > I can include("script.jl") to read in the whole file but I would 
>>> like to 
>>> > > do it piecemeal if possible. 
>>> > > 
>>> > > Thanks 
>>>
>>> ​
>

[julia-users] Re: looking for an equivalent function to "independence_test" in the R coin package

[Taylor Maxwell]

Depending on what parts of their package you are interested in, to 
implement a small subset of their functionality I recently wrote some code 
in both Julia and R (to make sure I had an idea what they were doing to 
replicate their results)  based on the equations in page 258 of:
Hothorn T, Hornik K, Van De Wiel MA, Zeileis A, 2006. A lego system for 
conditional inference. The American Statistician 60:257-263.

The base "independencetest" function I wrote only implements scenarios 
where you have ytrans (one or more transforms of a dependent variable) and 
one set of factors as a matrix of indicator variable (xtrans) although I 
think incorporating more on the "x" side would not be hard.  It does not 
take into account multiple "blocks" and so far I have not needed to look 
into how hard that would be.

On top of this I wrote functions to do the some of the transformations that 
convert y to ytrans as input into the independenttest function (which is 
surprisingly simple).  The transformations I wrote make it possible to 
implement the Kruskal-Wallis test for location, the Fligner-Killeen test 
for scale, and a multisample Lepage test that incorporates those two tests 
jointly (the reason I did this), all with the quadratic test (Cquad).  It 
would be trivial to write the transformations for many of the other 
nonparametric tests. I wrote it to do both a parametric test as well one 
using Monte Carlo. It appears to work well and reasonably fast with the 
little testing I have done.

If you are interested I could send you the code.


On Thursday, February 11, 2016 at 5:27:26 PM UTC-7, Albert Cardona wrote:
>
> The formula operator is there via the GLM package [1], but I can't find 
> the equivalent of the independent_test function [2].
>
> [1] https://github.com/JuliaStats/GLM.jl 
> [2] http://www.inside-r.org/packages/cran/coin/docs/IndependenceTest
>
>
> Thanks in advance for any pointers.
>
> Albert
>